A fast method for estimating epicenter location based on strong earthquake acceleration records

Abstract

The present invention discloses a method for quickly estimating the epicenter position based on strong vibration acceleration records, which belongs to the field of geophysical technology and is used for estimating the epicenter position of an earthquake, including extracting the P wave arrival time of the earthquake, constructing the observation equation of the P wave arrival time of the station, preliminarily selecting quasi-observations and preliminarily calculating the true error estimate, using the random sample consensus algorithm RANSAC to fit the eigenvector and identify the inner and outer points, reselecting the inner points as quasi-observations and recalculating the true error estimate, and outputting earthquake elements. The present invention introduces a more efficient and rapid quasi-observation selection scheme to optimize the accuracy of the gross error quasi-observation verification method, and realizes real-time and rapid estimation of the epicenter position of a large earthquake based on strong vibration observation means; realizes more accurate and reliable determination of the epicenter position of an earthquake, and can provide technical support for rapid earthquake warning and post-earthquake emergency response.

Description

Technical Field

[0001] The invention discloses a method for quickly estimating the epicenter position based on strong vibration acceleration records, and belongs to the technical field of geophysics. Background Art

[0002] Among the methods for calculating the epicenter position using strong earthquake records, linear positioning methods represented by the traditional Geiger method and nonlinear positioning methods such as grid search methods have been widely used. However, the Geiger method relies on P-wave arrival information. The sensitivity of the P-wave arrival extraction algorithm to the triggering of observation stations at different epicentral distances is different, and the accuracy of the extracted arrival time is also different. At the same time, it may be affected by the gross error of the P-wave arrival time, thereby reducing the accuracy of the epicenter position calculation. Generally speaking, the observations containing gross errors are less than 1% to 10% of the total data volume. Therefore, most of the observation data can be considered as normal observations, and these normal observations are called quasi-accurate observations QAO. The traditional quasi-accurate verification method QUAD locates the gross errors by adding the condition that the norm of the true error estimate of the quasi-accurate observation is extremely small, constructs a rank-deficient equation group between the true error and the observation, and solves the true error estimate. The gross errors existing in a certain type of observations are accurately located from the true error estimate. However, the traditional quasi-observation selection method includes two steps: preliminary selection and reselection. Multiple true error indicators are introduced in the reselection process, and the category of the observation value needs to be judged, which is cumbersome. It is necessary to introduce a more efficient and rapid quasi-observation selection scheme to optimize the accuracy of the gross error quasi-observation verification method. In addition, the robust estimation method can be constructed in the form of an IGG3 equivalent weight function. Based on iterative calculation, the influence of the gross error of the observation value on the parameter estimation result is suppressed. However, if most of the observation values ​​are contaminated by gross errors, the accuracy of the epicenter estimation result is poor. Summary of the invention

[0003] The purpose of the present invention is to provide a method for quickly estimating the epicenter position based on strong vibration acceleration records, so as to solve the problem in the prior art that the traditional Geiger linear positioning method is affected by the gross error of the P wave arrival time, resulting in inaccurate epicenter position, and the IGG3 anti-error estimation method has poor accuracy of epicenter estimation results when most observation values ​​are contaminated by gross errors.

[0004] A method for quickly estimating the epicenter position based on strong earthquake acceleration records, comprising:

[0005] The first step is to obtain the acceleration records of each station, use the short-time average value STA, long-time average value LTA and Akaike information criterion AIC to extract the P wave arrival time of the earthquake, and determine the number of stations triggered by the P wave arrival. If it is greater than or equal to 4, proceed to the next step, otherwise return to extract the P wave arrival time of the earthquake;

[0006] The second step is to obtain the P-wave arrival time of the first-touch station as the initial value of the earthquake occurrence time, use the centroid position of all triggering stations as the first initial value of the P-wave arrival observation equation, use the theoretical P-wave propagation speed as the second initial value of the P-wave arrival observation equation, and construct the station P-wave arrival observation equation;

[0007] The third step is to preliminarily select the quasi-observation And initially calculate the true error estimate , solve the estimated value of the seismic elements after correcting the gross errors, according to Constructing feature vectors , using the random sample consensus algorithm RANSAC fitting And mark the internal and external points, check the internal points as the quasi-observation points And recalculate the true error estimate , judge the change of the RANSAC fitting internal point, perform test statistics, and if it does not exceed the test statistics threshold, construct the robust equivalent weight to estimate the seismic elements , is the sequence number of acceleration records. If it exceeds the test statistic threshold, it will be returned for reselection and recalculation. and Is the absolute value of the difference less than or equal to the initial value of the proposed observation If yes, then output earthquake elements, otherwise execute and , return to conduct primary elections and preliminary calculations.

[0008] Extracting the arrival time of the P wave of an earthquake includes introducing a characteristic function and calculating the ratio of STA and LTA to preliminarily extract the arrival time of the P wave phase in the seismic signal:

[0009] ;

[0010] ;

[0011] In the formula, It's an earthquake signal The characteristic function of yes The recorded value of vertical acceleration at the moment, yes STA at the moment, yes LTA at the moment, is the number of record points contained in the short window, is the number of record points contained in the long window, and are the time values ​​obtained by LTA and STA respectively. When the trigger threshold is exceeded, it indicates that an abnormal signal has appeared and an earthquake event is determined to have occurred.

[0012] AIC is used to further extract the arrival time of the P-wave phase in the seismic signal:

[0013] ;

[0014] In the formula, is the AIC function expression, represents the variance calculation function, It means An array of acceleration values ​​for the interval, is the number of acceleration recording points contained in the entire recording window. The time corresponding to the minimum point of the AIC function is the arrival time of the P-wave phase in the seismic signal.

[0015] After extracting the arrival time of the earthquake's P wave, the spatial distance calculation formula from each station to the epicenter is:

[0016] ;

[0017] ;

[0018] In the formula, represents the number of stations, for The epicentral distance of the station, for The location of the station, is the initial position of the epicenter in the Earth-centered Earth-fixed coordinate system, for The P-wave arrival time of the trigger at the station is equivalent to the earthquake arrival time function, For the moment of earthquake, is the velocity of seismic wave propagation;

[0019] Eliminate the earthquake occurrence time by using the difference method between stations :

[0020] ;

[0021] In the formula, is the total number of triggered stations, and the error equation is obtained by linearization: :

[0022] ;

[0023] ;

[0024] ;

[0025] ;

[0026] In the formula, is the design matrix, is the observation vector, To correct the number, , and is the correction for the epicenter position, is the correction for the seismic P-wave velocity.

[0027] After judging that the number of stations triggered by the P wave is greater than or equal to 4, the centroid position of the station is calculated as , using the least squares principle to solve :

[0028] ;

[0029] ;

[0030] In the formula, is the weight matrix extracted from seismic waves at each station, is the identity matrix, is the order of the identity matrix, we get The corrected epicenter position and P-wave velocity are then obtained, and the corrected epicenter position is converted from the Earth-centered Earth-fixed coordinate system to the geodetic coordinate system to obtain the geodetic longitude and latitude and geodetic height of the epicenter position.

[0031] The observation equation of the station P wave arrival time is:

[0032] ;

[0033] In the formula, for dimensional coefficient matrix, for The true value vector of the seismic element to be estimated. yes The observed value of the P wave arrival time, yes True error vector, construct adjustment factor matrix :

[0034] ;

[0035] In the formula, for dimensional observation weight matrix, the relationship between the true error and the observation value is:

[0036] ;

[0037] right Make an estimate to obtain the estimated value of the gross error and its position in the observation sequence. Make an adjustment based on the measurement of the equipment to obtain the absolute value of the residual. A pseudo observation, , No. The quasi-observation corresponds to the observation vector , the corresponding true error is , the corresponding weight is , the corresponding coefficient matrix is ,

[0038] Remaining There are gross errors in the observations, The true error corresponding to the observation value with gross error is , introduce additional norm minimum constraint:

[0039] ;

[0040] Find a definite solution for the true error estimate:

[0041] ;

[0042] ;

[0043] ;

[0044] ;

[0045] In the formula, is the residual vector of the fitting, and yes The two components of For observations with gross errors, The observation vector corresponding to the quasi-observation, , are two parameter matrices.

[0046] Solving the estimation of seismic elements after gross error correction includes obtaining the true error estimation:

[0047] ;

[0048] When the observed value contains gross errors, the true error estimation has the characteristics of grouping, and the values ​​greater than the threshold are The observations are judged to contain gross errors. Assume that the observation sequence is A rough error, indivual -dimensional unit vector:

[0049] ;

[0050] ;

[0051] In the formula, Corresponding to An observation with gross error, Middle The components are 1 and the rest are 0, and the gross error is expressed as , rewrite the observation equation as:

[0052] ;

[0053] ;

[0054] In the formula, for dimensional coefficient matrix, After separating gross errors , obtained by the least squares criterion Valuation :

[0055] ;

[0056] Estimation of seismic elements after gross error correction for:

[0057] .

[0058] In RANSAC, it is assumed that the data set contains normal values ​​and abnormal values. The normal values ​​are recorded as inliers and the abnormal values ​​are recorded as outliers. The objective function of RANSAC is:

[0059] ;

[0060] ;

[0061] In the formula, is the best fitting function, is the indicator function, when the data point for The value is 1 when it is an interior point, otherwise it is 0. Represents the residual function, expressed as the Euclidean distance from a point to a straight line. For the data set;

[0062] The identification of internal and external points includes, for a given and mathematical models , the minimum sample set MSS for each sampling is recorded as , the number of samples in MSS is ,exist Randomly select 1 , and according to The sample points in the model are used to calculate the model parameters and fit the mathematical model ; For other sample points in the data set, calculate the residual between the sample and the fitted model, and set the residual threshold to If the difference is less than the threshold, it is an internal point. If the difference is greater than the threshold, it is an external point. The number of internal points is recorded.

[0063] Repeat the internal and external point identification, calculate the number of internal points and model parameters in this cycle, if the number of internal points this time is greater than the previous number of internal points, save the number of internal points and model parameters calculated this time; otherwise, keep the previous calculation parameters until the maximum number of iterations is met :

[0064] ;

[0065] In the formula, for After iterations, the probability that all points are internal points in at least one sampling is set to 99%. Get the probability of the correct model for each iteration, represents the probability of failure of a single iteration;

[0066] Select the set of model parameters with the largest number of inliers as the optimal model for output, and calculate the inlier rate of the model at the same time :

[0067] ;

[0068] In the formula, is the number of external points, is the number of internal points;

[0069] The preliminary and re-selection methods are consistent with the proposed standard inspection method.

[0070] When checking, confirm , and select As the initial number of quasi-observations, calculate and of The value of the moment , check the quasi-observation, calculate :

[0071] ;

[0072] Using RANSAC Perform linear model fitting, and the fitting result is ,Will The internal and external point identifications are used as the group identification of the true error, and the internal point set corresponding to RANSAC is selected The corresponding measurement component is the quasi-observation;

[0073] Recalculate based on the selected quasi-observation and , RANSAC is used again to determine the true error clustering characteristics. If the clustering identification of the internal and external points changes significantly, the quasi-observation is adjusted, otherwise proceed to the next step;

[0074] For non-quasi-accurate observations, the IGG3 equivalent weight function is constructed, and combined with the selected quasi-accurate observations, the estimated parameters are calculated using the quasi-accurate test. If the difference between the estimated parameters of the two iterations is less than the threshold or reaches the maximum number of iterations, the filtering ends at the current moment; otherwise, the number of iterations is increased by 1, and ;

[0075] like The current moment estimation ends; otherwise, re-execute the check.

[0076] The final calculated earthquake occurrence time is:

[0077] .

[0078] Compared with the prior art, the present invention has the following beneficial effects: the present invention introduces a more efficient and rapid quasi-observation selection scheme to optimize the accuracy of the gross error quasi-observation verification method, and realizes real-time rapid estimation of the epicenter of a large earthquake based on strong vibration observation means; it realizes a more accurate and reliable determination of the epicenter of an earthquake, and can provide technical support for rapid earthquake warning and post-earthquake emergency response. BRIEF DESCRIPTION OF THE DRAWINGS

[0079] Figure 1 It is a technical flow chart of the present invention;

[0080] Figure 2 Schematic diagram of X- and Y-coordinates of earthquake events and stations distribution;

[0081] Figure 3 Schematic diagram of the Y and Z coordinates of earthquake events and station distribution;

[0082] Figure 4 Schematic diagram of X- and Z-coordinates of earthquake events and station distribution;

[0083] Figure 5 This is a schematic diagram of the P wave arrival time extraction results of the first earthquake event;

[0084] Figure 6 This is a schematic diagram of the P wave arrival time extraction results of the second earthquake event;

[0085] Figure 7 This is a schematic diagram of the P-wave arrival time extraction results for the third earthquake event;

[0086] Figure 8 Schematic diagram of the P-wave arrival time extraction results for the fourth earthquake event. DETAILED DESCRIPTION

[0087] In order to make the purpose, technical solution and advantages of the present invention clearer, the technical solution of the present invention is described clearly and completely below. Obviously, the described embodiments are part of the embodiments of the present invention, but not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0088] The technical flow chart of the present invention is as follows: Figure 1 As shown, a method for quickly estimating the epicenter position based on strong vibration acceleration records includes:

[0089] The first step is to obtain the acceleration records of each station, use the short-time average value STA, long-time average value LTA and Akaike information criterion AIC to extract the P wave arrival time of the earthquake, and determine the number of stations triggered by the P wave arrival. If it is greater than or equal to 4, proceed to the next step, otherwise return to extract the P wave arrival time of the earthquake;

[0090] The second step is to obtain the P-wave arrival time of the first-touch station as the initial value of the earthquake occurrence time, use the centroid position of all triggering stations as the first initial value of the P-wave arrival observation equation, use the theoretical P-wave propagation speed as the second initial value of the P-wave arrival observation equation, and construct the station P-wave arrival observation equation;

[0091] The third step is to preliminarily select the quasi-observation And initially calculate the true error estimate , solve the estimated value of the seismic elements after correcting the gross errors, according to Constructing feature vectors , using the random sample consensus algorithm RANSAC fitting And mark the internal and external points, check the internal points as the quasi-observation points And recalculate the true error estimate , judge the change of the RANSAC fitting internal point, perform test statistics, and if it does not exceed the test statistics threshold, construct the robust equivalent weight to estimate the seismic elements , is the sequence number of acceleration records. If it exceeds the test statistic threshold, it will be returned for reselection and recalculation. and Is the absolute value of the difference less than or equal to the initial value of the proposed observation If yes, then output earthquake elements, otherwise execute and , return to conduct primary elections and preliminary calculations.

[0092] Extracting the arrival time of the P wave of an earthquake includes introducing a characteristic function and calculating the ratio of STA and LTA to preliminarily extract the arrival time of the P wave phase in the seismic signal:

[0093] ;

[0094] ;

[0095] In the formula, It's an earthquake signal The characteristic function of yes The recorded value of vertical acceleration at the moment, yes STA at the moment, yes LTA at the moment, is the number of record points contained in the short window, is the number of record points contained in the long window, and are the time values ​​obtained by LTA and STA respectively. When the trigger threshold is exceeded, it indicates that an abnormal signal has appeared and an earthquake event is determined to have occurred.

[0096] AIC is used to further extract the arrival time of the P-wave phase in the seismic signal:

[0097] ;

[0098] In the formula, is the AIC function expression, represents the variance calculation function, It means An array of acceleration values ​​for the interval, is the number of acceleration recording points contained in the entire recording window. The time corresponding to the minimum point of the AIC function is the arrival time of the P-wave phase in the seismic signal.

[0099] After extracting the arrival time of the earthquake's P wave, the spatial distance calculation formula from each station to the epicenter is:

[0100] ;

[0101] ;

[0102] In the formula, represents the number of stations, for The epicentral distance of the station, for The location of the station, is the initial position of the epicenter in the Earth-centered Earth-fixed coordinate system, for The P-wave arrival time of the trigger at the station is equivalent to the earthquake arrival time function, For the moment of earthquake, is the velocity of seismic wave propagation;

[0103] Eliminate the earthquake occurrence time by using the difference method between stations :

[0104] ;

[0105] In the formula, is the total number of triggered stations, and the error equation is obtained by linearization: :

[0106] ;

[0107] ;

[0108] ;

[0109] ;

[0110] In the formula, is the design matrix, is the observation vector, To correct the number, , and is the correction for the epicenter position, is the correction for the seismic P-wave velocity.

[0111] After judging that the number of stations triggered by the P wave is greater than or equal to 4, the centroid position of the station is calculated as , using the least squares principle to solve :

[0112] ;

[0113] ;

[0114] In the formula, is the weight matrix extracted from seismic waves at each station, is the identity matrix, is the order of the identity matrix, we get The corrected epicenter position and P-wave velocity are then obtained, and the corrected epicenter position is converted from the Earth-centered Earth-fixed coordinate system to the geodetic coordinate system to obtain the geodetic longitude and latitude and geodetic height of the epicenter position.

[0115] The observation equation of the station P wave arrival time is:

[0116] ;

[0117] In the formula, for dimensional coefficient matrix, for The true value vector of the seismic element to be estimated. yes The observed value of the P wave arrival time, yes True error vector, construct adjustment factor matrix :

[0118] ;

[0119] In the formula, for dimensional observation weight matrix, the relationship between the true error and the observation value is:

[0120] ;

[0121] right Make an estimate to obtain the estimated value of the gross error and its position in the observation sequence. Make an adjustment based on the measurement of the equipment to obtain the absolute value of the residual. A pseudo observation, , No. The quasi-observation corresponds to the observation vector , the corresponding true error is , the corresponding weight is , the corresponding coefficient matrix is ,

[0122] Remaining There are gross errors in the observations, The true error corresponding to the observation value with gross error is , introduce additional norm minimum constraint:

[0123] ;

[0124] Find a definite solution for the true error estimate:

[0125] ;

[0126] ;

[0127] ;

[0128] ;

[0129] In the formula, is the residual vector of the fitting, and yes The two components of For observations with gross errors, The observation vector corresponding to the quasi-observation, , are two parameter matrices.

[0130] Solving the estimation of seismic elements after gross error correction includes obtaining the true error estimation:

[0131] ;

[0132] When the observed value contains gross errors, the true error estimation has the characteristics of grouping, and the values ​​greater than the threshold are The observations are judged to contain gross errors. Assume that the observation sequence is A rough error, indivual -dimensional unit vector:

[0133] ;

[0134] ;

[0135] In the formula, Corresponding to An observation with gross error, Middle The components are 1 and the rest are 0, and the gross error is expressed as , rewrite the observation equation as:

[0136] ;

[0137] ;

[0138] In the formula, for dimensional coefficient matrix, After separating gross errors , obtained by the least squares criterion Valuation :

[0139] ;

[0140] Estimation of seismic elements after gross error correction for:

[0141] .

[0142] In RANSAC, it is assumed that the data set contains normal values ​​and abnormal values. The normal values ​​are recorded as inliers and the abnormal values ​​are recorded as outliers. The objective function of RANSAC is:

[0143] ;

[0144] ;

[0145] In the formula, is the best fitting function, is the indicator function, when the data point for The value is 1 when it is an interior point, otherwise it is 0. Represents the residual function, expressed as the Euclidean distance from a point to a straight line. For the data set;

[0146] The identification of internal and external points includes, for a given and mathematical models , the minimum sample set MSS for each sampling is recorded as , the number of samples in MSS is ,exist Randomly select 1 , and according to The sample points in the model are used to calculate the model parameters and fit the mathematical model ; For other sample points in the data set, calculate the residual between the sample and the fitted model, and set the residual threshold to If the difference is less than the threshold, it is an internal point. If the difference is greater than the threshold, it is an external point. The number of internal points is recorded.

[0147] Repeat the internal and external point identification, calculate the number of internal points and model parameters in this cycle, if the number of internal points this time is greater than the previous number of internal points, save the number of internal points and model parameters calculated this time; otherwise, keep the previous calculation parameters until the maximum number of iterations is met :

[0148] ;

[0149] In the formula, for After iterations, the probability that all points are internal points in at least one sampling is set to 99%. Get the probability of the correct model for each iteration, represents the probability of failure of a single iteration;

[0150] Select the set of model parameters with the largest number of inliers as the optimal model for output, and calculate the inlier rate of the model at the same time :

[0151] ;

[0152] In the formula, is the number of external points, is the number of internal points;

[0153] The preliminary and re-selection methods are consistent with the proposed standard inspection method.

[0154] When checking, confirm , and select As the initial number of quasi-observations, calculate and of The value of the moment , check the quasi-observation, calculate :

[0155] ;

[0156] Using RANSAC Perform linear model fitting, and the fitting result is ,Will The internal and external point identifications are used as the group identification of the true error, and the internal point set corresponding to RANSAC is selected The corresponding measurement component is the quasi-observation;

[0157] Recalculate based on the selected quasi-observation and , RANSAC is used again to determine the true error clustering characteristics. If the clustering identification of the internal and external points changes significantly, the quasi-observation is adjusted, otherwise proceed to the next step;

[0158] For non-quasi-accurate observations, the IGG3 equivalent weight function is constructed, and combined with the selected quasi-accurate observations, the estimated parameters are calculated using the quasi-accurate test. If the difference between the estimated parameters of the two iterations is less than the threshold or reaches the maximum number of iterations, the filtering ends at the current moment; otherwise, the number of iterations is increased by 1, and ;

[0159] like The current moment estimation ends; otherwise, re-execute the check.

[0160] The final calculated earthquake occurrence time is:

[0161] .

[0162] The embodiment of the present invention constructs a strong vibration acceleration data set of 23 earthquake events with magnitudes between 6.0 and 9.0 worldwide, and unifies the time history direction, recording dimension and file format of the acceleration record. In view of the fact that most earthquakes occur in the shallow crust, and the destructive power of earthquakes decreases rapidly with the increase of the focal depth. Therefore, the epicenter estimation method is mainly applied to the positioning of shallow source earthquakes with strong destructiveness. The initial depth of the earthquake is set to 10km, and its uncertainty index is set to 10km. For different regions, different initial depth prior information can be given according to the actual situation of historical earthquakes in the region. In order to test the effectiveness of the gross error quasi-verification filtering algorithm, synthetic seismic event data is used for testing. Using an area of ​​100km×100km, 50 randomly uniformly distributed stations and 16 randomly distributed seismic events are set, and the focal depth of each event is within the range of 5km to 15km. It is assumed that the seismic wave velocity model in the selected rupture area obeys one-dimensional propagation and varies with depth. The initial depth of the earthquake event is set at 10km, with an uncertainty of 5km, and the values ​​of the other initial earthquake elements follow the above simulation method. Each time, 1% to 10% of the P-wave arrival error observations are added in a random sample mode in 50 stations. These gross error observations caused by false triggering or misjudgment of the P-wave arrival are subject to a normal distribution with a variance of 1s. The traditional least squares Geiger method, the epicenter estimation method based on IGG3 robust estimation, the traditional gross error quasi-precision verification method and the new method are used to solve the earthquake elements of 16 synthetic earthquake cases. The epicenter results calculated by the epicenter estimation method adopted in the technical solution of the present invention have an average accuracy improvement rate of 11%, 7% and 5% compared with the Geiger least squares, IGG3 robust estimation and traditional quasi-precision verification methods, respectively. The accuracy of obtaining the earthquake occurrence time is improved by 41%, 23% and 2% compared with the least squares, IGG3 robust estimation and traditional quasi-precision verification methods, respectively.

[0163] In order to verify the accuracy of the positioning results of the proposed method, the empirical P-wave propagation velocity and the earthquake epicenter estimation results and the earthquake occurrence time are used as reference values ​​in the retrospective analysis of real earthquake cases. The P-wave arrival time, epicenter position and earthquake occurrence time of each station are calculated using the four earthquake events collected in the constructed fusion co-seismic deformation dataset. The epicenter results calculated by the epicenter estimation method adopted in this technical solution have an average accuracy improvement rate of 67%, 56% and 27% compared with the Geiger least squares, IGG3 robust estimation and traditional quasi-precision verification methods, respectively. The accuracy of the earthquake occurrence time acquisition is improved by 68%, 67% and 58% compared with the least squares, IGG3 robust estimation and traditional quasi-precision verification methods, respectively. The influence of the gross error in the P-wave phase arrival time identification on the epicenter location results is effectively reduced, and the accuracy of the epicenter location results is improved.

[0164] For the convenience of representation, the traditional least squares Geiger method, the epicenter estimation method based on IGG3 robust estimation, the traditional gross error quasi-accuracy verification method, and the new method are abbreviated as LS, IGG3, Quad, and RANSAC, respectively. The residual sequences obtained by the four schemes are plotted to preliminarily explore the ability of the quasi-accuracy verification method to handle gross errors. From the residual sequences of events 1, 5, 11, 14, 15, and 16, it can be seen that the residuals of the IGG3 and LS schemes still contain the influence of the simulated arrival time gross errors, while the Quad and RANSAC methods can accurately locate the location of the gross errors of each station. The residuals of the non-quasi-accuracy station observations repaired by the quasi-accuracy verification no longer have the characteristics of step and grouping. The residuals of most events are distributed within ±1m, and the peak values ​​of a few residual values ​​containing gross errors exceed ±2m, with the maximum amplitude of about ±4m.

[0165] Table 1 shows the average accuracy of earthquake element estimation of each scheme in all 16 random synthetic events. It can be seen that the positioning results of the traditional least squares method are slightly lower than the accuracy of the other three schemes, mainly because the P-wave arrival time observations of the stations containing gross errors are not eliminated or down-weighted. The IGG3 method uses a three-segment equal weight function to unweight and down-weight the observations containing gross errors and suspected gross errors. The two quasi-accuracy verification methods given by the Quad and RANSAC methods can effectively locate and repair gross errors. The RANSAC method reduces the initial selection principle of the traditional Quad method, and performs anti-error equivalent weighting processing on only non-quasi-accurate observations, which can effectively avoid the possible down-weighting of quasi-accurate observations by the IGG3 method. Its epicenter extraction accuracy in the X, Y and Z directions is 0.2432km, 0.3655km and 0.9378km. For the earthquake onset time indicator, most of the error results are distributed within 2s, and RANSAC has the smallest deviation, which is only 1.08s.

[0166] Table 1. Estimation accuracy of earthquake elements in each scheme in the simulation experiment

[0167] ;

[0168] The STA / LTA and STA / LTA+AIC methods were used to extract the P-wave arrival time of each station for the four earthquake events. Before processing the P-wave arrival time, the original acceleration records were corrected for zero bias, that is, all acceleration records of the four events were subtracted from the average acceleration of each earthquake 5 seconds before the event. The short window of the STP / LTP+AIC method used was set to 1s, the long window was set to 10s, the threshold for picking the P-wave arrival time was set to 4, and the P-wave arrival time was roughly picked. , set the AIC method interval to Under the assumed uniform P-wave propagation velocity model, the processing results of the PphasePicker software and the theoretical P-wave arrival time observation values ​​calculated based on the epicenter distance and the theoretical P-wave velocity are used as references to give the P-wave arrival time results picked up by the STA / LTA+AIC method at each station. In order to determine the number of stations where the P-wave is correctly identified, the identification accuracy is defined as:

[0169] ;

[0170] In the formula, To identify the correct number of stations, is the number of triggered stations, and Ratio is the recognition accuracy.

[0171] All P-wave arrival time observations are arranged in ascending order according to the epicentral distance. All P-wave arrival time observations are converted from moment observations to travel time observations starting from the earthquake occurrence time. The zero value of the ordinate indicates the earthquake occurrence time. The STA / LTA and STA / LTA+AIC methods have poor filtering effects on the noise caused by site effects at different stations, which causes false triggering due to noise advance or lag beyond the empirical threshold, which is quite different from the extraction results of PphasePicker and theoretical P-wave arrival time. With the increase of the epicentral distance, the P-wave energy decays faster, and the P-wave travel time deviation misjudged by the fixed empirical threshold gradually prolongs.

[0172] The Ratio values ​​of the P wave arrival time extraction results of all earthquake cases are shown in Table 2. The Ratio value of the P wave arrival time extraction of the No. 1 earthquake is 0.71s, the Ratio of No. 4 is the highest at 0.97s, the P wave of No. 2 contains obvious gross errors, and its Ratio value is the lowest at 0.48s. Although there are fewer available stations for the No. 3 event, the availability of the P wave arrival time observation value is high, and the Ratio value is 0.93s. All four earthquake cases involved in the solution have inaccurate P wave arrival times. These outliers include some gross errors caused by site effects and misjudged observations suspected of gross errors.

[0173] Table 2. P wave automatic picking results (unit: seconds)

[0174] ;

[0175] The earthquake events and station distribution of the present invention are three-dimensional data, wherein the schematic diagram of the X coordinate and the Y coordinate is shown in FIG2 , and the schematic diagram of the Y coordinate and the Z coordinate is shown in FIG2 . Figure 3 As shown, the schematic diagram of X coordinate and Z coordinate is as follows Figure 4 As shown in Figure 2. The P-wave arrival time extracted by the STA / LTA+AIC method was used as the observation value, and the epicenter positions of the four events were estimated using the LS method, IGG3 method, QUAD method, and RANSAC method. The schematic diagram of the P-wave arrival time extraction result of the first earthquake event is shown in Figure 2. Figure 5 As shown in the figure, the P wave arrival time extraction result diagram of the second earthquake event is shown in the figure. Figure 6 As shown in the figure, the P wave arrival time extraction result diagram of the third earthquake event is shown in the figure. Figure 7 The schematic diagram of the extraction result of the P wave arrival time of the fourth earthquake event is shown in Figure 8 As shown in the figure. Among the four events, the epicenter deviation of No. 2 is the largest. All the strong earthquake stations of this event are located on the same side of the main shock, and the P-wave arrival time observation values ​​between adjacent stations fluctuate greatly. Although there are relatively few stations for No. 3 earthquake, the epicenter deviation of the earthquake is significantly smaller than that of No. 2. This is because the stations of No. 3 earthquake are distributed in the land area northeast of the coastline and are relatively evenly concentrated. Therefore, the uneven distribution of stations will lead to a large deviation in the epicenter estimation results. Although the stations of No. 1 are distributed on the same side of the main shock, the estimated epicenter position of No. 1 earthquake is more accurate. This is because the local strong earthquake stations are densely distributed, and the 100 triggered stations used provide more redundant observations for the least squares method. The deviation of No. 4 earthquake is the smallest. This is because the P-wave arrival time observation values ​​of this earthquake have fewer gross errors. Therefore, after reducing the influence of gross errors, a larger number of triggered stations will improve the accuracy of epicenter positioning.

[0176] In order to quantitatively describe the results of epicenter estimation, the epicenter position published by USGS was used as the reference value, and the estimated epicenter position was subtracted from it. The epicenter deviations of four earthquake events obtained by four methods were statistically analyzed. As shown in Table 3, due to the site effect, noise influence and different magnitudes of different earthquakes, the deviations of the epicenter position estimation vary greatly. The average epicenter deviation was used as the indicator for the final accuracy evaluation. The average epicenter deviation of the four earthquakes calculated by the LS method was 32.50 km. The epicenter results obtained by the IGG3 method were slightly better than those of the LS method, with an average epicenter deviation of 24.71 km. The epicenter accuracy obtained by the traditional Quad method and RANSAC method using gross error quasi-calibration was significantly improved compared with the IGG3 and LS methods, with average epicenter deviations of 14.75 km and 10.64 km. The results show that the average accuracy improvement rates of the epicenter results calculated by the new method compared with the traditional LS, IGG3 and Quad methods are 67%, 56% and 27%, respectively.

[0177] Table 3. Statistics of epicenter estimation accuracy (unit: kilometers)

[0178] ;

[0179] The feasibility of the new method can be further verified by using the LS method and the proposed improved method to obtain the main shock time. Table 4 gives the accuracy statistics of the earthquake time of the four earthquake events. Compared with the theoretical earthquake time published by USGS, the earthquake time obtained by the LS method in the four earthquake events was 5.41s, and the average deviation of the IGG3 method was 4.11s. The average deviation of the earthquake time obtained by the Quad method was 2.45s, and the minimum average earthquake time deviation of the RANSAC method was 1.75s. The results show that the earthquake time acquisition accuracy of the earthquake element estimation method based on random sample consistency and gross error quasi-precision verification is improved by 68%, 67% and 58% compared with the least squares, IGG3 robust estimation and traditional quasi-precision verification methods, respectively.

[0180] Table 4. Statistics of average earthquake occurrence time accuracy (unit: seconds)

[0181] ;

[0182] The epicenter results calculated by the epicenter estimation method adopted in the technical solution of the present invention have an average accuracy improvement rate of 67%, 56% and 27% compared with Geiger least squares, IGG3 robust estimation and traditional quasi-precision calibration methods. The accuracy of obtaining the time of earthquake occurrence has been improved by 68%, 67% and 58% respectively compared with the least squares, IGG3 robust estimation and traditional quasi-precision calibration methods. The influence of gross errors in the identification of the P-wave phase arrival time on the epicenter positioning results is effectively reduced, and the accuracy of the epicenter positioning results is improved. The application effect is good for earthquakes with dense and uniform distribution of stations near the epicenter, and the epicenter deviation is within 2km. However, the epicenter estimation results of earthquakes with sparse and uneven stations still deviate from the reference value by more than 29km. This is because the four earthquake cases used are distributed in different regions of the world, and the arrival accuracy of each earthquake case station under different sites, piers and environmental factors will be slightly different, and the contribution rate of the P-wave phase arrival observation value of each earthquake case station to the estimation of the epicenter position at different epicenter distances is also slightly different, which reduces the applicability of the new method. In addition, this technical solution assumes that the propagation speed of seismic waves in all directions is equal, but the actual propagation speed varies due to different geological structures, so the epicenter location contains errors caused by the uniform propagation assumption. In summary, in practical applications, the anisotropy of propagation speed must also be considered, and the trigger station closer to the epicenter is used to accurately estimate the epicenter position, and the stations farther from the epicenter are eliminated or downgraded.

[0183] The above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit the same. Although the present invention has been described in detail with reference to the aforementioned embodiments, a person skilled in the art should understand that the technical solutions described in the aforementioned embodiments may still be modified, or some or all of the technical features may be replaced by equivalents, and these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for quickly estimating the epicenter position based on strong vibration acceleration records, characterized in that: include: The first step is to obtain the acceleration records of each station, use the short-time average value STA, long-time average value LTA and Akaike information criterion AIC to extract the P wave arrival time of the earthquake, and determine the number of stations triggered by the P wave arrival. If it is greater than or equal to 4, proceed to the next step, otherwise return to extract the P wave arrival time of the earthquake; The second step is to obtain the P-wave arrival time of the first-touch station as the initial value of the earthquake occurrence time, use the centroid position of all triggering stations as the first initial value of the P-wave arrival observation equation, use the theoretical P-wave propagation speed as the second initial value of the P-wave arrival observation equation, and construct the station P-wave arrival observation equation; The third step is to preliminarily select the quasi-observation And initially calculate the true error estimate , solve the estimated value of the seismic elements after correcting the gross errors, according to Constructing feature vectors , using the random sample consensus algorithm RANSAC fitting And mark the internal and external points, check the internal points as the quasi-observation points And recalculate the true error estimate , judge the change of the RANSAC fitting internal point, perform test statistics, and if it does not exceed the test statistics threshold, construct the robust equivalent weight to estimate the seismic elements , is the sequence number of acceleration records. If it exceeds the test statistic threshold, it will be returned for reselection and recalculation. and Is the absolute value of the difference less than or equal to the initial value of the proposed observation If yes, then output earthquake elements, otherwise execute and , return to conduct primary elections and preliminary calculations.

2. A method for rapid estimation of epicenter position based on strong vibration acceleration records according to claim 1, characterized in that: Extracting the arrival time of the P wave of an earthquake includes introducing a characteristic function and calculating the ratio of STA and LTA to preliminarily extract the arrival time of the P wave phase in the seismic signal: ； ； In the formula, It's an earthquake signal The characteristic function of yes The recorded value of vertical acceleration at the moment, yes STA at the moment, yes LTA at the moment, is the number of record points contained in the short window, is the number of record points contained in the long window, and are the time values ​​obtained by LTA and STA respectively. When the trigger threshold is exceeded, it indicates that an abnormal signal has appeared and an earthquake event is determined to have occurred.

3. The method for rapid estimation of epicenter position based on strong vibration acceleration records according to claim 1 is characterized in that: AIC is used to further extract the arrival time of the P-wave phase in the seismic signal: ； In the formula, is the AIC function expression, represents the variance calculation function, It means An array of acceleration values ​​for the interval, is the number of acceleration recording points contained in the entire recording window. The time corresponding to the minimum point of the AIC function is the arrival time of the P-wave phase in the seismic signal.

4. The method for rapid estimation of epicenter position based on strong vibration acceleration records according to claim 1 is characterized in that: After extracting the arrival time of the earthquake's P wave, the spatial distance calculation formula from each station to the epicenter is: ； ； In the formula, represents the number of stations, for The epicentral distance of the station, for The location of the station, is the initial position of the epicenter in the Earth-centered Earth-fixed coordinate system, for The P-wave arrival time of the trigger at the station is equivalent to the earthquake arrival time function, For the moment of earthquake, is the velocity of seismic wave propagation; Eliminate the earthquake occurrence time by using the difference method between stations : ； In the formula, is the total number of triggered stations, and the error equation is obtained by linearization: : ； ； ； ； In the formula, is the design matrix, is the observation vector, To correct the number, , and is the correction for the epicenter position, is the correction for the seismic P-wave velocity.

5. The method for rapid estimation of epicenter position based on strong vibration acceleration records according to claim 4 is characterized in that: After judging that the number of stations triggered by the P wave is greater than or equal to 4, the centroid position of the station is calculated as , using the least squares principle to solve : ； ； In the formula, is the weight matrix extracted from seismic waves at each station, is the identity matrix, is the order of the identity matrix, we get The corrected epicenter position and P-wave velocity are then obtained, and the corrected epicenter position is converted from the Earth-centered Earth-fixed coordinate system to the geodetic coordinate system to obtain the geodetic longitude and latitude and geodetic height of the epicenter position.

6. The method for rapid estimation of epicenter position based on strong vibration acceleration records according to claim 1, characterized in that: The observation equation of the station P wave arrival time is: ； In the formula, for dimensional coefficient matrix, for The true value vector of the seismic element to be estimated. yes The observed value of the P wave arrival time, yes True error vector, construct adjustment factor matrix : ； In the formula, for dimensional observation weight matrix, the relationship between the true error and the observation value is: ； right Make an estimate to obtain the estimated value of the gross error and its position in the observation sequence. Make an adjustment based on the measurement of the equipment to obtain the absolute value of the residual. A pseudo observation, , No. The quasi-observation corresponds to the observation vector , the corresponding true error is , the corresponding weight is , the corresponding coefficient matrix is , Remaining There are gross errors in the observations, The true error corresponding to the observation value with gross error is , introduce additional norm minimum constraint: ； Find a definite solution for the true error estimate: ； ； ； ； In the formula, is the residual vector of the fitting, and yes The two components of For observations with gross errors, The observation vector corresponding to the quasi-observation, , are two parameter matrices.

7. The method for rapid estimation of epicenter position based on strong vibration acceleration records according to claim 6, characterized in that: Solving the estimation of seismic elements after gross error correction includes obtaining the true error estimation: ； When the observed value contains gross errors, the true error estimation has the characteristics of grouping, and the values ​​greater than the threshold are The observations are judged to contain gross errors. Assume that the observation sequence is A rough error, indivual -dimensional unit vector: ； ； In the formula, Corresponding to An observation with gross error, Middle The components are 1 and the rest are 0, and the gross error is expressed as , rewrite the observation equation as: ； ； In the formula, for dimensional coefficient matrix, After separating gross errors , obtained by the least squares criterion Valuation : ； Estimation of seismic elements after gross error correction for: 。 8. The method for rapid estimation of epicenter position based on strong vibration acceleration records according to claim 1, characterized in that: In RANSAC, it is assumed that the data set contains normal values ​​and abnormal values. The normal values ​​are recorded as inliers and the abnormal values ​​are recorded as outliers. The objective function of RANSAC is: ； ； In the formula, is the best fitting function, is the indicator function, when the data point for The value is 1 when it is an interior point, otherwise it is 0. Represents the residual function, expressed as the Euclidean distance from a point to a straight line. For the data set; The identification of internal and external points includes, for a given and mathematical models , the minimum sample set MSS for each sampling is recorded as , the number of samples in MSS is ,exist Randomly select 1 , and according to The sample points in the model are used to calculate the model parameters and fit the mathematical model ; For other sample points in the data set, calculate the residual between the sample and the fitted model, and set the residual threshold to If the difference is less than the threshold, it is an internal point. If the difference is greater than the threshold, it is an external point. The number of internal points is recorded. Repeat the internal and external point identification, calculate the number of internal points and model parameters in this cycle, if the number of internal points this time is greater than the previous number of internal points, save the number of internal points and model parameters calculated this time; otherwise, keep the previous calculation parameters until the maximum number of iterations is met : ； In the formula, for After iterations, the probability that all points are internal points in at least one sampling is set to 99%. Get the probability of the correct model for each iteration, represents the probability of failure of a single iteration; Select the set of model parameters with the largest number of inliers as the optimal model for output, and calculate the inlier rate of the model at the same time : ； In the formula, is the number of external points, is the number of internal points; The preliminary and re-selection methods are consistent with the proposed standard inspection method.

9. The method for rapid estimation of epicenter position based on strong vibration acceleration records according to claim 1, characterized in that: When checking, confirm , and select As the initial number of quasi-observations, calculate and of The value of the moment , check the quasi-observation, calculate : ； Using RANSAC Perform linear model fitting, and the fitting result is ,Will The internal and external point identifications are used as the group identification of the true error, and the internal point set corresponding to RANSAC is selected The corresponding measurement component is the quasi-observation; Recalculate based on the selected quasi-observation and , RANSAC is used again to determine the true error clustering characteristics. If the clustering identification of the internal and external points changes significantly, the quasi-observation is adjusted, otherwise proceed to the next step; For non-quasi-accurate observations, the IGG3 equivalent weight function is constructed, and combined with the selected quasi-accurate observations, the estimated parameters are calculated using the quasi-accurate test. If the difference between the estimated parameters of the two iterations is less than the threshold or reaches the maximum number of iterations, the filtering ends at the current moment; otherwise, the number of iterations is increased by 1, and ; like Then the current moment estimation ends; Otherwise, re-select.

10. The method for rapid estimation of epicenter position based on strong vibration acceleration records according to claim 4, characterized in that: The final calculated earthquake occurrence time is: 。

Images