A focal mechanism inversion method based on correlation degree fitting of polarity and amplitude ratio

Abstract

The application provides a focal mechanism inversion method based on polarity and amplitude ratio correlation fitting, and belongs to the field of earthquake correlation. The method firstly quantitatively analyzes systematic deviation of regional medium difference and path attenuation to actual amplitude ratio data, and then jointly inverts the focal mechanism based on polarity and amplitude ratio data by using the HCN algorithm. The algorithm generates multiple data sets by randomly adding error disturbances to original data, and performs initial screening on selectable solutions in the whole space of focal mechanism parameters based on polarity data of P, SV and SH waves. The result with the highest correlation degree of theoretical and actual amplitude ratio curves is selected by secondary screening. Finally, the average tensor of the secondary screening results of multiple error disturbance data sets is calculated to determine the final focal mechanism inversion result, and the fitting degree and dispersion of the result with the actual observation data are evaluated. The application can effectively deal with the influence of regional medium difference and path attenuation on amplitude ratio data, improve the accuracy of focal mechanism inversion, and provide a reliable basis for earthquake research.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of earthquake correlation, and more particularly relates to a focal mechanism inversion method based on polarity and amplitude ratio correlation fitting. BACKGROUND

[0002] The focal mechanism of an earthquake is key information for revealing the regional stress state, fault activity characteristics and the mechanism of earthquake origin, and the inversion accuracy directly affects the understanding of the above problems. For moderate earthquakes, methods such as gCAP based on waveform fitting can better invert the focal mechanism thereof; but for small earthquakes (e.g. M≤3.5), waveform fitting inversion methods are difficult to give reliable results due to factors such as low signal-to-noise ratio of waveforms, non-optimized effective signal frequency band, etc. However, small earthquakes have high occurrence frequency, wide spatial distribution, and are sensitive to local stress field changes, and often have good potential in the fine study of the above deep-seated geoscience problems (Hu Xingping et al., 2021), therefore, it is very important to develop a stable and reliable small earthquake focal mechanism inversion method. The constraint based on polarity and amplitude ratio data is an important way to invert the focal mechanism, especially the focal mechanism of small earthquakes; compared with the inversion method based on polarity alone, such as FPFIT, the introduction of amplitude ratio data can make up for the lack of quantity and distribution of polarity data, and provide more constraints for focal mechanism inversion, and so far it is still an important way for many scholars to study the focal mechanism of small earthquakes. At present, the classic algorithms for inverting the focal mechanism using polarity and amplitude ratio data include FOCMEC using the polarity of P wave, SV wave and SH wave and their amplitude ratio, and HASH and its Python improved version SKHASH using the polarity of P wave and the amplitude ratio of PS.

[0003] However, the actual observed seismic wave amplitude is largely affected by factors such as regional medium differences and path attenuation, and there is a systematic deviation; if not corrected directly, its reliability cannot be guaranteed, and even academic controversy has arisen as to whether the amplitude ratio is beneficial in the practical application of focal mechanism inversion.

[0004] In response to this, subsequent scholars have proposed different improved algorithms and techniques from different ideas, including relative focal mechanism inversion algorithms such as relative radiation pattern joint inversion (REFOC) of adjacent earthquakes, relative focal mechanism inversion (RFMI) of adjacent small earthquakes using reliable main earthquake source mechanism, and amplitude ratio correction techniques such as using the average theory and observed amplitude ratio of DAS dense observation to correct the amplitude ratio, using the linear transformation given by the analysis of a large amount of amplitude ratio data to correct the amplitude ratio. However, the above improved algorithms and techniques also have certain limitations. For example, the REFOC algorithm relies on the clustered distribution of earthquakes, and is difficult to apply to scattered earthquakes; the RFMI algorithm truly solves the deviation of small earthquakes relative to the main event, and the reliability of the absolute mechanism of small earthquakes needs to be further clarified; at the same time, the calculation cost of relative algorithms such as REFOC and RFMI is relatively high, which may be restricted by computing power in large-scale application. The use of average values of dense observations such as DAS to correct amplitude ratio reduces the applicability when observation stations are dispersed; the correction method using linear transformation given by analysis of a large amount of amplitude ratio data still needs to be verified when applied in different regions.

[0005] It can be said that the previous improved algorithms and techniques have many limitations such as requiring more prior information, requiring high computing power, and having low universality when dealing with systematic deviations of seismic wave amplitude caused by regional medium differences and path attenuation and other factors, which makes it difficult to effectively support the research and development of small earthquake source mechanism and other related earth science problems. SUMMARY

[0006] In response to this, the present application starts from the basic theory of quantitative seismology, quantitatively analyzes the systematic deviations of amplitude ratio data caused by regional medium differences and path attenuation and other factors, and develops an HCN algorithm for inverting the source mechanism using the polarity and amplitude ratio data of P, SV and SH waves based on MATLAB language. In the process of source mechanism inversion, a new amplitude ratio data fitting method is adopted, that is, the correlation degree fitting of the amplitude ratio change curve between different stations of the theoretical and actual observation, which can eliminate the systematic deviations caused by regional medium differences and propagation path attenuation in principle, so as to maximize the reliability and accuracy of the source mechanism inversion results.

[0007] In order to achieve the above purpose, the present application is realized by adopting the following technical solutions: the method comprises:

[0008] According to the theory of quantitative seismology, the systematic deviations of actual amplitude ratio data caused by regional medium differences and path attenuation are quantitatively analyzed;

[0009] The HCN algorithm is used to jointly invert the focal mechanism based on the polarity and amplitude ratio data of P, SV and SH waves. Firstly, random error disturbance is added to the original data to generate multiple data sets. The initial screening of the selectable solutions is based on the polarity data of P, SV and SH waves in the full space of the focal mechanism parameters. Then, the best fitting results of the theoretical and actual observed amplitude ratio curves are selected as the selectable solutions for the secondary screening. The union of the secondary screening results of multiple random error disturbances is taken to determine the final focal mechanism inversion result by calculating the average tensor. The quality of the final focal mechanism inversion result is comprehensively evaluated according to the polarity weighted fitting degree, the amplitude ratio curve correlation coefficient and the angular dispersion of the selectable solution space of the actual observed data and the result.

[0010] In one scheme, the system bias quantitative analysis is specifically:

[0011] According to the theory of quantitative seismology, the principle of amplitude ratio inversion of focal mechanism is that the radiation patterns of P and S waves in the wave displacement field generated by a double force couple point source are different, which can be expressed as:

[0012] (1)

[0013] Considering that the seismic wave will be attenuated related to the propagation path from the source to the station, it can be expressed as wherein, 、 、 represent the frequency, travel time and attenuation factor of the seismic wave respectively. In addition, considering the possible regional differences of the actual PS wave velocity ratio, the system bias caused by the two can be quantitatively expressed as:

[0014] (2)

[0015] If the commonly used logarithmic value (log value) in amplitude ratio fitting is used, formula (2) can be written as:

[0016] (3)

[0017] In the above formula, A P 、A SV 、A SH represent the amplitudes of P, SV and SH waves respectively, α、β represent the propagation velocities of P and S waves respectively, θ is the angle between the ray and the normal of the fault plane, is the angle between the projection of the ray in the fault plane and the fault slip direction, and the upper left superscript theoIndicates theoretical value, superscript. prac The actual observed value is represented by formulas (2) and (3). Formulas (2) and (3) quantitatively reflect that regional medium differences and path attenuation factors will cause systematic deviations between theoretical and actual amplitude ratio data.

[0018] In one approach, the HCN algorithm inversion process includes:

[0019] Step 1: Use polarity data to weighted fit and initially screen for possible solutions;

[0020] Step 2: From the initial screening of possible solutions, a secondary screening is conducted to select the results with the highest correlation to the amplitude ratio curve;

[0021] Step 3: Calculate the averaged tensor final solution and evaluate its quality.

[0022] In one approach, the initial screening of possible solutions using polarity data weighted fitting includes:

[0023] Set the disturbance errors at the source angle and azimuth angle as err ih =5° and err az =2°, the distance from the source and azimuth of the original observation data are respectively within ±err ih and ±err az Do within the scope N D Ten random perturbations are performed (default is 10 times) to generate the inversion (default is no perturbation for the first time). N D The perturbation data samples were analyzed. Using a spatial angle step of 5°, the focal mechanism parameters (fault strike 0°~355°, fault dip 0°~90°, fault slip angle 0°~355°) were divided into 72×19×72=98496 parameter combinations. For each focal mechanism parameter combination, the theoretical polarity of each station was calculated and compared with the actual observed polarity data to select the weighted polarity fit (…). PF Not less than the maximum polarity fit among all parameter combinations ( PF M Subtract the downward floating threshold ( DPF The parameter combination is used as the initial screening of possible solutions, i.e. PF ≥ PF M - RPF Due to the existence N D Given a sample of error-perturbed data, the initial screening solution could be: N D The union of the solutions from the initial screening.

[0024] In one approach, the secondary selection of the result with the highest correlation to the amplitude ratio curve from the initially screened alternatives specifically includes:

[0025] According to the polarity data fitting preliminary screening can be selected solution in each observation station (N stations, STA1..., STA N The three-term overall curve composed of the three theoretical amplitude ratios (A1, A2, A3) and the three curves composed of the actual observation are calculated for correlation degree, and the result with the highest comprehensive correlation degree is selected as the secondary screening result. A P / A SV , A P / A SH , A SV / A SH The weights of the three curves are equal, or different weights are set according to the actual amplitude reliability evaluation; for N D error disturbance data samples, the union of the N D independent calculation and the secondary screening of the selected solutions is selected as the result of this stage.

[0026] In one scheme, the calculation of the average tensor final solution and the evaluation of its quality include:

[0027] For N D random error disturbances, the union of the selected solutions is calculated, and the average tensor of these solutions is determined, so as to determine the final focal mechanism inversion result; then, the polarity weighted fitting degree of the inversion result and the actual observation data, the amplitude ratio curve correlation coefficient, and the angular dispersion of the selected solutions are calculated, so as to comprehensively evaluate the quality of the final focal mechanism inversion result.

[0028] The present application has the following advantages:

[0029] 1. The previous amplitude ratio inversion algorithm for focal mechanism all uses amplitude ratio residual statistics based screening to select the selected solution, compared with which, the HCN algorithm fitting the relative height of amplitude ratio between different stations is more theoretically reasonable.

[0030] 2. The correlation fitting used by the HCN algorithm has the characteristic of linear transformation invariance, which can eliminate the systematic deviation of regional medium wave velocity ratio difference and path attenuation in principle.

[0031] 3. The HCN algorithm considers the error disturbance of polarity picking and ray path, which can reduce the probability of losing the "true" solution, and can also improve the reliability of the final focal mechanism inversion result. At the same time, the step-by-step screening method is adopted for polarity and amplitude ratio data instead of comprehensive screening, which is more reasonable and stable.

[0032] 4、HCN algorithm has priori information, low algorithmic cost, strong universality, good convenience and many other technological advantages and practical advantages, and shows excellent inversion effect in different actual scenes, and provides a general, simple and efficient technical approach for accurate determination of focal mechanism solution in multiple scenes. BRIEF DESCRIPTION OF DRAWINGS

[0033] Figure 1 The method technical framework diagram of the application is shown in the figure;

[0034] Figure 2 The actual data inversion step example of the method of the application is shown in the figure;

[0035] Figure 3 The inversion results (a) in different scenes such as teleseism, regional earthquake, small and micro mine earthquake are shown in the figure;

[0036] Figure 4 The inversion results (b) in different scenes such as teleseism, regional earthquake, small and micro mine earthquake are shown in the figure. DETAILED DESCRIPTION

[0037] In order to facilitate the understanding of the present application, the present application will be described more fully below with reference to the accompanying drawings. The drawings show typical embodiments of the present application. However, the present application can be implemented in many different forms and is not limited to the embodiments described in the present application. On the contrary, the purpose of providing these embodiments is to make the disclosure of the present application more thorough and comprehensive.

[0038] Unless otherwise defined, all technical and scientific terms used in the present application have the same meaning as understood by those skilled in the art to which the present application belongs. The terms used in the present application in the specification of the present application are only for the purpose of describing the specific embodiments and are not intended to limit the present application. In order to facilitate the understanding of the present application, the present application will be described more fully below with reference to the accompanying drawings. The drawings show typical embodiments of the present application. However, the present application can be implemented in many different forms and is not limited to the embodiments described in the present application. On the contrary, the purpose of providing these embodiments is to make the disclosure of the present application more thorough and comprehensive.

[0039] The principle of the present application: the HCN algorithm first generates multiple data sets by randomly disturbing the original data, respectively fits the initial screening selectable solutions in the whole space of focal mechanism parameters based on the polarity data of P wave, SV wave and SH wave, then selects the best fitting result of the theoretical and actual observation amplitude ratio change curve as the selectable solution for secondary screening; the union of the secondary screening results of multiple random error disturbances is determined by calculating the average tensor to determine the final focal mechanism inversion result, and the quality of the final focal mechanism inversion result is comprehensively evaluated according to the polarity weighted fitting degree of the result and the actual observation data, the amplitude ratio curve correlation coefficient and the angular dispersion of the selectable solution space.

[0040] In addition to effectively dealing with regional medium differences and propagation path attenuation, the HCN algorithm can improve the accuracy of focal mechanism inversion. It has many technological and practical advantages, such as less prior information, low computational cost, strong universality, and good convenience, which can provide efficient and reliable technological support for the study of earthquake focal mechanisms and other related earth science problems.

[0041] As shown in Figure 1 , a focal mechanism inversion method based on polarity and amplitude ratio correlation fitting includes:

[0042] The principle of using amplitude ratio to invert focal mechanism is that the radiation patterns of P and S waves in the wave displacement field generated by a double force couple point source are different. By fitting the waveform amplitude ratio recorded by different stations, the spatial angle relationship between the exit point of the seismic ray backtracking to the focal sphere and the fault plane normal and slip direction can be inferred, so as to determine the focal mechanism. According to the theory of quantitative seismology, this principle can be expressed as:

[0043] (1)

[0044] In formula (1), A P 、A SV 、A SH A, B, and C represent the amplitudes of P, SV, and SH waves, respectively. α、β Vp and Vs represent the propagation velocities of P and S waves, respectively. θ is the angle between the ray and the fault plane normal. is the angle between the projection of the ray in the fault plane and the fault slip direction, and the upper left superscript theo represents the theoretical value.

[0045] Based on the basic principle of using polarity and amplitude ratio to invert focal mechanism expressed in formula (1), considering that seismic waves will attenuate related to the propagation path from the source to the station, it can be expressed as where, 、 、 A, B, and C represent the amplitudes of P, SV, and SH waves, respectively.

[0046] (2)

[0047] If the commonly used logarithmic value (log value) in amplitude ratio fitting is used, formula (2) can be written as:

[0048] (3)

[0049] The upper left subscript in formula (2), (3) prac represents the actual observation value, which quantitatively reflects that the regional medium difference and path attenuation and other factors will cause systematic deviation between the theoretical and actual amplitude ratio data. How to reasonably and effectively eliminate the influence of these error factors is the key to obtaining accurate and reliable focal mechanism. The previous inversion methods and technical improvements all have certain limitations.

[0050] The HCN algorithm uses the polarity and amplitude ratio data of P waves, SV waves and SH waves to jointly invert the focal mechanism, so as to maximize the utilization of polarity and amplitude ratio data. The steps include: ① using polarity data to weightedly fit the initial screening selectable solution; ② from the initial screening selectable solution, the result with the highest correlation degree of amplitude ratio curve is selected again; ③ the average tensor final solution is calculated and its quality is evaluated. Among them, step ② is an innovative solution proposed by the HCN algorithm for the above-mentioned systematic deviation, that is, a brand-new amplitude ratio data fitting method is adopted, that is, the residual fitting of the actual observation amplitude ratio and the theoretical value of different stations in the previous algorithm is changed into the correlation fitting of the observation amplitude ratio change curve and the theoretical curve between different stations; according to the standard of amplitude ratio screening focal mechanism, the residual statistics of actual observation amplitude ratio and theoretical value is changed from minimum to the highest correlation degree of amplitude ratio change curve and theoretical curve between different stations, so as to avoid the influence of abnormal observation to the maximum extent.

[0051] Figure 2 Taking the actual polarity and amplitude ratio data of the Kamchatka 8.7 earthquake event on July 30, 2025 (Beijing time) as an example, the focal mechanism inversion process of the HCN algorithm is introduced. The polarity and amplitude ratio data obtained by manual picking of original seismic data belong to the data preprocessing process, or can be completed by means of seismic signal picking technology. For the focal mechanism inversion process of the HCN algorithm, the polarity and amplitude ratio data obtained are used as the starting point.

[0052] Step 1, using polarity data to weightedly fit the initial screening selectable solution

[0053] The HCN algorithm considers that the errors of earthquake positioning, regional velocity structure and other factors will bring uncertainty to the position of the observation station backtracking to the focal sphere (reflected in the input data of off-source angle and azimuth), and sets the perturbation errors of off-source angle and azimuth as err ih =5° and err az =2°, respectively. The off-source angle and azimuth of the original observation data are respectively within ±err ih and ±err az range N DThe number of random perturbation inversions (default 10 times) (default 1st time without perturbation) generates N D error range is set based on previous research and applicant's previous research understanding.

[0054] According to a spatial angle step of 5°, the focal mechanism parameter full space (strike of fault plane 0°~355°, dip of fault plane 0°~90°, and fault slip angle 0°~355°) is divided into 72×19×72=98496 parameter combinations; this 5° step is a value balancing focal mechanism inversion accuracy and computing power requirements. For each focal mechanism parameter combination, the theoretical polarity of each station is calculated and compared with the actual observed polarity data, and the parameter combination with a weighted polarity fitting degree (WPF) not less than the maximum polarity fitting degree (MPF) of all parameter combinations minus a floating amplitude threshold (FT) (MPF-FT) is selected as a preliminary screening solution, i.e. PF ≥ PF M DPF PF ≥ PF M - RPF Since there are N D error perturbation data samples, the preliminary screening solution is the union of N D times of preliminary screening solutions.

[0055] The HCN algorithm sets DPF to 10%, because statistics show that the initial motion polarity picking error rate reaches 10%~20%. The polarity data weight in the weighted polarity fitting degree includes two layers of weighting. One is that the user can set the initial data weight according to the polarity picking quality in advance, i.e. assigning a value between -1 and 1, ± indicating the polarity positive and negative, and the absolute value being the set initial weight; the other is that the HCN algorithm dynamically adjusts the weight on the basis of the initial weight according to the P, SH, and SV radiation patterns, for example, the P wave theoretical amplitude near the nodal plane is very small, and its P wave polarity weight is reduced. The reason for adopting this weighting method is that the polarity with small theoretical amplitude is more likely to be misjudged under noise interference, and the position of station observation data back to the focal sphere has a certain error; adopting dynamic weight can avoid the excessive control of polarity data on the focal mechanism nodal plane and the interference with the reasonable evaluation of the final result polarity fitting degree.

[0056] Figure 2 ​​① The distribution of polarity data and the initial screening results of focal mechanism on focal sphere: + and - represent the picking polarity of P wave in vertical direction, positive (upward) and negative (downward) respectively; F and B represent the picking polarity of SV wave in radial direction, positive (away from the source) and negative (towards the source) respectively; R and L represent the picking polarity of SH wave in tangential direction, positive (right, as seen from the earthquake to the station) and negative (left, as seen from the earthquake to the station) respectively; the yellow thin line represents the optional solution obtained based on polarity weighted fitting initial screening. It can be seen that due to the number and distribution of polarity data on the focal sphere, the constraint of the optional solution is poor, and the initial screening results are very discrete in the full space of focal mechanism.

[0057] Step 2, secondary screening of the result with the highest correlation degree of amplitude ratio curve from the initial screening optional solution

[0058] According to the fitting of the initial screening optional solution based on polarity data, the 3 theoretical amplitude ratios of each observation station (for example, N stations, STA1..., STA N A P / A SV , A P / A SH , A SV / A SH Set the secondary screening standard to select the optional solution that meets the amplitude ratio fitting standard. In the past algorithm, the difference between the N station theoretical amplitude ratio and the actual observation amplitude ratio is calculated, and then different statistical methods are used to calculate the statistical minimum value of the N residual errors for secondary screening; while the HCN algorithm is essentially different, which is to calculate the correlation degree of the 3 theoretical amplitude ratios of each observation station forming 3 overall curves (which can also be regarded as 3 arrays) and the 3 curves (or arrays) formed by the actual observation, and select the result with the highest comprehensive correlation degree as the secondary screening result. When calculating the comprehensive correlation degree, the HCN algorithm defaults that the weights of the 3 curves are equal, and also allows users to set different weights according to the actual amplitude reliability evaluation. For the N D For each error disturbance data sample, the HCN algorithm selects the N D The union of the optional solutions independently calculated and secondary screened is taken as the result of this stage.

[0059] Figure 2 ② The blue thin line in the middle represents the secondary screening optional solution, which is the co-averaging of the original data N D ​=10 times of ray error perturbation, the union of the results of the second screening respectively. It can be seen that even if the initial screening result is poor, the amplitude ratio curve correlation can well constrain the selected solution, and N D The selected solution obtained by the second screening of the secondary ray random error perturbation is very concentrated, reflecting the stability of the second screening based on the amplitude ratio curve correlation.

[0060] Step 3, calculate the average tensor final solution and evaluate its quality

[0061] According to N D The union of the selected solutions obtained by the second screening of the secondary perturbation, the average tensor of these selected solutions is calculated, and the final focal mechanism inversion result is determined; then, the polarity weighted fitting degree, the amplitude ratio curve correlation coefficient, and the angular dispersion of the selected solution space of the result and the actual observation data are calculated, and the quality of the final focal mechanism inversion result is comprehensively evaluated. According to a large number of inversion applications and tests of the HCN algorithm, the focal mechanism inversion result with polarity weighted fitting degree ≥0.8, amplitude ratio curve correlation coefficient ≥0.5, and angular dispersion of selected solution space ≤25° is considered to be a very reliable high-quality result.

[0062] If there are reference results given by other methods, such as the focal mechanism given by waveform fitting method, regional stress field results, etc., the HCN algorithm will also calculate the spatial angle difference between the inversion result and the reference result. For result output, the HCN algorithm will output the parameter file of the focal mechanism selected solution at different stages of screening, and also through the built-in result drawing function module, the input data used in the inversion, the selected solution screening at different stages, and the data fitting of the final inversion result are displayed according to the general paradigm of seismology.

[0063] Figure 2 ③The inversion result of the HCN algorithm for the Kamchatka earthquake focal mechanism and the data fitting of the final result are shown. The left figure is a focal sphere projection diagram, in which the green thick line segment surface represents the final focal mechanism inversion result, and the overall fitting of the final result and the input data is displayed as a numerical value on both sides of the focal sphere diagram, and the fitting of each data is displayed in different forms of symbols; in the left figure, different symbols represent P, SV, and SH polarity, and the colors red and black represent whether the theoretical and actual polarity matches or not; the right figure shows the fitting of the theoretical and actual observation amplitude ratio of each station.

[0064] For this Kamchatka 8.7 magnitude earthquake, the Global Centroid Moment Tensor Research Center (GCMT) also gives the focal mechanism research result (red thick line segment surface in the left figure), which can be considered as a reliable result; it can be seen that the inversion result given by the HCN algorithm is quite close to the GCMT result. At the same time, Figure 2 ③The right figure shows that the PS amplitude ratio (A P / A SV , A P / A SH There is a large systematic deviation between the theoretical value and the actual observation value of the amplitude ratio (the difference between the theoretical and actual log values is more than 1, that is, the difference between the theoretical and actual PS amplitude ratios is more than 10 times), which is likely to be caused by the long propagation distance and large amplitude attenuation of teleseismic observation. If the traditional amplitude ratio constraint method is used, it must be reasonably corrected by various technical means, otherwise it is difficult to constrain a reliable and close-to-true focal mechanism solution; and the HCN algorithm is based on the correlation fitting of the amplitude ratio change curve between the theoretical and actual observation values, and it shows very strong stability and accuracy without correction in the face of such amplitude attenuation influence.

[0065] Compared with the same kind of method and improved technology based on the inversion of focal mechanism of polarity and amplitude ratio data in the past, the main advantages of the method of the present application are reflected in the following aspects.

[0066] (1) Compared with the traditional fitting residual statistics, the HCN algorithm fitting of the relative high and low of the amplitude ratio between different stations is more theoretically reasonable. The algorithms for constraining the focal mechanism by using the amplitude ratio in the past are all based on the residual statistics of the actual observation amplitude ratio and the theoretical value of different stations to select the solution, for example, FOCMEC adopts the limited fitting residual threshold and the number of records exceeding the threshold, and HASH and SKHASH (Skoumal et al., 2024) adopt the minimum sum of absolute residual values. Since the theoretical basis for constraining the focal mechanism by using the amplitude ratio is that the amplitude ratio of different stations is controlled by the position of the focal sphere, there must be a correlation between the theoretical and actual observation amplitude ratio, and this correlation is also supported by a large amount of amplitude ratio measurement data. Therefore, in principle, the correlation fitting of the theoretical and actual amplitude ratio used by the HCN algorithm can more accurately match the relative high and low difference of the amplitude ratio between different stations. In contrast, even the selected solution that meets the residual statistics fitting standard in the past algorithm may still have a large uncertainty in the relative high and low of the amplitude ratio between different stations. For example, assuming that the residual of the theoretical and actual amplitude ratio log values of two stations is 0.4, the relative high and low of the amplitude ratio of the two stations may differ by 6.3 times at most, 0.8 which is enough to cause a wrong judgment of the "true" position of the two stations on the focal sphere, so that it is difficult to constrain a more reasonable result close to the truth.

[0067] (2) For the amplitude ratio data which is greatly affected by regional medium differences and path attenuation and other factors, the most prominent theoretical advantage of the HCN algorithm is that the correlation fitting has linear transformation invariance, that is, This linear transformation invariance shows that the HCN algorithm is fundamentally consistent with the technical improvement proposed by Trugman to correct the amplitude ratio using linear transformation, while fundamentally avoiding the uncertainty of transforming parameters when applied in different regions. Returning to the basic principles, formulas (2) and (3) show that when using the S / P amplitude ratio to constrain the focal mechanism in different regions, even if the overall wave velocity ratio in different regions changes significantly (for example, from 1.4 to 2.0), the regional difference in wave velocity ratio only brings a constant amplitude scaling (or constant amplitude shift in log value) of the amplitude ratio, which has no effect on the HCN algorithm that uses the correlation fitting of the theoretical and actual amplitude ratio. The same is true for the path attenuation factor. The systematic overall attenuation of seismic wave amplitude does not affect the correlation fitting of the HCN algorithm. Therefore, it can be said that due to the innovative conversion to amplitude ratio correlation fitting, the HCN algorithm can fundamentally eliminate regional wave velocity ratio differences and systematic deviations of path attenuation, leaving only local differences between different stations that are often high-order small quantities relative to systematic deviations.

[0068] Taking the case of Kamchatka seismic data as an example, assume that the PS theoretical amplitude ratio of 5 observation stations corresponding to the true focal mechanism is [22, 42, 62, 82, 102]; if the path attenuation of P and S waves from the source to the observation station is significantly different, the PS wave velocity ratio is significantly different from the theoretical value, and other factors, the actual observed amplitude ratio may be approximately [1, 2, 3, 4, 5]. If the previous amplitude ratio fitting method is used, the residual error of the theoretical and actual amplitude ratio of each station is large, and the overall similarity coefficient is also low, which is obviously difficult to screen out reasonable results without reasonable correction. For such systematic deviations that can change the results, existing various improved algorithms and techniques, including REFOC, RFMI, and other relative focal mechanism inversion algorithms, and amplitude ratio correction techniques, have many limitations such as requiring more prior information, requiring high computational cost, and having lower universality, while the HCN algorithm is based on the correlation fitting of the theoretical and actual observed amplitude ratio curve, which fundamentally avoids the interference of such systematic deviations. This principle advantage is evident in the Kamchatka earthquake inversion.

[0069] (3) The HCN algorithm reasonably considers various data errors, which can reduce the probability of losing the "true" solution and improve the reliability of the final focal mechanism inversion result. For polarity picking errors, the HCN algorithm sets a polarity fitting under-float amplitude threshold ( DPF) is 10%, which is higher than the corresponding range of many algorithms of the same kind, and will naturally increase the dispersion of the preliminary screening results. The reason for this setting is that previous studies have shown that the initial polarity picking error rate is 10%~20%, so the results within this range of decline may be real results. If the same kind of algorithm is set too low, it may cause the "false" concentration of optional solutions and the loss of "true" solutions. On the other hand, the confidence in the reliability of the HCN algorithm is that even if the preliminary screening solution is very dispersed, the HCN algorithm can still reasonably and efficiently constrain the reliable focal mechanism results through the amplitude ratio correlation fitting of the secondary screening. At the same time, the HCN algorithm dynamically adjusts the polarity data weight based on the P, SH, and SV radiation patterns on the basis of user-defined initial weights, which can avoid excessive control of polarity data on the focal mechanism plane and interfere with the reasonable evaluation of the polarity fitting degree of the final result. In addition, according to the research accumulation, the off-source angle and azimuth error disturbance are reasonably set according to the ray path, and the union of the optional solutions selected by multiple disturbance tests further reduces the probability of losing "true" solutions and improves the reliability of the final results and the rationality of their evaluation.

[0070] In addition, compared with the step-by-step screening method used by the HCN algorithm, some algorithms of the same kind may consider a comprehensive screening method, that is, to calculate the normalized sum of multiple data (such as polarity and amplitude ratio) fitting degrees to get a comprehensive fitting degree, and then select the result with the highest comprehensive fitting degree. However, different data fitting degrees have great differences in sensitivity to the difference between optional solutions and real solutions, so the comprehensive fitting degree obtained by normalizing the sum of different data fitting degrees cannot achieve the effect of joint constraint of multiple data, and the stability of the results screened by it is also difficult to guarantee. It is precisely because of the different sensitivities of polarity and amplitude ratio data to the difference between optional solutions and real solutions that the HCN algorithm uses a step-by-step screening method, and the stability and rationality of this method have been verified in the development and testing process of our algorithm.

[0071] (4) HCN algorithm does not require prior information such as reliable nearby focal mechanism, regional accurate wave velocity structure, etc. It does not depend on the conditions such as clustered distribution of earthquakes, dense observation, etc. It also does not involve large-scale operation of earthquake clusters. The parameter setting has scientific basis and does not need to be changed in most cases. These all show that it has many technological advantages such as low requirement for prior information, low requirement for computing power, strong universality, etc. At the same time, HCN algorithm is developed based on MATLAB language, which integrates all chain function modules from data reading, inversion calculation, disturbance test to result drawing. It can realize the process from input data to result display with one key, has excellent user-friendly degree and convenient use. It uses matrix operation instead of a large number of loops to significantly improve the efficiency of focal mechanism inversion. It also supports parallel computing, which can realize task distribution through the multi-process module of MATLAB to quickly improve the computing task of a large number of focal mechanism inversions. These technological and practical advantages make HCN algorithm can be used as a classic improved algorithm of focal mechanism, support the application and promotion of small earthquake focal mechanism research, and play a greater value in other major earth science problems.

[0072] (5) HCN algorithm is applied to different seismic research scenarios such as teleseisms, regional earthquakes, and small and micro mine earthquakes. Two different types of seismic activities, strike-slip and dip-slip, are selected for each scenario. HCN algorithm shows excellent accuracy and robustness in the inversion of actual seismic data in these different scenarios and different mechanisms. For example, Figure 3 and Figure 4The actual results of the inversion of the two teleseismic event data (Myanmar earthquake and Kamchatka earthquake), the two regional seismic data (Gaoyou earthquake in Jiangsu and Changning aftershock in Sichuan), and the two small mine earthquake data (Yanzhou coal mine earthquake in Shandong and dry hot rock mine earthquake in Gonghe in Qinghai) are shown from top to bottom. The left figure is the distribution of the polarity data and the focal mechanism inversion results on the focal sphere: + and - represent the picking polarity of P wave in the vertical direction, positive (upward) and negative (downward), respectively; F and B represent the picking polarity of SV wave in the radial direction, positive (away from the focus) and negative (towards the focus), respectively; R and L represent the picking polarity of SH wave in the tangential direction, positive (right, as viewed from the earthquake to the station) and negative (left, as viewed from the earthquake to the station), respectively; the colors of these symbols, red and black, represent matching or not matching the final results, respectively; the yellow thin line represents the optional solution obtained based on the polarity weighted fitting primary screening, the blue thin line represents the optional solution selected in the secondary screening, the green thick line represents the final determined focal mechanism inversion result, and the red thick line represents the reference focal mechanism solution. The right figure is the variation curve of the amplitude ratio between the theoretical and actual observation at different stations. The above results are 10 sample results considering the random disturbance of the ray path error. It can be seen that even if the polarity primary screening result is low in constraint due to limited data, the HCN algorithm can still accurately constrain the approximate true focal mechanism result (for example, taking the results of the Global Centroid Moment Tensor Research Center or the waveform fitting algorithm as the reference) through the correlation degree fitting of the amplitude ratio curve, and it performs very stably in the error disturbance test of the ray path. In addition, compared with the limitation of the waveform fitting algorithm that is difficult to apply to the focal mechanism inversion of small earthquakes, the HCN algorithm obviously has stronger universality; at the same time, the HCN algorithm has low demand for prior information of the research area and the earthquake itself and small demand for computer algorithm power, which provides a general, simple and efficient technical approach for accurate determination of the focal mechanism of earthquakes in various scenarios.

[0073] Those skilled in the art can understand that all or part of the processes in the above-mentioned embodiment methods can be completed by a computer program instructing related hardware. The program can be stored in a computer readable storage medium, and when executed, can include the processes of the above-mentioned embodiments of each method. The storage medium can be a magnetic disc, an optical disc, a read-only memory (ROM) or a random access memory (RAM), etc.

[0074] It should be understood that the above detailed description of the technical solutions of the present application by means of preferred embodiments is illustrative rather than limiting. Those skilled in the art can modify the technical solutions recorded in each embodiment on the basis of the description of the present application, or make equivalent substitutions for part of the technical features; and these modifications or substitutions do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of each embodiment of the present application.

Claims

1. A focal mechanism inversion method based on polarity and amplitude ratio correlation fitting, characterized in that: The method comprises: According to the quantitative seismology theory, the regional medium difference and the systematic deviation of the path attenuation to the actual amplitude ratio data are quantitatively analyzed; The HCN algorithm is used to jointly invert the focal mechanism based on the polarity and amplitude ratio data of P waves, SV waves and SH waves. First, a plurality of data sets are generated by randomly disturbing the original data, and the initial screening selectable solutions are fitted based on the polarity data of P waves, SV waves and SH waves in the whole space of the focal mechanism parameters. Then, the fitting results with the best correlation degree between the theoretical and actual observation amplitude ratio curves are selected as the selectable solutions for secondary screening. The union of the secondary screening results of a plurality of random error disturbances is taken, and the final focal mechanism inversion result is determined by calculating the average tensor. According to the polarity weighted fitting degree, the amplitude ratio curve correlation coefficient and the angular dispersion of the selectable solution space of the result and the actual observation data, the quality of the final focal mechanism inversion result is comprehensively evaluated.

2. The method according to claim 1, wherein the method is characterized by: The quantitative analysis of the systematic deviation comprises: According to the quantitative seismology theory, the principle of amplitude ratio inversion focal mechanism is that the radiation patterns of P waves and S waves in the wave displacement field generated by a double force couple point source are different, which can be expressed as: （1） Considering that the seismic wave propagates from the source to the station with a path-dependent attenuation, which can be expressed as where, 、 、 respectively represent the frequency, travel time, and attenuation factor of the seismic wave; in addition, considering that the actual PS wave velocity ratio can have regional differences, the system bias caused by the two can be quantitatively expressed as: （2） If the logarithmic value (log value) commonly used in amplitude ratio fitting is used, formula (2) can be written as: （3） In the above formula, A P 、A SV 、A SH respectively represent the amplitudes of P, SV and SH waves, α、β respectively represent the P and S wave propagation velocities, θ is the angle between the ray and the normal of the fault plane, is the angle between the projection of the ray in the fault plane and the fault slip direction, and the superscript theo represents the theoretical value, and the superscript prac represents the actual observed value; the formulas (2) and (3) quantitatively reflect that the regional medium difference and the path attenuation factor will cause systematic deviation between the theoretical and actual amplitude ratio data.

3. The method of claim 1, wherein the method is characterized by: The inversion process of the HCN algorithm comprises: Step 1, using the polarity data to weight and fit the initial screening selectable solutions; Step 2, secondary screening the result with the highest amplitude ratio curve correlation degree from the initial screening selectable solutions; Step 3, calculating the average tensor final solution and evaluating its quality.

4. The method according to claim 3, wherein the method is characterized by: The use of polarity data to weight and fit the initial screening selectable solutions comprises: The perturbation errors of the take-off angle and the azimuth angle are respectively err ih = 5° and err az = 2°, the take-off angle and the azimuth angle of the original observation data are respectively within ±err ih and ±err az , and the N D first (default 10) random perturbation inversion (default the first one is not perturbed) generates N D perturbed data samples; according to the spatial angle step of 5°, the whole space (strike of fault plane 0°~355°, dip of fault plane 0°~90°, slip angle of fault plane 0°~355°) of the source mechanism parameters is divided into 72×19×72=98496 kinds of parameter combinations; for each source mechanism parameter combination, the theoretical polarity of each station is calculated and compared with the actual observation polarity data, and the parameter combination with the weighted polarity fitting degree (R) PF not less than the maximum polarity fitting degree (R PF M ) minus the floating amplitude threshold (ΔR DPF ) of all parameter combinations is selected as the primary screening optional solution, that is, PF ≥ PF M - RPF; Since there are N D error perturbation data samples, the primary screening optional solution is the union of N D primary screening solutions.

5. The method of claim 3, wherein the method is characterized by: The secondary screening of the result with the highest amplitude ratio curve correlation degree from the initial screening selectable solutions comprises: According to the polarity data fitting preliminary screening can be selected solution in each observation station (N stations, STA1..., STA N ) of the 3 theoretical amplitude ratio (A1, A2, A3) A P / A SV , A P / A SH , A SV / A SH ) of the 3 overall curve, and the actual observation of the 3 curve correlation calculation, select the highest overall correlation results as the secondary screening results; in the calculation of the overall correlation, the weight of the 3 curve is equal, or according to the actual amplitude reliability evaluation set different weight; for N D error disturbance data samples, select the union of the N D independent calculation and secondary screening of the optional solution as the result of this stage.

6. The method of claim 3, wherein the method is characterized by: The calculation of the average tensor final solution and the evaluation of its quality comprises: To N D The optional solutions selected by the two secondary random error disturbances are taken together, the average tensor of these optional solutions is calculated, and thus the final focal mechanism inversion result is determined. Then, the polarity weighted fitting degree, the amplitude ratio curve correlation coefficient, and the angular dispersion of the optional solution space of the inversion result and the actual observation data are calculated, and thus the quality of the final focal mechanism inversion result is comprehensively evaluated.

Images