A microseismic source positioning method and system based on station contribution degree

Abstract

The application discloses a microseismic source positioning method and system based on station contribution degree, relates to the technical field of microseismic intelligent monitoring and source positioning, and first collects P-wave arrival time data recorded by stations and station spatial coordinates and obtains a target microseismic event, then obtains initial source spatial coordinates and an initial time of earthquake occurrence based on the target microseismic event by using a given wave velocity model and calculates arrival time residuals of each station; the microseismic source positioning function with lower data input requirement, clear weight source and strong positioning result feedback capability is realized, the station weight can be automatically obtained according to the current arrival time residual and the spatial geometric condition, the station weight can be further modified according to the actual influence of a single station on the final confidence ellipse area, the function of introducing the station azimuth coverage and the propagation distance balance information in addition to the residual weight can explicitly adjust the geometric influence caused by the uneven distribution of stations, and the positioning iteration stability under the bad initial value is improved.

Description

Technical Field

[0001] This invention relates to the field of intelligent microseismic monitoring and source location technology, specifically to a method and system for microseismic source location based on station contribution. Background Technology

[0002] In scenarios such as monitoring rockbursts and seismic events in coal mines, monitoring hydraulic fracturing in oil and gas reservoirs, and evaluating the stability of tunnels and underground engineering projects, multiple microseismic sensors or seismic geophones are typically deployed around the target area. When a fracture occurs within the rock mass, different stations record the arrival time of the P-wave. Based on the station coordinates, P-wave propagation time, and wave velocity model, the spatial location of the seismic source and the time of origin can be inverted.

[0003] Currently, most existing microseismic source location methods employ the Gauss-Newton method, least squares method, or weighted least squares method. Station weights are typically set based on empirical values, pickup quality, or arrival residuals. After location is established, the source parameter covariance matrix is ​​calculated using the design matrix and observation residuals, and a plane confidence ellipse at a specific confidence level is constructed to quantify the uncertainty range of the source plane location. However, during engineering monitoring, some monitoring systems do not directly provide the P-wave pickup standard deviation or pickup uncertainty for each station. If these parameters are used as weights in the calculation... The requirement for specific inputs increases the data preparation requirements, and relying solely on the residuals at the time of arrival makes it difficult to explicitly reflect the differences in azimuth coverage and propagation distance of a station relative to the current seismic source. Furthermore, the ordinary weighted least squares method may produce a large parameter correction step size when the initial seismic source deviation is large, and the confidence ellipse is given only after the location is completed. This results in the lack of a feedback mechanism to use the confidence domain results to evaluate the role of a single station, which seriously affects the effect of microseismic intelligent monitoring and the accuracy of seismic source location. Therefore, it is necessary to design a microseismic source location method and system based on station contribution. Summary of the Invention

[0004] The purpose of this invention is to overcome the shortcomings of existing technologies and to better and more effectively address the common problems in existing microseismic source location methods during engineering monitoring, such as the fact that some monitoring systems do not directly provide the P-wave picking standard deviation or picking uncertainty for each station. If these parameters are used as necessary inputs for weight calculation, it increases the data preparation requirements. Furthermore, evaluating stations solely based on the time-to-time residuals makes it difficult to explicitly reflect the differences in azimuth coverage and propagation distance between the stations and the current seismic source. Additionally, the ordinary weighted least squares method may produce large parameter correction steps when the initial source deviation is large. Moreover, providing the confidence ellipse only after location is complete lacks a feedback mechanism for using the confidence domain results to evaluate the role of individual stations. This paper addresses the problem of microseismic intelligent monitoring effectiveness and source location accuracy. It proposes a microseismic source location method and system based on station contribution, which achieves microseismic source location with low data input requirements, clear weight sources, and strong feedback capabilities. This system not only allows station weights to be automatically obtained based on current arrival residuals and spatial geometric conditions, but also enables further correction of station weights based on the actual impact of individual stations on the final confidence ellipse area. Furthermore, it incorporates station azimuth coverage and propagation distance balance information beyond residual weights to explicitly adjust the geometric impact caused by uneven station spatial distribution, thus improving the stability of location iteration under unfavorable initial values.

[0005] To achieve the above objectives, the technical solution adopted by the present invention is as follows: A method for locating microseismic sources based on station contribution includes the following steps: Step A: Collect P-wave arrival time data and station spatial coordinates recorded by the stations and obtain the target microseismic event. Then, based on the target microseismic event, use a given wave velocity model to obtain the initial source spatial coordinates and initial earthquake time, and calculate the arrival time residuals of each station. Step B involves linearizing the arrival equation at the current epicenter location and constructing a design matrix. Then, based on the design matrix, the residual consistency factor, azimuth coverage factor, and distance equilibrium factor are calculated and fused to form the station basic weights. Step C: Based on the station's basic weights, construct a weighted least squares equation with damping terms using the design matrix and the current station weight matrix, and solve it to obtain the source parameter correction amount. Then, based on the source parameter correction amount, use a weighted residual objective function to determine whether the parameter correction is effective and adaptively adjust the damping factor to iteratively update the source spatial coordinates and the time of origin. Step D: Calculate the parametric covariance matrix based on the corrected source parameters, design matrix, and current station weight matrix. Then, extract the plane coordinate covariance submatrix from the parametric covariance matrix and obtain the major axis, minor axis, area, and azimuth of the confidence ellipse under the preset confidence level after eigenvalue decomposition. Step E: Reset the current station weight to zero and re-execute steps C and D to obtain the confidence ellipse area after removing the current station. Then, construct the area contribution of the current station based on the change in the confidence ellipse area after including all stations and removing the current station. Step F: Determine whether the current station increases the area of ​​the current confidence ellipse based on the area contribution of the current station and obtain the judgment result. Then, based on the judgment result, reduce the weight of the current station according to the area contribution and obtain the station weight matrix for the next round through smooth update. Step G: Repeat steps B to F until the relative rate of change of the confidence ellipse area, the change of the hypocenter plane position, and the maximum change of the station weights simultaneously meet the preset convergence conditions. The final hypocenter spatial location, earthquake time, hypocenter plane confidence ellipse parameters, area contribution of each station, and final station weights will be obtained.

[0006] The aforementioned microseismic source location method based on station contribution includes step A, which involves acquiring P-wave arrival time data and station spatial coordinates recorded by the stations to obtain the target microseismic event. Then, based on the target microseismic event, a given wave velocity model is used to obtain the initial source spatial coordinates and initial occurrence time, and the arrival time residuals for each station are calculated. The specific steps are as follows. Step A1: Collect P-wave arrival time data and station spatial coordinates recorded by the stations and obtain the target microseismic event. Specifically, collect P-wave arrival time data and station spatial coordinates recorded by no less than 6 valid stations. Step A2: Based on the target microseismic event, the initial source spatial coordinates and initial occurrence time are obtained using a given wave velocity model, and the arrival time residuals of each station are calculated. Specifically, the P-wave arrival time observation equation for the i-th station is shown in formula (1). ; (1) in, For the measured arrival time of the P-wave at the i-th station, The moment of occurrence of the microseismic event. Let be the propagation distance from the epicenter to the i-th station. To preset the P-wave propagation speed, For the time error term, ( , , () represents the coordinates of the earthquake source. , , (i) represents the coordinates of the i-th station.

[0007] The aforementioned microseismic source location method based on station contribution, in step B, involves linearizing the arrival equation at the current source location and constructing a design matrix. Then, based on the design matrix, the residual consistency factor, azimuth coverage factor, and distance equilibrium factor are calculated and fused to form the station's basic weights. The specific steps are as follows. Step B1 involves linearizing the arrival time equations at the current epicenter location and constructing a design matrix. Specifically, this involves performing a first-order linearization of the arrival time equations at the current epicenter location, as shown in formula (2). ; ; (2) in, Let i be the arrival time residual of the i-th station. The arrival time of the theoretical P-wave at the i-th station, calculated based on the current source parameters and a given wave velocity model. , and Let be the partial derivative of the P-wave travel time with respect to the spatial coordinates of the source. , , )and These are the coordinates of the earthquake source to be determined and the correction amount for the time of earthquake occurrence, respectively. Design an N×4 matrix, where N is the total number of stations; Step B2 involves calculating the residual consistency factor, azimuth coverage factor, and distance equilibrium factor based on the design matrix, and then fusing them to form the station's basic weights. Specifically, the residual consistency factor is calculated using the deviation of each station's arrival residual from the median level of the residuals of all participating stations. The azimuth coverage factor is calculated using the azimuth coverage of each station relative to the current seismic source. The distance equilibrium factor is calculated using the deviation of the propagation distance from the median propagation distance of the station group. The specific steps are as follows. Step B21: The residual consistency factor is calculated by the deviation of the arrival time residuals of each station from the median level of the residuals of all participating stations. Specifically, the median of the absolute values ​​of the arrival time residuals of all participating positioning stations is calculated, and the residual consistency factor is obtained based on the magnitude of the residual of the i-th station relative to this median, as shown in formula (3). ; (3) in, This represents the median of the absolute values ​​of the arrival residuals from all participating positioning stations. To perform the median calculation, For the arrival time residual of the Nth station, Let be the residual consistency factor of the i-th station, and let be the residual consistency factor of the i-th station. Follow | |Increases and decreases, To prevent small positive time quantities where the denominator of the residual scale is zero; Step B22: The azimuth coverage factor is calculated based on the azimuth coverage of the stations relative to the current seismic source. Specifically, the stations are sorted according to their plane azimuth relative to the current seismic source, and the average angular interval between the adjacent stations on both sides of the i-th station is calculated to obtain the azimuth coverage factor, as shown in formula (4). ; ; (4) in, Let be the plane azimuth angle of the i-th station relative to the current seismic source. The source plane obtained in the kth round of outer iteration Direction coordinates The source plane obtained in the kth round of outer iteration Direction coordinates Let be the local average azimuth interval corresponding to the i-th station. For the first The plane azimuth of each station relative to the current epicenter. For the first The plane azimuth of each station relative to the current epicenter. For azimuth coverage factor, This is a truncation function. and These are the lower and upper limits of the azimuth coverage factor, respectively; Step B23: The distance balance factor is calculated based on the deviation between the propagation distance and the median propagation distance of the station group. Specifically, the distance balance factor is obtained based on the deviation between the propagation distance of the i-th station and the median propagation distance of all participating stations, as shown in formula (5). (5) in, As a distance balance factor, The median distance of the propagation distance of the positioning stations. For distance scale parameters; Step B24 involves fusing the residual consistency factor, azimuth coverage factor, and distance balance factor to form the station's basic weights, as shown in formula (6). (6) in, This represents the base weight of the i-th station in the k-th outer iteration. Let be the residual consistency factor of the i-th station in the k-th outer iteration. The azimuth coverage factor of the i-th station is expressed by an index. Adjusted weighting terms The distance balance factor of the i-th station is expressed by an exponent. Adjusted weighting terms Let be the residual consistency factor of the j-th station in the k-th outer iteration. The azimuth coverage factor of the j-th station is expressed by the index. Adjusted weighting terms The distance balance factor of the j-th station is expressed by the index. The adjusted weight term, where j is the station number participating in the normalized summation, and These are the adjustment indices for the azimuth coverage factor and the distance balance factor, respectively.

[0008] The aforementioned microseismic source location method based on station contribution includes step C, which involves constructing a weighted least squares equation with a damping term using the design matrix and the current station weight matrix based on the station's basic weights, solving it to obtain the source parameter correction amount, and then using a weighted residual objective function based on the source parameter correction amount to determine whether the parameter correction is effective and adaptively adjusting the damping factor to iteratively update the source spatial coordinates and origin time. The specific steps are as follows. Step C1: Based on the station's basic weights, construct a weighted least squares equation with damping terms using the design matrix and the current station weight matrix, and solve it to obtain the source parameter correction, as shown in formula (7). ; ; ; ; (7) in, For the inner layer The weighted matrix of the next positioning iteration For the first The station weight matrix used in the outer layer iteration, and in the current round of inner layer source location iteration. remain unchanged. The diagonal damping scale matrix is... This is an operation to extract the diagonal elements of a matrix and construct a diagonal matrix. For the inner layer The source parameter correction vector obtained from the second location iteration. The damping factor, For the inner layer The time-delayed residual vector corresponding to the next positioning iteration For the inner layer Source parameter vector after +1 location iterations For the inner layer Source parameter vector before the next localization iteration. It is a parameter vector composed of the spatial coordinates of the earthquake source and the time of occurrence. For transpose; Step C2: Based on the source parameter correction, a weighted residual objective function is used to determine whether the parameter correction is effective and adaptively adjust the damping factor to iteratively update the source spatial coordinates and the time of origin. Specifically, when the updated weighted residual objective function decreases, the source parameter update is accepted and the damping factor is reduced; otherwise, the damping factor is increased and the solution is recalculated, as shown in formula (8). ; (8) in, The objective function is the weighted residual. For the inner layer The damping factor used in the +1 iteration, c↓ is the proportionality coefficient used to decrease the damping factor when the objective function decreases, and c↑ is the proportionality coefficient used to increase the damping factor when the objective function does not decrease. , The weighted residual objective function value is calculated using the updated source parameters.

[0009] The aforementioned microseismic source location method based on station contribution includes step D, which involves calculating the parametric covariance matrix based on the corrected source parameters, design matrix, and current station weight matrix. Then, the planar coordinate covariance submatrix is ​​extracted from the parametric covariance matrix and, after eigenvalue decomposition, the major axis, minor axis, area, and azimuth of the confidence ellipse at a pre-set confidence level are obtained. The specific steps are as follows. Step D1: Calculate the parameter covariance matrix based on the corrected source parameters, design matrix, and current station weight matrix, as shown in formula (9). ; (9) in, For weighted residual variance estimation, The number of unknown parameters and =4, Let covariance be the parameter. This is the station weight matrix used in the current covariance calculation; Step D2 involves extracting the planar coordinate covariance submatrix from the parameter covariance matrix and then performing eigenvalue decomposition to obtain the major axis, minor axis, area, and azimuth of the confidence ellipse at the preset confidence level, as shown in formula (10). ; ; ; ; ; (10) in, The covariance submatrix of the plane coordinates. The epicenter Orientation standard deviation The epicenter Orientation standard deviation The epicenter direction and The covariance term of orientation positioning error, The first submatrix of the covariance of the plane coordinates 1 eigenvector For the eigenvector The corresponding eigenvalues, and All are planar coordinate covariance submatrices eigenvalues, and , and These are the major axis radius and minor axis radius of the confidence ellipse, respectively. For a degree of freedom of 2 and a confidence level of The critical value of the chi-square distribution, The area of ​​the confidence ellipse. The azimuth angle of the confidence ellipse.

[0010] The aforementioned microseismic source location method based on station contribution involves step E, where the current station weight is reset to zero, and steps C and D are re-executed to obtain the area of ​​the confidence ellipse excluding the current station. Then, the area contribution of the current station is constructed based on the change in the area of ​​the confidence ellipse including all stations and excluding the current station. The specific steps are as follows. Step E1: Reset the current station weight to zero and re-execute steps C and D to obtain the area of ​​the confidence ellipse after removing the current station. Specifically, the area contribution of the i-th station is shown in formula (11). (11) in, Let i be the area contribution of the i-th station. For the first The confidence ellipse area obtained from all participating stations is used in the outer layer iteration. The area of ​​the confidence ellipse obtained by repositioning only the i-th station after resetting its rights to zero; Step E2: Construct the area contribution of the current station based on the change in the area of ​​the confidence ellipse after including all stations and excluding the current station. Specifically, if the area contribution of the i-th station... If the value is greater than 0, it indicates that the area of ​​the confidence ellipse decreases after removing the i-th station, meaning that the i-th station has an adverse amplifying effect on the current positioning uncertainty range. If the area contribution of the i-th station is... If the value is ≤0, then no additional contribution feedback weight reduction will be implemented for the i-th station.

[0011] The aforementioned microseismic source location method based on station contribution includes step F, which involves determining whether the current station increases the area of ​​the current confidence ellipse based on its area contribution and obtaining the determination result. Then, based on the determination result, the weight of the current station is reduced according to its area contribution, and the station weight matrix for the next round is obtained through smooth update. The specific steps are as follows. Step F1: Determine whether the current station increases the area of ​​the current confidence ellipse based on the area contribution of the current station and obtain the determination result. The determination result includes area increase and area decrease. If the area contribution of the i-th station indicates that the station increases the area of ​​the current confidence ellipse, then the area increases; otherwise, the area decreases. Step F2 involves reducing the current station weight based on the area contribution according to the judgment result and obtaining the station weight matrix for the next round through smooth update. Specifically, when the judgment result is an increase in area, the current station weight is reduced based on the area contribution and the station weight matrix for the next round is obtained through smooth update, as shown in formula (12). ; (12) in, The weight of the i-th station candidate is obtained after weight reduction based on area contribution feedback. As the lower limit of the weight, This represents the base weight of the i-th station in the k-th outer iteration. The feedback weighting coefficient for area contribution. The weight of the i-th station used in the (k+1)-th outer iteration after smoothing update. This is the weighted smoothing coefficient. Let be the weight of the i-th station in the k-th outer iteration.

[0012] The aforementioned microseismic source location method based on station contribution, in step G, steps B to F are repeated until the relative change rate of the confidence ellipse area, the change in the hypocenter plane position, and the maximum change in the station weight simultaneously satisfy the preset convergence conditions. The final hypocenter spatial location, time of origin, plane confidence ellipse parameters, area contribution of each station, and final station weight are obtained as shown in formula (13). ; ; (13) in, and These are the areas of the confidence ellipses obtained in the (k+1)th and kth rounds of outer iteration, respectively. The threshold for the relative change in the area of ​​the confidence ellipse. and The source planes obtained from the (k+1)th outer iteration are respectively direction and Direction coordinates The threshold for the change in the plane position of the seismic source. This is the threshold for changes in station weights.

[0013] A microseismic source location system based on station contribution includes a microseismic event acquisition module, a station basic weight generation module, a source parameter update module, a plane confidence ellipse calculation module, an area contribution calculation module, a smoothing update module, and a convergence judgment module. The microseismic event acquisition module is used to acquire P-wave arrival time data and station spatial coordinates recorded by the stations and obtain the target microseismic event. Then, based on the target microseismic event, a given wave velocity model is used to obtain the initial source spatial coordinates and the initial earthquake time, and the arrival time residuals of each station are calculated. The station basic weight generation module is used to linearize the arrival equation at the current source location and construct a design matrix. Then, based on the design matrix, it calculates the residual consistency factor, azimuth coverage factor and distance balance factor and integrates them to form the station basic weight. The source parameter update module is used to construct a weighted least squares equation with damping terms based on the station's basic weights, using the design matrix and the current station weight matrix, and solve it to obtain the source parameter correction amount. Then, based on the source parameter correction amount, a weighted residual objective function is used to determine whether the parameter correction is effective and adaptively adjust the damping factor to iteratively update the source spatial coordinates and the time of earthquake occurrence. The planar confidence ellipse calculation module is used to calculate the parameter covariance matrix based on the corrected source parameters, design matrix and current station weight matrix, and then extract the planar coordinate covariance submatrix from the parameter covariance matrix and obtain the major axis, minor axis, area and azimuth of the confidence ellipse under the preset confidence level after eigenvalue decomposition. The area contribution calculation module is used to reset the current station weight to zero and re-execute the source parameter update module and the plane confidence ellipse calculation module to obtain the area of ​​the confidence ellipse after removing the current station. Then, the area contribution of the current station is constructed based on the change in the area of ​​the confidence ellipse after including all stations and removing the current station. The smooth update module is used to determine whether the current station increases the area of ​​the current confidence ellipse based on the area contribution of the current station and obtain the judgment result. Then, based on the judgment result, the weight of the current station is reduced according to the area contribution and the station weight matrix for the next round is obtained through smooth update. The convergence judgment module is used to re-execute the station basic weight generation module, the source parameter update module, the plane confidence ellipse calculation module, the area contribution calculation module, and the smoothing update module until the relative change rate of the confidence ellipse area, the change in the source plane position, and the maximum change in the station weight simultaneously meet the preset convergence conditions. The module will then obtain the final source spatial position, earthquake time, plane confidence ellipse parameters, area contribution of each station, and final station weight.

[0014] The beneficial effects of this invention are: (1) This invention obtains the target microseismic event by collecting P-wave arrival time data and station spatial coordinates recorded by the stations, and then obtains the initial source spatial coordinates and initial earthquake time based on the target microseismic event using a given wave velocity model, and calculates the arrival time residuals of each station. Then, at the current source location, the arrival time equation is linearized and a design matrix is ​​constructed. Based on the design matrix, the residual consistency factor, azimuth coverage factor and distance balance factor are calculated and fused to form the station basic weights. This effectively realizes the microseismic source location method and system with low data input requirements, clear weight sources and strong location result feedback capabilities. It not only enables the station weights to be automatically obtained based on the current arrival time residuals and spatial geometric conditions, but also enables the station weights to be further corrected based on the actual impact of a single station on the final confidence ellipse area.

[0015] (2) This invention constructs a weighted least squares equation with damping terms based on the station's basic weights using the design matrix and the current station weight matrix, and solves it to obtain the source parameter correction amount. Then, based on the source parameter correction amount, a weighted residual objective function is used to determine whether the parameter correction is effective and adaptively adjusts the damping factor to iteratively update the source spatial coordinates and the time of origin. Then, based on the corrected source parameters, the design matrix, and the current station weight matrix, the parameter covariance matrix is ​​calculated. Then, the plane coordinate covariance submatrix is ​​extracted from the parameter covariance matrix and obtained by eigenvalue decomposition to obtain the major axis, minor axis, area, and azimuth of the confidence ellipse under the preset confidence level. This effectively realizes that the microseismic source location method and system do not require the standard deviation or uncertainty of the P-wave picking at each station as a necessary input. It can not only directly use the station coordinates, P-wave arrival time, and location residuals to construct residual consistency weights, but also reduce the dependence on additional data fields of the monitoring system.

[0016] (3) This invention obtains the area of ​​the confidence ellipse after removing the current station by resetting the current station weight to zero and re-executing the source parameter update module and the plane confidence ellipse calculation module. Then, it constructs the area contribution of the current station based on the change in the area of ​​the confidence ellipse after including all stations and removing the current station. Subsequently, it judges whether the current station increases the area of ​​the current confidence ellipse based on the area contribution of the current station and obtains the judgment result. Based on the judgment result, it reduces the weight of the current station according to the area contribution and obtains the station weight matrix for the next round through smooth update. This effectively realizes that the microseismic source location method and system can explicitly adjust the geometric influence caused by uneven spatial distribution of stations by introducing station azimuth coverage and propagation distance balance information in addition to residual weight. Moreover, it can improve the stability of the location iteration under poor initial value by using the damping term to control the source parameter correction step size.

[0017] (4) This invention re-executes the station basic weight generation module, the source parameter update module, the plane confidence ellipse calculation module, and the area contribution calculation module until the relative change rate of the confidence ellipse area, the change in the source plane position, and the maximum change in the station weight simultaneously meet the preset convergence conditions. The final source spatial location, earthquake time, plane confidence ellipse parameters, area contribution of each station, and final station weight are obtained. This effectively realizes that the microseismic source location method and system has a unified index for evaluating the contribution of a single station by using the confidence ellipse area as a single station contribution evaluation. It forms a closed loop through single station virtual elimination, area change evaluation, and feedback weight reduction. At the same time, the exponential bounded weight reduction and weight smoothing of the area contribution can avoid the complete deletion of stations at one time. The triple convergence criterion of area-location-weight can reduce the possibility of premature termination by a single area criterion. It can output not only the source location and confidence ellipse, but also the area contribution of each station and the final weight. This provides a reference for the quality diagnosis and layout optimization of monitoring system stations. Attached Figure Description

[0018] Figure 1 This is an overall flowchart of a microseismic source location method based on station contribution according to the present invention; Figure 2 This is a schematic diagram illustrating the change of the residual consistency factor of the present invention with the normalized residual; Figure 3 This is a schematic diagram illustrating the balance between station azimuth coverage and propagation distance in this invention; Figure 4 This is a schematic diagram of the planar confidence ellipse parameters of the present invention; Figure 5 This is a schematic diagram of the single-station virtual elimination and area contribution calculation of the present invention; Figure 6 This is a schematic diagram of the area contribution feedback weighting and iterative convergence process of the present invention. Detailed Implementation

[0019] The present invention will now be further described with reference to the accompanying drawings.

[0020] like Figure 1 As shown, the present invention provides a microseismic source location method based on station contribution, comprising the following steps: Step A involves acquiring P-wave arrival time data and station spatial coordinates recorded by various stations to obtain the target microseismic event. Then, based on the target microseismic event, a given wave velocity model is used to obtain the initial source spatial coordinates and initial occurrence time, and the arrival time residuals for each station are calculated. The specific steps are as follows. Step A1: Collect P-wave arrival time data and station spatial coordinates recorded by the stations and obtain the target microseismic event. Specifically, collect P-wave arrival time data and station spatial coordinates recorded by no less than 6 valid stations. Step A2: Based on the target microseismic event, the initial source spatial coordinates and initial occurrence time are obtained using a given wave velocity model, and the arrival time residuals of each station are calculated. Specifically, the P-wave arrival time observation equation for the i-th station is shown in formula (1). ; (1) in, For the measured arrival time of the P-wave at the i-th station, The moment of occurrence of the microseismic event. Let be the propagation distance from the epicenter to the i-th station. To preset the P-wave propagation speed, For the time error term, ( , , () represents the coordinates of the earthquake source. , , Let be the coordinates of the i-th station; The initial source location can be obtained using the geometric center of the station, a coarse-grid search, the initial location results of existing microseismic systems, or other available initial values. This embodiment preferably uses a coarse-grid search to obtain the initial source coordinates. , , and the time of the earthquake The theoretical arrival time of each station is calculated at the initial source location, and the initial residual is obtained from the difference between the measured arrival time and the theoretical arrival time.

[0021] Step B involves linearizing the arrival equations at the current epicenter location and constructing a design matrix. Then, based on the design matrix, the residual consistency factor, azimuth coverage factor, and distance equilibrium factor are calculated and fused to form the station's basic weights. The specific steps are as follows: Step B1 involves linearizing the arrival time equations at the current epicenter location and constructing a design matrix. Specifically, this involves performing a first-order linearization of the arrival time equations at the current epicenter location, as shown in formula (2). ; ; (2) in, Let i be the arrival time residual of the i-th station. The arrival time of the theoretical P-wave at the i-th station, calculated based on the current source parameters and a given wave velocity model. , and Let be the partial derivative of the P-wave travel time with respect to the spatial coordinates of the source. , , )and These are the coordinates of the earthquake source to be determined and the correction amount for the time of earthquake occurrence, respectively. Design an N×4 matrix, where N is the total number of stations; Step B2 involves calculating the residual consistency factor, azimuth coverage factor, and distance equilibrium factor based on the design matrix, and then fusing them to form the station's basic weights. Specifically, the residual consistency factor is calculated using the deviation of each station's arrival residual from the median level of the residuals of all participating stations. The azimuth coverage factor is calculated using the azimuth coverage of each station relative to the current seismic source. The distance equilibrium factor is calculated using the deviation of the propagation distance from the median propagation distance of the station group. The specific steps are as follows. like Figure 2 As shown, where Figure 2 The horizontal axis represents the normalized time residual. ,and Figure 2The black dots in the figure represent the sampling points of the residual consistency factor corresponding to different normalized time residuals; in step B21, the residual consistency factor is calculated by the degree of deviation of the time residual of each station from the median level of the residuals of all participating stations. Specifically, the median of the absolute values ​​of the time residuals of all participating positioning stations is calculated and the residual consistency factor is obtained according to the magnitude of the residual of the i-th station relative to the median, as shown in formula (3). ; (3) in, This represents the median of the absolute values ​​of the arrival residuals from all participating positioning stations. To perform the median calculation, For the arrival time residual of the Nth station, Let be the residual consistency factor of the i-th station, and let be the residual consistency factor of the i-th station. Follow | |Increases and decreases, To prevent small positive time quantities where the denominator of the residual scale is zero; This invention does not require inputting the standard deviation of P-wave pickup at each station. To evaluate whether the arrival time of a particular station is consistent with the current positioning result, the median of the absolute values ​​of the residuals of all participating positioning stations is first calculated. Then calculate the residual consistency factor. When the arrival residual of the i-th station | When the level is similar to the overall residual level, Within the normal range; when the residual is reached | When the residual consistency factor is significantly greater than that of most stations, the residual consistency factor is... The median decreases automatically. Since the median is insensitive to a small number of outlier residuals, it is not necessary to know the station picking errors in advance; a relatively consistent assessment can be formed using observations from multiple stations within the event itself.

[0022] like Figure 3 As shown, in step B22, the azimuth coverage factor is calculated based on the azimuth coverage of the stations relative to the current seismic source. Specifically, the stations are sorted according to their plane azimuth angles relative to the current seismic source, and the average angular interval between the adjacent stations on both sides of the i-th station is calculated to obtain the azimuth coverage factor, as shown in formula (4). ; ; (4) in, Let be the plane azimuth angle of the i-th station relative to the current seismic source. The source plane obtained in the kth round of outer iteration Direction coordinates The source plane obtained in the kth round of outer iteration Direction coordinates Let be the local average azimuth interval corresponding to the i-th station. For the first The plane azimuth of each station relative to the current epicenter. For the first The plane azimuth of each station relative to the current epicenter. For azimuth coverage factor, This is a truncation function. and These are the lower and upper limits of the azimuth coverage factor, respectively; Residuals alone cannot reflect the spatial geometry of a station. From Figure 3 As can be seen from this, the present invention calculates the plane azimuth of each station with the current seismic source as the center. Then sort them, and calculate the local azimuth interval based on the angular interval of adjacent stations. This invention normalizes the angular interval into an azimuth coverage factor. When stations are sparse in a certain direction, the independent geometric constraint value of stations in that direction is high; when multiple stations are clustered in similar azimuths, the azimuth coverage factor is high. Relative decrease.

[0023] Step B23: The distance balance factor is calculated based on the deviation between the propagation distance and the median propagation distance of the station group. Specifically, the distance balance factor is obtained based on the deviation between the propagation distance of the i-th station and the median propagation distance of all participating stations, as shown in formula (5). (5) in, As a distance balance factor, The median distance of the propagation distance of the positioning stations. For distance scale parameters; This embodiment utilizes the median propagation distance of all participating stations. Construct distance equilibrium factor This invention maintains a high baseline value for stations with propagation distances close to the typical level of the station group, while moderately suppressing stations with distances significantly far from the overall array.

[0024] Step B24 involves fusing the residual consistency factor, azimuth coverage factor, and distance balance factor to form the station's basic weights, as shown in formula (6). (6) in, This represents the base weight of the i-th station in the k-th outer iteration. Let be the residual consistency factor of the i-th station in the k-th outer iteration. The azimuth coverage factor of the i-th station is expressed by an index. Adjusted weighting terms The distance balance factor of the i-th station is expressed by an exponent. Adjusted weighting terms Let be the residual consistency factor of the j-th station in the k-th outer iteration. The azimuth coverage factor of the j-th station is expressed by the index. Adjusted weighting terms The distance balance factor of the j-th station is expressed by the index. The adjusted weight term, where j is the station number participating in the normalized summation, and These are the adjustment indices for the azimuth coverage factor and the distance balance factor, respectively. In the first outer iteration, the basic weights of each station obtained in step B are used as the initial values ​​for the current station weights, i.e. Then construct the current station weight matrix. ; Will , and Multiplicative fusion is performed to obtain the basic weights of the outer layer in the kth round. ;in and Used to adjust the influence intensity of azimuth coverage and distance balance, the average value after normalization of the basic weights is 1. Preferred Take a value of 0.3 to 1.0. Take a value of 0.1 to 0.8; Take a value of 0.3 to 0.7. Take a value of 1.5 to 3.0.

[0025] Step C involves constructing a weighted least squares equation with a damping term based on the station's basic weights using the design matrix and the current station weight matrix, and solving it to obtain the source parameter correction. Then, based on the source parameter correction, a weighted residual objective function is used to determine the effectiveness of the parameter correction and adaptively adjust the damping factor to iteratively update the source spatial coordinates and origin time. The specific steps are as follows. Step C1: Based on the station's basic weights, construct a weighted least squares equation with damping terms using the design matrix and the current station weight matrix, and solve it to obtain the source parameter correction, as shown in formula (7). ; ; ; ; (7) in, For the inner layer The weighted matrix of the next positioning iteration For the first The station weight matrix used in the outer layer iteration, and in the current round of inner layer source location iteration. remain unchanged. The diagonal damping scale matrix is... This is an operation to extract the diagonal elements of a matrix and construct a diagonal matrix. For the inner layer The source parameter correction vector obtained from the second location iteration. The damping factor, For the inner layer The time-delayed residual vector corresponding to the next positioning iteration For the inner layer Source parameter vector after +1 location iterations For the inner layer Source parameter vector before the next localization iteration. It is a parameter vector composed of the spatial coordinates of the earthquake source and the time of occurrence. For transpose; Step C2: Based on the source parameter correction, a weighted residual objective function is used to determine whether the parameter correction is effective and adaptively adjust the damping factor to iteratively update the source spatial coordinates and the time of origin. Specifically, when the updated weighted residual objective function decreases, the source parameter update is accepted and the damping factor is reduced; otherwise, the damping factor is increased and the solution is recalculated, as shown in formula (8). ; (8) in, The objective function is the weighted residual. For the inner layer The damping factor used in the +1 iteration, c↓ is the proportionality coefficient used to decrease the damping factor when the objective function decreases, and c↑ is the proportionality coefficient used to increase the damping factor when the objective function does not decrease. , The weighted residual objective function value is calculated using the updated source parameters; Repeat steps C1 and C2 to perform source parameter correction, objective function judgment, and damping factor adjustment until the source spatial location correction and the earthquake time correction meet the preset inner-layer convergence conditions or reach the preset maximum number of inner-layer iterations. Then, use the current source parameters as the source location result of the k-th round of outer-layer iteration.

[0026] like Figure 4As shown, step D involves calculating the parametric covariance matrix based on the corrected source parameters, design matrix, and current station weight matrix. Then, the plane coordinate covariance submatrix is ​​extracted from the parametric covariance matrix and, after eigenvalue decomposition, the major axis, minor axis, area, and azimuth of the confidence ellipse at the preset confidence level are obtained. The specific steps are as follows. Step D1: Calculate the parameter covariance matrix based on the corrected source parameters, design matrix, and current station weight matrix, as shown in formula (9). ; (9) in, For weighted residual variance estimation, The number of unknown parameters and =4, Let covariance be the parameter. This is the station weight matrix used in the current covariance calculation; Step D2 involves extracting the planar coordinate covariance submatrix from the parameter covariance matrix and then performing eigenvalue decomposition to obtain the major axis, minor axis, area, and azimuth of the confidence ellipse at the preset confidence level, as shown in formula (10). ; ; ; ; ; (10) in, The covariance submatrix of the plane coordinates. The epicenter Orientation standard deviation The epicenter Orientation standard deviation The epicenter direction and The covariance term of orientation positioning error, The first submatrix of the covariance of the plane coordinates 1 eigenvector For the eigenvector The corresponding eigenvalues, and All are planar coordinate covariance submatrices eigenvalues, and , and These are the major axis radius and minor axis radius of the confidence ellipse, respectively. For a degree of freedom of 2 and a confidence level of The critical value of the chi-square distribution, The area of ​​the confidence ellipse. The azimuth angle of the confidence ellipse.

[0027] like Figure 5 As shown, in step E, the current station weight is reset to zero, and steps C and D are re-executed to obtain the confidence ellipse area after removing the current station. Then, the area contribution of the current station is constructed based on the change in the confidence ellipse area after including all stations and removing the current station. The specific steps are as follows. Step E1: Reset the current station weight to zero and re-execute steps C and D to obtain the area of ​​the confidence ellipse after removing the current station. Specifically, the area contribution of the i-th station is shown in formula (11). (11) in, Let i be the area contribution of the i-th station. For the first The confidence ellipse area obtained from all participating stations is used in the outer layer iteration. The area of ​​the confidence ellipse obtained by repositioning only the i-th station after resetting its rights to zero; To evaluate the actual effect of a single station on the current uncertainty range of source location, this invention employs a single-station virtual elimination strategy. For the i-th station, the current station weight matrix is ​​copied, and only the weight of the i-th station is temporarily set to zero, while the weights of other stations remain unchanged. Then, the damped weighted source location and confidence ellipse calculations are re-executed to obtain... This operation is for mathematical evaluation only and does not represent on-site physical removal or permanent decommissioning of the station.

[0028] Step E2: Construct the area contribution of the current station based on the change in the area of ​​the confidence ellipse after including all stations and excluding the current station. Specifically, if the area contribution of the i-th station... If the value is greater than 0, it indicates that the area of ​​the confidence ellipse decreases after removing the i-th station, meaning that the i-th station has an adverse amplifying effect on the current positioning uncertainty range. If the area contribution of the i-th station is... If the value is ≤0, then no additional contribution feedback weight reduction will be implemented for the i-th station.

[0029] like Figure 6 As shown, in step F, it is determined whether the current station increases the area of ​​the current confidence ellipse based on its area contribution, and the determination result is obtained. Then, based on the determination result, the weight of the current station is reduced according to its area contribution, and the station weight matrix for the next round is obtained through smoothing update. The specific steps are as follows. Step F1: Determine whether the current station increases the area of ​​the current confidence ellipse based on the area contribution of the current station and obtain the determination result. The determination result includes area increase and area decrease. If the area contribution of the i-th station indicates that the station increases the area of ​​the current confidence ellipse, then the area increases; otherwise, the area decreases. Step F2 involves reducing the current station weight based on the area contribution according to the judgment result and obtaining the station weight matrix for the next round through smooth update. Specifically, when the judgment result is an increase in area, the current station weight is reduced based on the area contribution and the station weight matrix for the next round is obtained through smooth update, as shown in formula (12). ; (12) in, The weight of the i-th station candidate is obtained after weight reduction based on area contribution feedback. As the lower limit of the weight, This represents the base weight of the i-th station in the k-th outer iteration. The feedback weighting coefficient for area contribution. The weight of the i-th station used in the (k+1)-th outer iteration after smoothing update. This is the weighted smoothing coefficient. Let be the weight of the i-th station in the k-th outer iteration.

[0030] For the The outer layer iteration is based on the base weights of that round. Construct candidate weights with area contribution, and then combine the candidate weights with the first... The current weight actually used in the round Perform smooth fusion, resulting in The current station weight will be used for the next round of outer layer iteration.

[0031] To prevent abrupt changes in weights due to a single area contribution evaluation, the target weights and current weights are smoothly merged. Weight smoothing coefficient. The smaller the value, the slower the weight changes; weight smoothing coefficient The larger the value, the faster the response to the current event's contribution evaluation, and the better. Take a value of 0.2 to 0.7.

[0032] Step G: Repeat steps B to F until the relative rate of change of the confidence ellipse area, the change of the hypocenter plane position, and the maximum change of the station weights simultaneously meet the preset convergence conditions. The final hypocenter spatial location, time of origin, plane confidence ellipse parameters, area contribution of each station, and final station weights are obtained as shown in formula (13). ; ; (13) in, and These are the areas of the confidence ellipses obtained in the (k+1)th and kth rounds of outer iteration, respectively. The threshold for the relative change in the area of ​​the confidence ellipse. and The source planes obtained from the (k+1)th outer iteration are respectively direction and Direction coordinates The threshold for the change in the plane position of the seismic source. The threshold for changes in station weights; Only when the relative rate of change of the confidence ellipse area is less than the relative rate of change of the confidence ellipse area threshold. The change in the plane position of the source of two adjacent earthquakes is less than the threshold for the change in the plane position of the source. And the maximum change in the weight of all stations is less than the threshold for the change in station weight. The outer iteration terminates only when the relative change threshold of the confidence ellipse area is reached. The threshold for changes in the plane position of the seismic source is set at 1% to 5%. The threshold for station weight changes is set to 0.5–5 m. Take a value of 0.01 to 0.05.

[0033] Figure 6 The convergence threshold marked in the figure corresponds to the preset convergence condition in step G, namely the confidence ellipse area relative change threshold. Threshold for change in the plane position of the seismic source Threshold for changes in station weights . Figure 6 The normalized convergence index and threshold line in the formula are used to illustrate the iterative convergence process. The actual calculation still judges the above three changes according to formula (13), and terminates the outer iteration only when all three conditions are met.

[0034] A microseismic source location system based on station contribution includes a microseismic event acquisition module, a station basic weight generation module, a source parameter update module, a plane confidence ellipse calculation module, an area contribution calculation module, a smoothing update module, and a convergence judgment module. The microseismic event acquisition module is used to acquire P-wave arrival time data and station spatial coordinates recorded by the stations and obtain the target microseismic event. Then, based on the target microseismic event, a given wave velocity model is used to obtain the initial source spatial coordinates and the initial earthquake time, and the arrival time residuals of each station are calculated. The station basic weight generation module is used to linearize the arrival equation at the current source location and construct a design matrix. Then, based on the design matrix, it calculates the residual consistency factor, azimuth coverage factor and distance balance factor and integrates them to form the station basic weight. The source parameter update module is used to construct a weighted least squares equation with damping terms based on the station's basic weights, using the design matrix and the current station weight matrix, and solve it to obtain the source parameter correction amount. Then, based on the source parameter correction amount, a weighted residual objective function is used to determine whether the parameter correction is effective and adaptively adjust the damping factor to iteratively update the source spatial coordinates and the time of earthquake occurrence. The planar confidence ellipse calculation module is used to calculate the parameter covariance matrix based on the corrected source parameters, design matrix and current station weight matrix, and then extract the planar coordinate covariance submatrix from the parameter covariance matrix and obtain the major axis, minor axis, area and azimuth of the confidence ellipse under the preset confidence level after eigenvalue decomposition. The area contribution calculation module is used to reset the current station weight to zero and re-execute the source parameter update module and the plane confidence ellipse calculation module to obtain the area of ​​the confidence ellipse after removing the current station. Then, the area contribution of the current station is constructed based on the change in the area of ​​the confidence ellipse after including all stations and removing the current station. The smooth update module is used to determine whether the current station increases the area of ​​the current confidence ellipse based on the area contribution of the current station and obtain the judgment result. Then, based on the judgment result, the weight of the current station is reduced according to the area contribution and the station weight matrix for the next round is obtained through smooth update. The convergence judgment module is used to re-execute the station basic weight generation module, the source parameter update module, the plane confidence ellipse calculation module, the area contribution calculation module, and the smoothing update module until the relative change rate of the confidence ellipse area, the change in the source plane position, and the maximum change in the station weight simultaneously meet the preset convergence conditions. The module will then obtain the final source spatial position, earthquake time, plane confidence ellipse parameters, area contribution of each station, and final station weight.

[0035] In summary, the microseismic source location method and system based on station contribution of the present invention first acquires P-wave arrival time data and station spatial coordinates recorded by stations to obtain the target microseismic event. Then, based on the target microseismic event, a given wave velocity model is used to obtain the initial source spatial coordinates and initial occurrence time, and the arrival time residuals of each station are calculated. Next, the arrival time equation is linearized at the current source location, and a design matrix is ​​constructed. Then, based on the design matrix, the residual consistency factor, azimuth coverage factor, and distance equilibrium factor are calculated and fused to form the station basic weights. Subsequently, based on the station basic weights, a weighted least squares equation with damping terms is constructed using the design matrix and the current station weight matrix, and the source parameter correction is obtained by solving the equation. Finally, a weighted least squares equation is used based on the source parameter correction. The residual objective function determines the effectiveness of parameter correction and adaptively adjusts the damping factor to iteratively update the source spatial coordinates and origin time. Then, based on the corrected source parameters, design matrix, and current station weight matrix, the parameter covariance matrix is ​​calculated. The plane coordinate covariance submatrix is ​​extracted from the parameter covariance matrix and, after eigenvalue decomposition, the major axis, minor axis, area, and azimuth of the confidence ellipse at a preset confidence level are obtained. Next, the current station weight is reset to zero, and the source parameter update module and plane confidence ellipse calculation module are re-executed to obtain the area of ​​the confidence ellipse after removing the current station. Then, the area contribution of the current station is constructed based on the change in the area of ​​the confidence ellipse including all stations and excluding the current station. Finally, the area contribution of the current station is used to determine the current... The system determines whether a station increases the area of ​​the current confidence ellipse and obtains a judgment result. Based on the judgment result, the weight of the current station is reduced according to the area contribution, and the station weight matrix for the next round is obtained through smooth updates. Finally, the station basic weight generation module, the source parameter update module, the plane confidence ellipse calculation module, and the area contribution calculation module are re-executed until the relative change rate of the confidence ellipse area, the change in the source plane position, and the maximum change in the station weight simultaneously meet the preset convergence conditions. The system then obtains the final source spatial location, earthquake time, plane confidence ellipse parameters, area contribution of each station, and final station weight. This effectively realizes that the microseismic source location method and system have low data input requirements, clear weight sources, and strong location result feedback capabilities. The source location function not only enables the station weights to be automatically obtained based on the current arrival time residuals and spatial geometric conditions, but also allows the station weights to be further corrected based on the actual impact of a single station on the final confidence ellipse area. It also achieves the function of not requiring the standard deviation or uncertainty of P-wave picking for each station as a necessary input. It can directly construct residual consistency weights using station coordinates, P-wave arrival time, and location residuals, and reduces the dependence on additional data fields of the monitoring system. It also achieves the function of explicitly adjusting the geometric impact caused by uneven spatial distribution of stations by introducing station azimuth coverage and propagation distance balance information in addition to residual weights. Furthermore, it can improve the stability of location iteration under adverse initial values ​​by controlling the source parameter correction step size using the damping term.This invention uses the area of ​​the confidence ellipse as a unified indicator for evaluating the contribution of individual stations. It forms a closed loop through virtual removal of individual stations, evaluation of area changes, and feedback weight reduction. Furthermore, the exponential bounded weight reduction and weight smoothing of area contribution prevent the complete deletion of stations at once. The triple convergence criterion of area, location, and weight reduces the possibility of premature termination based on a single area criterion. It can output not only the source location and confidence ellipse but also the area contribution and final weight of each station, providing a reference for the quality diagnosis and deployment optimization of monitoring system stations.

[0036] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of the present invention is defined by the appended claims and their equivalents.

Claims

1. A method for locating microseismic sources based on station contribution, characterized in that: Includes the following steps, Step A: Collect P-wave arrival time data and station spatial coordinates recorded by the stations and obtain the target microseismic event. Then, based on the target microseismic event, use a given wave velocity model to obtain the initial source spatial coordinates and initial earthquake time, and calculate the arrival time residuals of each station. Step B involves linearizing the arrival equation at the current epicenter location and constructing a design matrix. Then, based on the design matrix, the residual consistency factor, azimuth coverage factor, and distance equilibrium factor are calculated and fused to form the station basic weights. Step C: Based on the station's basic weights, construct a weighted least squares equation with damping terms using the design matrix and the current station weight matrix, and solve it to obtain the source parameter correction amount. Then, based on the source parameter correction amount, use a weighted residual objective function to determine whether the parameter correction is effective and adaptively adjust the damping factor to iteratively update the source spatial coordinates and the time of origin. Step D: Calculate the parametric covariance matrix based on the corrected source parameters, design matrix, and current station weight matrix. Then, extract the plane coordinate covariance submatrix from the parametric covariance matrix and obtain the major axis, minor axis, area, and azimuth of the confidence ellipse under the preset confidence level after eigenvalue decomposition. Step E: Reset the current station weight to zero and re-execute steps C and D to obtain the confidence ellipse area after removing the current station. Then, construct the area contribution of the current station based on the change in the confidence ellipse area after including all stations and removing the current station. Step F: Determine whether the current station increases the area of ​​the current confidence ellipse based on the area contribution of the current station and obtain the judgment result. Then, based on the judgment result, reduce the weight of the current station according to the area contribution and obtain the station weight matrix for the next round through smooth update. Step G: Repeat steps B to F until the relative rate of change of the confidence ellipse area, the change of the hypocenter plane position, and the maximum change of the station weights simultaneously meet the preset convergence conditions. The final hypocenter spatial location, earthquake time, hypocenter plane confidence ellipse parameters, area contribution of each station, and final station weights will be obtained.

2. The microseismic source location method based on station contribution as described in claim 1, characterized in that: Step A involves acquiring P-wave arrival time data and station spatial coordinates recorded by various stations to obtain the target microseismic event. Then, based on the target microseismic event, a given wave velocity model is used to obtain the initial source spatial coordinates and initial occurrence time, and the arrival time residuals for each station are calculated. The specific steps are as follows. Step A1: Collect P-wave arrival time data and station spatial coordinates recorded by the stations and obtain the target microseismic event. Specifically, collect P-wave arrival time data and station spatial coordinates recorded by no less than 6 valid stations. Step A2: Based on the target microseismic event, the initial source spatial coordinates and initial occurrence time are obtained using a given wave velocity model, and the arrival time residuals of each station are calculated. Specifically, the P-wave arrival time observation equation for the i-th station is shown in formula (1). ； （1） in, For the measured arrival time of the P-wave at the i-th station, The moment of occurrence of the microseismic event. Let be the propagation distance from the epicenter to the i-th station. To preset the P-wave propagation speed, For the time error term, ( , , () represents the coordinates of the earthquake source. , , (i) represents the coordinates of the i-th station.

3. The microseismic source location method based on station contribution as described in claim 2, characterized in that: Step B involves linearizing the arrival equations at the current epicenter location and constructing a design matrix. Then, based on the design matrix, the residual consistency factor, azimuth coverage factor, and distance equilibrium factor are calculated and fused to form the station's basic weights. The specific steps are as follows: Step B1 involves linearizing the arrival time equations at the current epicenter location and constructing a design matrix. Specifically, this involves performing a first-order linearization of the arrival time equations at the current epicenter location, as shown in formula (2). ； ； （2） in, Let i be the arrival time residual of the i-th station. The arrival time of the theoretical P-wave at the i-th station, calculated based on the current source parameters and a given wave velocity model. , and Let be the partial derivative of the P-wave travel time with respect to the spatial coordinates of the source. , , )and These are the coordinates of the earthquake source to be determined and the correction amount for the time of earthquake occurrence, respectively. Design an N×4 matrix, where N is the total number of stations; Step B2 involves calculating the residual consistency factor, azimuth coverage factor, and distance equilibrium factor based on the design matrix, and then fusing them to form the station's basic weights. Specifically, the residual consistency factor is calculated using the deviation of each station's arrival residual from the median level of the residuals of all participating stations. The azimuth coverage factor is calculated using the azimuth coverage of each station relative to the current seismic source. The distance equilibrium factor is calculated using the deviation of the propagation distance from the median propagation distance of the station group. The specific steps are as follows. Step B21: The residual consistency factor is calculated by the degree of deviation of the arrival time residuals of each station from the median level of the residuals of all participating stations. Specifically, the median of the absolute values ​​of the arrival time residuals of all participating positioning stations is calculated, and the residual consistency factor is obtained based on the magnitude of the residual of the i-th station relative to this median, as shown in formula (3). ； （3） in, This represents the median of the absolute values ​​of the arrival residuals from all participating positioning stations. To perform the median calculation, For the arrival time residual of the Nth station, Let be the residual consistency factor of the i-th station, and let be the residual consistency factor of the i-th station. Follow | |Increases and decreases, To prevent small positive time quantities where the denominator of the residual scale is zero; Step B22: The azimuth coverage factor is calculated based on the azimuth coverage of the stations relative to the current seismic source. Specifically, the stations are sorted according to their plane azimuth relative to the current seismic source, and the average angular interval between the adjacent stations on both sides of the i-th station is calculated to obtain the azimuth coverage factor, as shown in formula (4). ； ； （4） in, Let be the plane azimuth angle of the i-th station relative to the current seismic source. The source plane obtained in the kth round of outer iteration Direction coordinates The source plane obtained in the kth round of outer iteration Direction coordinates Let be the local average azimuth interval corresponding to the i-th station. For the first The plane azimuth of each station relative to the current epicenter. For the first The plane azimuth of each station relative to the current epicenter. This is the azimuth coverage factor. This is a truncation function. and These are the lower and upper limits of the azimuth coverage factor, respectively; Step B23: The distance balance factor is calculated based on the deviation between the propagation distance and the median propagation distance of the station group. Specifically, the distance balance factor is obtained based on the deviation between the propagation distance of the i-th station and the median propagation distance of all participating stations, as shown in formula (5). （5） in, As a distance balance factor, The median distance of the propagation distance of the positioning stations. For distance scale parameters; Step B24 involves fusing the residual consistency factor, azimuth coverage factor, and distance balance factor to form the station's basic weights, as shown in formula (6). （6） in, This represents the base weight of the i-th station in the k-th outer iteration. Let be the residual consistency factor of the i-th station in the k-th outer iteration. The azimuth coverage factor of the i-th station is expressed by an index. Adjusted weighting terms The distance balance factor of the i-th station is expressed by an exponent. Adjusted weighting terms Let be the residual consistency factor of the j-th station in the k-th outer iteration. The azimuth coverage factor of the j-th station is expressed by the index. Adjusted weighting terms The distance balance factor of the j-th station is expressed by the index. The adjusted weight term, where j is the station number participating in the normalized summation, and These are the adjustment indices for the azimuth coverage factor and the distance balance factor, respectively.

4. The microseismic source location method based on station contribution as described in claim 3, characterized in that: Step C involves constructing a weighted least squares equation with a damping term based on the station's basic weights using the design matrix and the current station weight matrix, and solving it to obtain the source parameter correction. Then, based on the source parameter correction, a weighted residual objective function is used to determine the effectiveness of the parameter correction and adaptively adjust the damping factor to iteratively update the source spatial coordinates and origin time. The specific steps are as follows. Step C1: Based on the station's basic weights, construct a weighted least squares equation with damping terms using the design matrix and the current station weight matrix, and solve it to obtain the source parameter correction, as shown in formula (7). ； ； ； ； （7） in, For the inner layer The weighted matrix of the next positioning iteration For the first The station weight matrix used in the outer layer iteration, and in the current round of inner layer source location iteration. remain unchanged. The diagonal damping scale matrix is... This is an operation to extract the diagonal elements of a matrix and construct a diagonal matrix. For the inner layer The source parameter correction vector obtained from the second location iteration. The damping factor, For the inner layer The time-delayed residual vector corresponding to the next positioning iteration For the inner layer Source parameter vector after +1 location iterations For the inner layer Source parameter vector before the next localization iteration. It is a parameter vector composed of the spatial coordinates of the earthquake source and the time of occurrence. For transpose; Step C2 involves using a weighted residual objective function based on the source parameter correction to determine the effectiveness of the parameter correction and adaptively adjusting the damping factor to iteratively update the source spatial coordinates and the time of origin. Specifically, when the updated weighted residual objective function decreases, the source parameter update is accepted and the damping factor is reduced; otherwise, the damping factor is increased and the solution is recalculated, as shown in formula (8). ； （8） in, The objective function is the weighted residual. For the inner layer The damping factor used in the +1 iteration, c↓ is the proportionality coefficient used to decrease the damping factor when the objective function decreases, and c↑ is the proportionality coefficient used to increase the damping factor when the objective function does not decrease. , The weighted residual objective function value is calculated using the updated source parameters.

5. The microseismic source location method based on station contribution as described in claim 4, characterized in that: Step D involves calculating the parametric covariance matrix based on the corrected source parameters, design matrix, and current station weight matrix. Then, the plane coordinate covariance submatrix is ​​extracted from the parametric covariance matrix and, after eigenvalue decomposition, the major axis, minor axis, area, and azimuth of the confidence ellipse at the preset confidence level are obtained. The specific steps are as follows: Step D1: Calculate the parameter covariance matrix based on the corrected source parameters, design matrix, and current station weight matrix, as shown in formula (9). ； （9） in, For weighted residual variance estimation, The number of unknown parameters and =4, Let covariance be the parameter. This is the station weight matrix used in the current covariance calculation; Step D2 involves extracting the planar coordinate covariance submatrix from the parameter covariance matrix and then performing eigenvalue decomposition to obtain the major axis, minor axis, area, and azimuth of the confidence ellipse at the preset confidence level, as shown in formula (10). ； ； ； ； ； （10） in, The covariance submatrix of the plane coordinates. The epicenter Orientation standard deviation The epicenter Orientation standard deviation The epicenter direction and The covariance term of the orientation positioning error, The first submatrix of the covariance of the plane coordinates 1 eigenvector For the eigenvector The corresponding eigenvalues, and All are planar coordinate covariance submatrices eigenvalues, and , and These are the major axis radius and minor axis radius of the confidence ellipse, respectively. For a degree of freedom of 2 and a confidence level of The critical value of the chi-square distribution, The area of ​​the confidence ellipse. The azimuth angle of the confidence ellipse.

6. The microseismic source location method based on station contribution as described in claim 5, characterized in that: Step E: Reset the current station weight to zero and re-execute steps C and D to obtain the confidence ellipse area after removing the current station. Then, construct the area contribution of the current station based on the change in the confidence ellipse area after including all stations and removing the current station. The specific steps are as follows: Step E1: Reset the current station weight to zero and re-execute steps C and D to obtain the area of ​​the confidence ellipse after removing the current station. Specifically, the area contribution of the i-th station is shown in formula (11). （11） in, Let i be the area contribution of the i-th station. For the first The confidence ellipse area obtained from all participating stations is used in the outer layer iteration. The area of ​​the confidence ellipse obtained by repositioning only the i-th station after resetting its rights to zero; Step E2: Construct the area contribution of the current station based on the change in the area of ​​the confidence ellipse after including all stations and excluding the current station. Specifically, if the area contribution of the i-th station... If the value is greater than 0, it indicates that the area of ​​the confidence ellipse decreases after removing the i-th station, meaning that the i-th station has an adverse amplifying effect on the current positioning uncertainty range. If the area contribution of the i-th station is... If the value is ≤0, then no additional contribution feedback weight reduction will be implemented for the i-th station.

7. The microseismic source location method based on station contribution as described in claim 6, characterized in that: Step F involves determining whether the current station increases the area of ​​the current confidence ellipse based on its area contribution and obtaining the determination result. Then, based on the determination result, the weight of the current station is reduced according to its area contribution, and the station weight matrix for the next round is obtained through smoothing updates. The specific steps are as follows: Step F1: Determine whether the current station increases the area of ​​the current confidence ellipse based on the area contribution of the current station and obtain the determination result. The determination result includes area increase and area decrease. If the area contribution of the i-th station indicates that the station increases the area of ​​the current confidence ellipse, then the area increases; otherwise, the area decreases. Step F2 involves reducing the current station weight based on the area contribution according to the judgment result and obtaining the station weight matrix for the next round through smooth update. Specifically, when the judgment result is an increase in area, the current station weight is reduced based on the area contribution and the station weight matrix for the next round is obtained through smooth update, as shown in formula (12). ； （12） in, The weight of the i-th station candidate is obtained after weight reduction based on area contribution feedback. As the lower limit of the weight, This represents the base weight of the i-th station in the k-th outer iteration. The feedback weighting coefficient for area contribution. The weight of the i-th station used in the (k+1)-th outer iteration after smoothing update. This is the weighted smoothing coefficient. Let be the weight of the i-th station in the k-th outer iteration.

8. The microseismic source location method based on station contribution as described in claim 7, characterized in that: Step G: Repeat steps B to F until the relative rate of change of the confidence ellipse area, the change of the hypocenter plane position, and the maximum change of the station weights simultaneously meet the preset convergence conditions. The final hypocenter spatial location, time of origin, plane confidence ellipse parameters, area contribution of each station, and final station weights are obtained as shown in formula (13). ； ； （13） in, and These are the areas of the confidence ellipse obtained in the (k+1)th and kth rounds of outer iteration, respectively. The threshold for the relative change in the area of ​​the confidence ellipse. and The source planes obtained from the (k+1)th outer iteration are respectively direction and Direction coordinates The threshold for the change in the plane position of the seismic source. This is the threshold for changes in station weights.

9. A microseismic source location system based on station contribution, wherein the specific location process of the microseismic source location system is based on the microseismic source location method according to any one of claims 1-8, characterized in that: It includes a microseismic event acquisition module, a station basic weight generation module, a source parameter update module, a plane confidence ellipse calculation module, an area contribution calculation module, a smoothing update module, and a convergence judgment module. The microseismic event acquisition module is used to acquire P-wave arrival time data and station spatial coordinates recorded by the stations and obtain the target microseismic event. Then, based on the target microseismic event, a given wave velocity model is used to obtain the initial source spatial coordinates and the initial earthquake time, and calculate the arrival time residuals of each station. The station basic weight generation module is used to linearize the arrival equation at the current source location and construct a design matrix. Then, based on the design matrix, it calculates the residual consistency factor, azimuth coverage factor and distance balance factor and integrates them to form the station basic weight. The source parameter update module is used to construct a weighted least squares equation with damping terms based on the station's basic weights, using the design matrix and the current station weight matrix, and solve it to obtain the source parameter correction amount. Then, based on the source parameter correction amount, a weighted residual objective function is used to determine whether the parameter correction is effective and adaptively adjust the damping factor to iteratively update the source spatial coordinates and the time of earthquake occurrence. The planar confidence ellipse calculation module is used to calculate the parameter covariance matrix based on the corrected source parameters, design matrix and current station weight matrix, and then extract the planar coordinate covariance submatrix from the parameter covariance matrix and obtain the major axis, minor axis, area and azimuth of the confidence ellipse under the preset confidence level after eigenvalue decomposition. The area contribution calculation module is used to reset the current station weight to zero and re-execute the source parameter update module and the plane confidence ellipse calculation module to obtain the area of ​​the confidence ellipse after removing the current station. Then, the area contribution of the current station is constructed based on the change in the area of ​​the confidence ellipse after including all stations and removing the current station. The smooth update module is used to determine whether the current station increases the area of ​​the current confidence ellipse based on the area contribution of the current station and obtain the judgment result. Then, based on the judgment result, the weight of the current station is reduced according to the area contribution and the station weight matrix for the next round is obtained through smooth update. The convergence judgment module is used to re-execute the station basic weight generation module, the source parameter update module, the plane confidence ellipse calculation module, the area contribution calculation module, and the smoothing update module until the relative change rate of the confidence ellipse area, the change in the source plane position, and the maximum change in the station weight simultaneously meet the preset convergence conditions. The module will then obtain the final source spatial position, earthquake time, plane confidence ellipse parameters, area contribution of each station, and final station weight.

Images