Method for calculating position of seismic station based on seismic phase observation report

Abstract

The invention discloses a method for calculating the position of a seismic station based on a seismic phase observation report. The method comprises the steps of seismic phase observation report data preprocessing, seismic station longitude and latitude information calculation, seismic station elevation information extraction and seismic station position information accuracy analysis. According to the method, the accurate latitude and longitude and elevation information of the station at the unknown position is obtained through the seismic event position, the multiple epicentral distance and inverse azimuth data and the surface relief grid data in the seismic phase observation report in combination with the steps of data preprocessing and accuracy analysis, so that the real position of the seismic station can be accurately obtained under the condition that the real position of the seismic station is difficult to directly obtain. And the richness of a data set in seismological research is improved as much as possible.

Description

Technical Field

[0001] The present invention belongs to the technical field of earthquake monitoring, and in particular relates to a method for calculating the position of a seismic station based on a seismic phase observation report. Background Art

[0002] Seismic station locations (including elevation information) are critical parameters in basic seismological research, and their accuracy directly determines the reliability of all subsequent related research. However, because seismic station locations are not fully public data, researchers often only have access to earthquake phase observation reports without the corresponding seismic station location information. This results in the abandonment of these observational data, reducing the data's usability and severely limiting the resolution and accuracy of subsequent related research.

[0003] Seismic phase observation reports contain three types of information: earthquake event location, epicentral distance, and back azimuth. Theoretically, the geographic coordinates of the station can be derived by jointly calculating these three parameters, providing an effective method for obtaining the latitude and longitude information of stations at unknown locations. However, a systematic solution has not yet been formed in this regard. There has been no significant progress in either the analytical calculation based on a single earthquake event or the joint solution algorithm for multiple events, which is not conducive to further improving the depth and breadth of seismological research. Therefore, it is very necessary to construct a complete system of accurate inversion calculation methods for earthquake station location parameters.

[0004] The patent application document with announcement number CN114721049A discloses a method for earthquake positioning with the azimuth of a virtual seismic station participating in an off-grid network. The method includes: determining the position of the virtual seismic station based on the rough epicenter position and the actual distribution of seismic observation stations; determining the back azimuth of each virtual seismic station to each earthquake event based on the position of each virtual seismic station and the rough epicenter coordinates of the earthquake; integrating the earthquake phase arrival time and back azimuth of the physical seismic station and the back azimuth of the virtual seismic station, and using a linear positioning method for positioning to complete the earthquake positioning with the azimuth of the virtual seismic station participating in the off-grid network.

[0005] The patent application document with publication number CN118707592A discloses a single-station earthquake location method based on staggered grid search. By designing multiple sets of grids that intersect each other in space, the spatial sampling rate of earthquake location is improved, making the location results more stable and accurate, and the reliability of the location results is proved by calculating the standard deviation of multiple sets of results.

[0006] The above patent application documents use different methods to locate earthquake positions using station data with known longitude and latitude positions, which falls within the scope of earthquake positioning work. However, the premise of earthquake positioning is to obtain accurate longitude and latitude information of the stations. However, sometimes due to the lack of information sources, the longitude and latitude information of some stations is difficult to obtain, which makes earthquake positioning difficult or the amount of data is too small. Unlike the technical problem to be solved in this application, which is how to determine the accurate longitude and latitude position information of earthquake stations from public observation report data when the longitude and latitude information of earthquake stations is difficult to obtain, the method of this application is a prerequisite for earthquake positioning. Summary of the Invention

[0007] In order to solve the above technical problems, the present invention provides a method for calculating the location of a seismic station based on seismic phase observation reports.

[0008] The present invention is achieved through the following technical solutions.

[0009] The present invention provides a method for calculating the location of a seismic station based on a seismic phase observation report, comprising the following steps:

[0010] A1: Preprocessing of seismic phase observation report data: To ensure the accuracy of seismic phase observation report data used in station location calculations, preprocessing is required to improve the accuracy of the basic data set and, in turn, ensure the credibility of the final calculation results. Set a longitude and latitude range, collect earthquake event locations, multiple epicentral distance and back-azimuth data, and surface relief grid data within this longitude and latitude range, and record known and unknown seismic stations within this longitude and latitude range.

[0011] A2: Calculate the longitude and latitude information of seismic stations: By sorting and extracting the longitude and latitude of earthquake events recorded by each seismic station in the seismic phase observation report and the corresponding epicenter distance and back azimuth parameters, a global inversion algorithm is used to obtain the longitude and latitude of each seismic station calculated for each earthquake event. Finally, the longitude and latitude information of the same seismic station obtained from multiple earthquake events is merged and sorted to obtain the accurate longitude and latitude information of each seismic station.

[0012] A3: Extracting seismic station elevation information: Based on the calculated seismic station longitude and latitude information, extract elevation information from the corresponding location in the surface relief grid data. If the station longitude and latitude are not on the grid point, interpolate the elevation information from the four nearby grid points.

[0013] A4: Accuracy analysis of seismic station location information: Based on the seismic stations whose location information has been disclosed in the seismic phase observation report, the calculation results are verified and analyzed for accuracy, so as to quantitatively determine the accuracy of the final calculation results.

[0014] Preferably, in step A1, the epicentral distance and phase arrival time data are extracted from the phase observation report, and the phase travel time is calculated by the difference between the phase arrival time and the earthquake occurrence time. Then, the theoretical time-distance curve is obtained by linear fitting, and different theoretical phase travel times are obtained according to different epicentral distances. Subsequently, the actual travel time data of the phase observation report is subtracted from the corresponding theoretical phase travel time data to obtain arrival time residual data. Finally, a suitable threshold is selected to remove data with excessive residuals, thereby completing the data preprocessing.

[0015] In step A1, the data can be further screened by limiting the minimum number of phases and station gap angles recorded by each station involved in the calculation, thereby further ensuring the accuracy of the final calculated seismic station position.

[0016] Preferably, the data is fitted with least squares by a linear fitting formula to obtain a theoretical time-distance curve, and the fitting formula is:

[0017]

[0018] t=ad+b (3)

[0019] Among them, d is the independent variable of the epicenter distance of the theoretical time-distance curve, t is the dependent variable of the earthquake phase travel time of the theoretical time-distance curve, a is the slope obtained by fitting, b is the intercept obtained by fitting, and x i is the epicenter distance of the earthquake event corresponding to the i-th phase in the data involved in the fitting, y i is the travel time of the ith phase, and n is the number of phases.

[0020] Preferably, in step A1, all the obtained travel time residual data are subjected to normal distribution statistics, with P wave travel time data being prioritized to obtain their mean and standard deviation, and 2 times the standard deviation is selected as the data screening threshold, thereby achieving adaptive threshold screening;

[0021] In the normal distribution statistics of the travel time residual data, data outside of 2 times the standard deviation are removed, and the normal distribution density function formula is:

[0022]

[0023] Where r is the travel time residual value, f(r) is the travel time residual normal distribution density function, μ is the mean, σ is the standard deviation, π is pi, and e is a natural constant.

[0024] Preferably, the global inversion algorithm used in step A2 is a grid search method, which assumes that the seismic station is located at each grid point, calculates the epicentral distance and back azimuth between each grid point and the earthquake event, and then compares them with the known information in the seismic phase observation report. The grid point with the smallest difference is the accurate seismic station location in the case of a single earthquake event;

[0025] The grid search method uses a variable grid strategy of "coarse grid first, then fine grid". It first performs a global search on a grid covering a larger range, gradually densifies the grid after obtaining the optimal position, and then performs a local search near the optimal position. After multiple rounds of iterations, a more accurate seismic station location is finally obtained.

[0026] Preferably, the calculation formula for the epicenter distance and the back azimuth is:

[0027] Δ=arccos(sin(φ grid ′)cos(λ grid )*sin(φ evt ′)cos(λ evt )+sin(φ grid ′)sin(λ grid )*sin(φ evt ′)sin(λ evt )+cos(φ grid ′)cos(φ evt ′)) (5)

[0028]

[0029]

[0030] where Δ and Bazi represent the epicentral distance and back azimuth between the grid point and the earthquake event location, respectively, and φ evt and λ evt represent the longitude and latitude of the earthquake event, φ grid and λ grid They represent the longitude and latitude of the grid point respectively, fla is the flattening of the earth, and π is pi.

[0031] Preferably, in step A2, the epicentral distance and back azimuth between each grid point and the latitude and longitude position of one earthquake event are calculated, and then compared with the known information, and the grid point with the smallest difference is selected as the optimal position of this round of iteration. The calculation formula is:

[0032]

[0033] Among them, R i is the gap value at the i-th grid point, and denote the back azimuth of the ith grid point and the actual seismic station relative to the earthquake event position, respectively. and They represent the epicenter distance of the i-th grid point and the actual seismic station relative to the earthquake event location, respectively.

[0034] Preferably, in step A2, since the positions of the same seismic station obtained from different earthquake events are different, a weighted average is calculated for each of the multiple longitude and latitude values ​​obtained, with the weight determined by the epicenter distance, to obtain the accurate position of the seismic station under the calculation of multiple earthquake events. The calculation formula is:

[0035]

[0036] in, The latitude and longitude position of the station calculated for the i-th earthquake event, w i is the weight of the i-th earthquake event, is the final latitude and longitude position of the seismic station, Dist max and Dist min are the maximum and minimum epicenter distances of all earthquake events, Dist i is the epicenter distance of the ith earthquake event, and n is the number of earthquake events.

[0037] Preferably, in step A3, the elevation values ​​between the grid points are obtained using a bilinear interpolation method, and the calculation formula is:

[0038]

[0039] Among them, H * (x, y) is the station elevation value at the (x, y) position after interpolation, H(x i ,y j ) are the four grid points around (x,y).

[0040] Preferably, in step A4, the three-dimensional spatial distance between the known station position and the corresponding station position obtained in step A2 and step A3 is calculated, and the error level is obtained by using the average value statistical method, so as to quantitatively estimate the accuracy of the final station position result.

[0041] The beneficial effects of the present invention are:

[0042] Through the operation of preprocessing the seismic phase observation report data, the residual screening threshold is adaptively adjusted according to the statistical relationship under different data discrete conditions, ensuring the accuracy and richness of the input seismic phase report data.

[0043] By adopting the "variable grid division + multi-round iterative search" strategy in calculating the longitude and latitude information of seismic stations, while achieving high-precision global search results, it also further reduces the calculation time, ensuring high precision and high efficiency in the process of obtaining the seismic station location.

[0044] By adopting the method of seismic station location information accuracy analysis, the final seismic station location error level can be quantitatively estimated and qualitatively judged, achieving clear control of the quality of the calculation results, which will help provide a quality control standard reference for subsequent related research using the data. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 It is a travel time residual distribution diagram and a screening range schematic diagram of the present invention;

[0046] Figure 2 is a schematic diagram of an example of grid point setting of the present invention;

[0047] Figure 3 It is a schematic flow chart of the method of the present invention;

[0048] Figure 4 This is a comparison chart of the gap between the station location information obtained in Example 2 and the known station information. DETAILED DESCRIPTION

[0049] The technical solution of the present invention is further described below, but the scope of protection claimed is not limited to the description.

[0050] Example 1:

[0051] like Figures 1 to 3 As shown, a method for calculating the location of a seismic station based on a seismic phase observation report includes the following steps:

[0052] Step 1: Preprocessing of seismic phase observation report data

[0053] This implementation uses the 2024 China Earthquake Networks seismic phase observation report data within the longitude and latitude range of 102°-106°E, 26°-30°N as an example. This data includes both 25 stations of the National Seismic Network whose station locations have been made public, and 93 stations of the provincial seismic networks that have not yet been made public.

[0054] Extract the phase arrival time and the corresponding earthquake occurrence time and epicenter distance data from the phase observation report, and convert them into the format of "phase travel time-epicenter distance", where the phase travel time is obtained by the difference between the phase arrival time and the earthquake occurrence time. This implementation selects the Pg phase for extraction and conversion because the Pg phase is usually the clearest and most reliable in the earthquake analysis process. After the format conversion, the above data is fitted with the least squares method using the linear fitting formula to obtain the theoretical time-distance curve. The fitting formula is as follows:

[0055]

[0056] t=ad+b (3)

[0057] Among them, d is the independent variable of the epicenter distance of the theoretical time-distance curve, t is the dependent variable of the earthquake phase travel time of the theoretical time-distance curve, a is the slope obtained by fitting, b is the intercept obtained by fitting, and x i is the epicenter distance of the earthquake event corresponding to the i-th phase in the data involved in the fitting, y i is the travel time of the ith phase, and n is the number of phases. After fitting, the slope a = 0.171 and the intercept b = 0.129 are obtained, so the theoretical time-distance curve formula is t = 0.171d + 0.129;

[0058] Input the epicenter distance data corresponding to each earthquake in the phase observation report into the theoretical time-distance curve to obtain the theoretical phase travel time, and then subtract it from the actual observed phase travel time extracted from the report to obtain the travel time residual data corresponding to each phase. Generally speaking, the travel time residual data satisfies the normal distribution with 0 as the mean. Data with excessive residuals may indicate a larger analysis error and need to be removed. Therefore, this implementation performs normal distribution statistics on the travel time residual data (such as Figure 2 ), remove the data beyond 2 times the standard deviation (about 4.6% of the total data), and the normal distribution density function is as follows:

[0059]

[0060] Where r is the travel time residual, f(r) is the normal distribution density function of the travel time residual, μ is the mean, σ is the standard deviation, π is the circumference of the circle, and e is a natural constant. After statistical analysis, the mean of the travel time residuals used in this implementation is 0, the standard deviation is 0.859, and 2 times the standard deviation is 1.717. This means that phase data with travel time residuals greater than 1.717 will be removed.

[0061] To ensure the credibility of input data, this implementation further restricts the minimum number of phases and maximum station gap angle recorded by each participating station. The minimum number of phases represents the station's ability to record earthquake events and reflects the quality of the station data. The more phases received, the better the station data quality. This implementation is set to 20. The maximum station gap angle represents the azimuth coverage of earthquake events recorded by the station and reflects the representativeness of the station data. The smaller the maximum station gap angle, the better the azimuth coverage and the more representative the station data. This implementation is set to 300°.

[0062] Step 2: Calculate the longitude and latitude information of the seismic station

[0063] This implementation uses a grid search method to obtain the longitude and latitude information of seismic stations. First, the longitude and latitude positions, epicenter distance, and back azimuth information of all earthquake events corresponding to each seismic station are extracted. Then, a search grid is divided within the study area, and each grid point is regarded as a possible location of a station. The epicenter distance and back azimuth between each grid point and each earthquake event are then calculated. Choosing different calculation methods will lead to different final results. The calculation formula used in this implementation is as follows:

[0064] △=arccos(sin(φ grid ′)cos(λ grid )*sin(φ evt ′)cos(λ evt )+sin(φ grid ′)sin(λ grid )*sin(φ evt ′)sin(λ evt )+cos(φ grid ′)cos(φ evt ′)) (5)

[0065]

[0066]

[0067] Where Δ and Bazi represent the epicentral distance and back azimuth between the grid point and the earthquake event location, respectively, and φ evt ,λ evt and φ grid ,λ grid Represent the longitude and latitude of the earthquake event and grid point respectively, and fla is the Earth's flattening. It should be noted that the density of the grid division directly determines the accuracy of the final seismic station location. That is, the smaller the grid density, the more accurate the final result, but the lower the calculation efficiency. The larger the grid density, the higher the calculation efficiency, but at the expense of the accuracy of the final result. Therefore, this implementation uses the strategy of "variable grid division + multiple rounds of iterative search" to calculate the station location. The specific steps are as follows:

[0068] 1) Determine the grid division mode according to the accuracy requirements. This implementation adopts the "1°→0.1°→0.01°→0.001°" gradually denser division mode for grid search;

[0069] 2) Perform multiple rounds of iterative search based on the divided grids. First, perform a global search on the grid points with a spacing of 1°×1°. Calculate the epicenter distance and back azimuth between each grid point and the latitude and longitude position of one of the earthquake events. Then compare it with the known information and select the grid point with the smallest difference as the optimal position for this round of iteration. Figure 3As shown in the figure, different shapes represent the search grid point positions at different grid spacings, and the shape color represents the difference between the epicenter distance and back azimuth information calculated at each grid point and the actual station information. The darker the color, the smaller the difference. The five-pointed star represents the calculated optimal station position. This implementation uses the minimum difference value represented by the sum of squared differences R as the optimal judgment criterion. The calculation formula is as follows:

[0070]

[0071] Among them, R i is the gap value at the i-th grid point, and denote the back azimuth of the ith grid point and the actual seismic station relative to the earthquake event position, respectively. and represent the epicentral distance of the ith grid point and the actual seismic station relative to the earthquake event location, respectively;

[0072] After the optimal grid point is determined in the previous round of iterations, a more dense and fine grid is divided with this point as the center and the previous round of grid spacing as the maximum edge distance to perform local grid search. That is, the grid spacing in the second round is 0.1°×0.1°, and the maximum edge distance from the grid center point is 1°. The grid spacing in the third round is 0.01°×0.01°, and the maximum edge distance from the grid center point is 0.1°. This process is repeated until the minimum grid spacing is reached, which is set to 0.001° in this implementation.

[0073] (3) The above calculated seismic station location is only obtained through the known information of one earthquake event. However, a seismic station can usually receive multiple earthquake events. Therefore, the latitude and longitude positions of the same station calculated by these earthquake events are weighted averaged to obtain the final latitude and longitude position of the seismic station:

[0074]

[0075] in, The latitude and longitude position of the station calculated for the i-th earthquake event, w i is the weight of the i-th earthquake event, is the final latitude and longitude position of the seismic station, Dist max and Dist min are the maximum and minimum epicenter distances of all earthquake events, Dist i is the epicenter distance of the ith earthquake event, and n is the number of earthquake events;

[0076] Step 3: Extract seismic station elevation information

[0077] This implementation uses global terrain relief data from SRTM15+V2.6. To improve computational efficiency and reduce memory usage, the global terrain relief data is clipped to the target area (102°-106°E, 26°-30°N) before extracting elevation information. Elevation information is then extracted at each seismic station.

[0078] It should be noted that the grid search spacing in this implementation is 0.001°, while the grid size of the SRTM15+V2.6 data is 15s (approximately 0.00417°). This results in some stations not being on the grid points of the terrain relief data and making it impossible to directly extract data. Therefore, this implementation uses the bilinear interpolation method to obtain the elevation information between the grid points. The calculation formula is as follows:

[0079]

[0080] Among them, H * (x, y) is the station elevation value at the (x, y) position after interpolation, H(x i ,y j ) are the four grid points around (x,y).

[0081] Step 4: Accuracy analysis of seismic station location information

[0082] This implementation obtained the latitude, longitude, and elevation information of all seismic stations in the study area, including some stations whose locations have been made public. The latitude, longitude, and elevation information of these known stations were compared with the corresponding information obtained in this implementation. Since the latitude, longitude, and elevation information are independent in the process of calculation, the accuracy analysis needs to be performed separately. The quantitative calculation formula for the comparison difference is as follows:

[0083]

[0084] Among them, R L and R H They represent the latitude and longitude difference and elevation difference of the stations respectively, x i ,y i ,H i represent the longitude, latitude and elevation of the i-th station obtained in this implementation, respectively, and x i ',y i ',y i 'represent the longitude, latitude and elevation values ​​of the i-th known position station, and n is the number of known position stations.

[0085] Based on the obtained gap values, this implementation defines the following criteria for quantitative accuracy assessment:

[0086] Table 1 Accuracy judgment criteria

[0087] Gap value Accuracy rating <![CDATA[R L <0.005° and R H <100m]]> S <![CDATA[0.005°<R L ≤0.01° or 100m <R H <200m]]> A <![CDATA[0.01°<R L ≤0.02° or 200m <R H <300m]]> B <![CDATA[0.02°<R L ≤0.05° or 300m <R H <500m]]> C <![CDATA[R L ≥0.05° and R H ≥500m]]> D

[0088] Example 2

[0089] To further demonstrate the accuracy and effectiveness of the present invention, this paper uses the 2024 China Earthquake Networks seismic phase observation report data within the latitude and longitude range of 102°-106°E, 26°-30°N, calculates the locations of all relevant seismic stations, and extracts the data of publicly available stations for comparison ( Figure 4 Among the 25 compared seismic stations, the average latitude and longitude difference is 0.0093°, and the average elevation is 46.7 m, which meets the A-level standard in Table 1. This shows that the results obtained in this implementation are relatively accurate and can be used for further research.

Claims

1. A method for calculating the location of a seismic station based on seismic phase observation reports, characterized in that: The following steps are involved: A1: Preprocessing of seismic phase observation report data: Set a longitude and latitude range, collect earthquake event locations, multiple epicenter distance and back azimuth data, and surface undulation grid data within the longitude and latitude range, and record known and unknown seismic stations within the longitude and latitude range; A2: Calculate the longitude and latitude information of seismic stations: By sorting and extracting the longitude and latitude of earthquake events recorded by each seismic station in the seismic phase observation report and the corresponding epicenter distance and back azimuth parameters, a global inversion algorithm is used to obtain the longitude and latitude of each seismic station calculated for each earthquake event. Finally, the longitude and latitude information of the same seismic station obtained from multiple earthquake events is merged and sorted to obtain the accurate longitude and latitude information of each seismic station. A3: Extracting seismic station elevation information: Based on the calculated seismic station longitude and latitude information, extract elevation information from the corresponding location in the surface relief grid data. If the station longitude and latitude are not on the grid point, interpolate the elevation information from the four nearby grid points. A4: Accuracy analysis of seismic station location information: Based on the seismic stations whose location information has been disclosed in the seismic phase observation report, the calculation results are verified and analyzed for accuracy, so as to quantitatively determine the accuracy of the final calculation results.

2. A method for calculating the location of a seismic station based on a seismic phase observation report according to claim 1, characterized in that: In step A1, the epicentral distance and phase arrival time data are extracted from the phase observation report, and the phase travel time is calculated by the difference between the phase arrival time and the earthquake occurrence time. Then, the theoretical time distance curve is obtained by linear fitting, and different theoretical phase travel times are obtained according to different epicentral distances. Subsequently, the actual travel time data of the phase observation report is subtracted from the corresponding theoretical phase travel time data to obtain travel time residual data. Finally, a suitable threshold is selected to remove data with excessive residuals, thereby completing the data preprocessing work.

3. A method for calculating the location of a seismic station based on a seismic phase observation report according to claim 2, characterized in that: The theoretical time-distance curve is obtained by performing least square fitting on the data through the linear fitting formula. The fitting formula is: t=ad+b (3) Among them, d is the independent variable of the epicenter distance of the theoretical time-distance curve, t is the dependent variable of the earthquake phase travel time of the theoretical time-distance curve, a is the slope obtained by fitting, b is the intercept obtained by fitting, and x i is the epicenter distance of the earthquake event corresponding to the i-th phase in the data involved in the fitting, y i is the travel time of the ith phase, and n is the number of phases.

4. The method for calculating the location of a seismic station based on a seismic phase observation report according to claim 2, wherein: In step A1, all the obtained travel time residual data are subjected to normal distribution statistics, with P wave travel time data being prioritized to obtain their mean and standard deviation, and 2 times the standard deviation is selected as the data screening threshold, thereby achieving adaptive threshold screening; In the normal distribution statistics of the travel time residual data, data outside the double standard deviation are removed. The normal distribution density function formula is: Where r is the travel time residual value, f(r) is the travel time residual normal distribution density function, μ is the mean, σ is the standard deviation, π is pi, and e is a natural constant.

5. The method for calculating the location of a seismic station based on a seismic phase observation report according to claim 1, wherein: The global inversion algorithm used in step A2 is a grid search method. By assuming that the seismic station is located at each grid point, the epicentral distance and back azimuth between each grid point and the earthquake event are calculated respectively. Then, the distance is compared with the known information in the seismic phase observation report. The grid point with the smallest difference is the accurate seismic station location in the case of a single earthquake event. The grid search method uses a variable grid strategy, first performing a global search on a grid covering a larger range, gradually encrypting the grid after obtaining the optimal position, and then performing a local search near the optimal position. After multiple rounds of iterations, a more accurate seismic station location is finally obtained.

6. A method for calculating the location of a seismic station based on seismic phase observation reports according to claim 5, characterized in that: The calculation formula of the epicenter distance and back azimuth is: Δ=arccos(sin(φ grid ')cos(λ grid )*sin(φ evt ')cos(λ evt ) +sin(φ grid ')sin(λ grid )*sin(φ evt ')sin(λ evt ) +cos(φ grid )')cos(φ evt ')) (5) where Δ and Bazi represent the epicentral distance and back azimuth between the grid point and the earthquake event location, respectively, and φ evt and λ evt represent the longitude and latitude of the earthquake event, φ grid and λ grid They represent the longitude and latitude of the grid point respectively, fla is the flattening of the earth, and π is pi.

7. The method for calculating the location of a seismic station based on seismic phase observation reports according to claim 5, characterized in that: In step A2, the epicentral distance and back azimuth between each grid point and the latitude and longitude position of one earthquake event are calculated, and then compared with the known information. The grid point with the smallest difference is selected as the optimal position for this round of iteration. The calculation formula is: Among them, R i is the gap value at the i-th grid point, and denote the back azimuth of the ith grid point and the actual seismic station relative to the earthquake event position, and They represent the epicentral distance of the i-th grid point and the actual seismic station relative to the earthquake event location, respectively.

8. The method for calculating the location of a seismic station based on seismic phase observation reports according to claim 5, characterized in that: In step A2, since the locations of the same seismic station obtained for different earthquake events are different, a weighted average is calculated for each of the multiple longitude and latitude values ​​obtained, with the weight determined by the epicenter distance, to obtain the accurate location of the seismic station under the calculation of multiple earthquake events. The calculation formula is: in, The latitude and longitude position of the station calculated for the i-th earthquake event, w i is the weight of the i-th earthquake event, is the final latitude and longitude position of the seismic station, Dist max and Dist min are the maximum and minimum epicenter distances of all earthquake events, Dist i is the epicenter distance of the ith earthquake event, and n is the number of earthquake events.

9. The method for calculating the location of a seismic station based on seismic phase observation reports according to claim 1, characterized in that: In step A3, the elevation values ​​between grid points are obtained using a bilinear interpolation method, and the calculation formula is: Among them, H * (x, y) is the station elevation value at the (x, y) position after interpolation, H(x i ,y j ) are the four grid points around (x,y).

10. The method for calculating the location of a seismic station based on seismic phase observation reports according to claim 1, characterized in that: In step A4, the three-dimensional spatial distance between the known station position and the corresponding station position obtained in steps A2 and A3 is calculated, and the error level is obtained using average value statistics, thereby quantitatively estimating the accuracy of the final station position result.

Images