Focus mechanism inversion method based on correlation fitting of polarity and amplitude ratio

Abstract

The invention provides a seismic source mechanism inversion method based on polarity and amplitude ratio correlation fitting, and belongs to the seismic correlation field, the method firstly quantitatively analyzes systematic deviation brought by regional medium difference and path attenuation to actual amplitude ratio data, and then adopts an HCN algorithm to jointly invert a seismic source mechanism based on polarity and amplitude ratio data. The algorithm performs random error disturbance on original data to generate a plurality of groups of data sets, and preliminarily screens selectable solutions in a full space of seismic source mechanism parameters based on polarity data fitting of P, SV and SH waves; screening the result with the highest correlation between the theoretical amplitude ratio curve and the actual amplitude ratio curve for the second time; and finally, calculating an average tensor of secondary screening results of the multiple groups of error disturbance data sets, determining a final seismic source mechanism inversion result, and evaluating a fitting degree and a dispersion degree of the final seismic source mechanism inversion result and actual observation data. The method can effectively cope with the influence of regional medium difference and path attenuation on amplitude ratio data, improves the inversion accuracy of a seismic source mechanism, and provides a reliable basis for seismic research.

Description

Technical Field

[0001] The present invention belongs to the field of earthquake-related fields, and more particularly relates to a focal mechanism inversion method based on polarity and amplitude ratio correlation fitting. Background Art

[0002] Earthquake focal mechanisms are key data for revealing regional stress states, fault activity characteristics, and earthquake genesis. The accuracy of their inversion directly impacts our understanding of these issues. For moderate-to-strong earthquakes, waveform-fitting methods such as gCAP can effectively invert their focal mechanisms. However, for small earthquakes (e.g., M ≤ 3.5), waveform-fitting inversion methods struggle to produce reliable results due to factors such as low waveform signal-to-noise ratios and incompatibility with the effective signal frequency band. However, small earthquakes, due to their high frequency, wide spatial distribution, and sensitivity to local stress field variations, often hold great potential for detailed investigation of these deep-seated Earth science questions (Hu Xingping et al., 2021). Therefore, developing robust and reliable methods for inverting the focal mechanisms of small earthquakes is crucial. Constraining them based on polarity and amplitude ratio data is an important approach for inverting focal mechanisms, especially for small earthquakes. Compared to purely polarity-based inversion methods, such as FPFIT, the inclusion of amplitude ratio data can compensate for the limited quantity and distribution of polarity data, providing more constraints for focal mechanism inversion. This approach remains a crucial approach for many researchers studying the focal mechanisms of small earthquakes. Currently, the classic algorithms for inverting focal mechanisms using polarity and amplitude ratio data include FOCMEC, which uses the polarity of P waves, SV waves, and SH waves and the amplitude ratio between them, and HASH, which uses the polarity of P waves and the amplitude ratio of PS waves, and its Python improved version SKHASH.

[0003] However, the actual observed seismic wave amplitude is largely affected by factors such as regional medium differences and path attenuation, and there are systematic deviations. If it is used directly without correction, its reliability is difficult to guarantee, and it has even caused academic controversy on whether the amplitude ratio is useful in the actual application of focal mechanism inversion.

[0004] To address this issue, subsequent researchers have proposed various improved algorithms and techniques based on different approaches. These include relative focal mechanism inversion algorithms such as the joint inversion of relative radiation patterns between adjacent earthquakes (REFOC), inversion of the relative focal mechanism of adjacent minor earthquakes (RFMI) using the reliable focal mechanism of the main earthquake, and amplitude ratio correction techniques such as correcting amplitude ratios using the average theory of dense DAS observations and observed amplitude ratios, and correcting amplitude ratios using linear transformations derived from analysis of large amounts of amplitude ratio data. However, each of these improved algorithms and techniques has limitations. For example, the REFOC algorithm relies on the clustered distribution of earthquakes and is not applicable to scattered earthquakes. The RFMI algorithm truly resolves the deviations of minor earthquakes relative to the main event, and the reliability of its absolute mechanisms for minor earthquakes needs further clarification. Furthermore, relative algorithms such as REFOC and RFMI are computationally expensive and may be limited by computing power for large-scale application. Correcting amplitude ratios using the average of dense observations such as DAS is less applicable when observation stations are dispersed. Correction methods using linear transformations derived from analysis of large amounts of amplitude ratio data require verification of the specific transformation parameters for application in different regions.

[0005] It can be said that previous improved algorithms and technologies have many limitations when dealing with systematic deviations in seismic wave amplitude caused by factors such as regional medium differences and path attenuation. These limitations include the need for a lot of prior information, high computing power requirements, and low universality. These limitations make it difficult to effectively support the research and development of small earthquake source mechanisms and other related earth science issues. Summary of the Invention

[0006] To address this issue, the present invention, based on the fundamental theory of quantitative seismology, quantitatively analyzes the systematic biases in amplitude ratio data caused by factors such as regional dielectric differences and path attenuation. Based on the MATLAB language, the HCN algorithm is developed to invert the focal mechanism by fitting the polarity and amplitude ratio data of P-waves, SV-waves, and SH-waves. This method employs a novel amplitude ratio data fitting approach during the focal mechanism inversion process, namely, fitting the correlation between theoretical and observed amplitude ratio curves between different stations. This method can, in principle, eliminate the systematic biases caused by regional dielectric differences and propagation path attenuation, thereby maximizing the reliability and accuracy of the focal mechanism inversion results.

[0007] In order to achieve the above object, the present invention is implemented by adopting the following technical solutions: the method comprises: Based on quantitative seismological theory, the systematic deviations of the actual amplitude ratio data caused by regional medium differences and path attenuation 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, the original data are perturbed with random errors to generate multiple data sets. The optional solutions are preliminarily screened based on the polarity data of P waves, SV waves, and SH waves in the entire focal mechanism parameter space. Then, the best fitting result of the correlation between the theoretical and actual observed amplitude ratio change curves is selected as the optional solution for secondary screening. The secondary screening results of multiple sets of random error perturbations are unioned, and the final focal mechanism inversion result is determined by calculating the average tensor. The quality of the final focal mechanism inversion result is comprehensively evaluated based on the polarity-weighted fitting degree of the result with the actual observed data, the correlation coefficient of the amplitude ratio curve, and the spatial angular dispersion of the optional solution.

[0008] In one embodiment, the quantitative analysis of system deviation is specifically as follows: According to quantitative seismology theory, the principle of using amplitude ratio to invert focal mechanism is based on the fact that P-waves and S-waves radiate differently in the wave displacement field generated by a double-couple point source, which can be expressed as: (1) Considering that the seismic wave will attenuate from the source to the station in accordance with the propagation path, it can be expressed as ,in, 、 、 represent the frequency, travel time, and attenuation factor of the seismic wave, respectively. In addition, considering the possible regional differences in the actual PS wave velocity ratio, the systematic deviation caused by these two factors can be quantitatively expressed as: (2) If the logarithmic value (log value) commonly used in amplitude ratio fitting is used, equation (2) can be written as: (3) In the above formula, A P 、A SV 、A SH Represent the amplitudes of P wave, SV wave and SH respectively. α、β denote the propagation speeds of P and S waves, respectively. θ is the angle between the ray and the normal to the fault plane, is the angle between the projection of the ray in the fault plane and the fault slip direction, theo Indicates theoretical value, upper left prac represents the actual observed value; 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.

[0009] In one embodiment, the HCN algorithm inversion process includes: Step 1: Use polarity data weighted fitting to preliminarily screen for optional solutions; Step 2: Secondarily screen the results with the highest correlation of amplitude ratio curve from the initially screened optional solutions; Step 3: Calculate the average tensor final solution and evaluate its quality.

[0010] In one embodiment, the method of using polarity data weighted fitting to initially screen for optional solutions includes: Set the disturbance errors of the source angle and azimuth angle to err ih =5° and err az =2°, the angle from the source and the azimuth of the original observation data are ±err ih ±err az Do within the scope N D times (default is 10 times) random perturbation inversion (default is no perturbation for the first time), generating N D disturbance data samples; with a spatial angle step of 5°, the entire focal mechanism parameter space (fault plane strike 0°~355°, fault plane dip 0°~90°, fault slip angle 0°~355°) is divided into 72×19×72=98496 parameter combinations; for each focal mechanism parameter combination, the theoretical polarity of each station is calculated and compared with the actual observed polarity data, and the weighted polarity fitting degree ( PF ) is not less than the maximum polarity fit among all parameter combinations ( PF M ) minus the downward amplitude threshold ( DPF ) parameter combination as the initial screening optional solution, that is PF ≥ PF M - RPF; Due to the existence N D error-perturbed data samples, the initial screening solution is N D The union of the first sieve solutions.

[0011] In one embodiment, the secondary screening of the results with the highest correlation of the amplitude ratio curves from the primary screening optional solutions specifically includes: Based on the polar data fitting, the optional solutions are preliminarily screened at each observation station (N stations, STA1...,STA N )'s three theoretical amplitude ratios ( A P / A SV , A P / ASH , A SV / A SH ) are correlated with the three curves composed of actual observations, and the result with the highest comprehensive correlation is selected as the secondary screening result; when calculating the comprehensive correlation, the weights of the three curves are equal, or different weights are set according to the actual amplitude reliability assessment; for N D error disturbance data samples, select this N D The result of this stage is the union of the optional solutions obtained through the second independent calculation and the second screening.

[0012] In one embodiment, computing the average tensor final solution and evaluating its quality includes: right N D The union of the optional solutions screened out by the random error perturbation is taken, and the average tensor of these optional solutions is calculated to determine the final focal mechanism inversion result; then, the polarity-weighted fitting degree of the inversion result with the actual observation data, the correlation coefficient of the amplitude ratio curve, and the spatial angular discreteness of the optional solution are calculated to comprehensively evaluate the quality of the final inversion result of the focal mechanism.

[0013] Beneficial effects of the present invention: 1. Previous algorithms that use amplitude ratios to invert focal mechanisms all use amplitude ratio residual statistics to screen optional solutions. In comparison, the HCN algorithm is more theoretically reasonable in fitting the relative high and low amplitude ratios between different stations.

[0014] 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 differences and path attenuation in principle.

[0015] 3. The HCN algorithm accounts for errors such as polarity picking and ray path perturbations, reducing the probability of missing the "true" solution and improving the reliability of the final focal mechanism inversion results. Furthermore, a step-by-step screening approach is used for polarity and amplitude ratio data, rather than a comprehensive approach, which is more rational and stable.

[0016] 4. The HCN algorithm has many technological and practical advantages, such as less prior information, low computing power cost, strong universality, and good convenience. It also shows excellent inversion effects in different actual scenarios, providing a universal, simple and efficient technical approach for the accurate determination of focal mechanism solutions in multiple scenarios. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 This is a technical framework diagram of the method of the present invention; Figure 2 This is an example of the actual data inversion steps of the method of the present invention; Figure 3 Inversion results for different scenarios such as teleseismic, regional earthquake, and small mine earthquake (a); Figure 4 (b) Inversion results for different scenarios such as teleseismic, regional earthquake, and small mine earthquake. DETAILED DESCRIPTION

[0018] To facilitate understanding of the present invention, the present invention will be described more fully below with reference to the accompanying drawings. The drawings illustrate exemplary embodiments of the present invention. However, the present invention may be implemented in many different forms and is not limited to the embodiments described herein. Rather, these embodiments are provided to provide a more thorough and comprehensive understanding of the present invention.

[0019] Unless otherwise defined, all technical and scientific terms used herein have the same meanings as those understood by those skilled in the art to which the present invention pertains. The terms used in the present specification are for the purpose of describing specific embodiments only and are not intended to limit the present invention. To facilitate understanding of the present invention, a more comprehensive description of the present invention will be provided below with reference to the accompanying drawings. Typical embodiments of the present invention are shown in the drawings. However, the present invention may be embodied in many different forms and is not limited to the embodiments described herein. Rather, these embodiments are provided to provide a more thorough and comprehensive understanding of the present invention.

[0020] The principle of the present invention is as follows: the HCN algorithm first performs random error perturbation on the original data to generate multiple sets of data sets, and then preliminarily screens the optional solutions based on the polarity data fitting of P waves, SV waves, and SH waves in the entire space of focal mechanism parameters, and then selects the best fitting result of the correlation between the theoretical and actual observed amplitude ratio change curves as the optional solution for secondary screening; the secondary screening results of multiple sets of random error perturbations are combined, and the final focal mechanism inversion result is determined by calculating the average tensor, and the quality of the final focal mechanism inversion result is comprehensively evaluated based on the polarity-weighted fitting degree of the result with the actual observed data, the correlation coefficient of the amplitude ratio curve, and the spatial angular dispersion of the optional solution.

[0021] In addition to being able to effectively cope with the influence of regional medium differences and propagation path attenuation, thereby improving the accuracy of focal mechanism inversion, the HCN algorithm also has many scientific and technological advantages and practical advantages, such as less prior information, low computing power cost, strong universality, and good convenience. It can provide efficient and reliable scientific and technological support for the research of earthquake focal mechanism and other related earth science issues.

[0022] like Figure 1 As shown in Figure 1, a focal mechanism inversion method based on polarity and amplitude ratio correlation fitting includes: The principle of using amplitude ratio to invert focal mechanism is based on the fact that P-waves and S-waves radiate differently in the wave displacement field generated by a double-couple point source. By fitting the amplitude ratios of the waveforms recorded at different stations, the spatial angular relationship between the exit points of the seismic rays from these stations back to the source sphere and the normal and slip directions of the fault plane can be inferred, thereby determining the focal mechanism. According to quantitative seismology theory, this principle can be expressed as: (1) In formula (1), A P 、A SV 、A SH Represent the amplitudes of P wave, SV wave and SH respectively. α、β denote the propagation speeds of P and S waves, respectively. θ is the angle between the ray and the normal to the fault plane, is the angle between the projection of the ray in the fault plane and the fault slip direction, theo Indicates theoretical value.

[0023] Based on the basic principle of using polarity and amplitude ratio to invert the focal mechanism expressed in formula (1), considering that the seismic wave will attenuate due to the propagation path when it propagates from the source to the station, it can be expressed as ,in, 、 、 represent the frequency, travel time, and attenuation factor of the seismic wave, respectively. In addition, considering the regional differences in the PS wave velocity ratio in practice, the systematic deviation caused by these two factors can be quantitatively expressed as: (2) If the logarithmic value (log value) commonly used in amplitude ratio fitting is used, equation (2) can be written as: (3) The upper left subscript in formula (2) and (3) prac These represent actual observations, quantitatively reflecting the systematic deviations between theoretical and actual amplitude ratio data caused by factors such as regional media differences and path attenuation. Effectively eliminating these errors is key to obtaining accurate and reliable focal mechanisms. Previous inversion methods and technical improvements have each had their limitations.

[0024] The HCN algorithm of the present invention uses polarity and amplitude ratio data of P waves, SV waves, and SH waves to jointly invert focal mechanisms, maximizing the utilization of polarity and amplitude ratio data. The algorithm's steps include: ① using polarity data weighted fitting to initially screen for possible solutions; ② secondary screening of the initially screened possible solutions for the results with the highest correlation in amplitude ratio curves; and ③ calculating the average tensor final solution and evaluating its quality. Step ② represents an innovative solution proposed by the HCN algorithm to address the aforementioned systematic deviations. It employs a novel amplitude ratio data fitting method—converting the residual fit between the actual observed amplitude ratios and theoretical values, commonly used in previous algorithms, into a correlation fit between the observed amplitude ratio variation curves and the theoretical variation curves across different stations. The criteria for selecting focal mechanisms based on amplitude ratios are converted from the statistical minimum of the residual between the actual observed amplitude ratios and the theoretical values ​​to the highest correlation between the observed amplitude ratio variation curves and the theoretical curves across different stations, minimizing the impact of abnormal observations.

[0025] Figure 2 Using the actual polarity and amplitude ratio data from the magnitude 8.7 Kamchatka earthquake on July 30, 2025 (Beijing time), as an example, we will describe the HCN algorithm's focal mechanism inversion process. Obtaining polarity and amplitude ratio data from raw seismic data by manually picking them is a data preprocessing step, or it can be accomplished using seismic signal picking technology. The HCN algorithm's focal mechanism inversion process begins with the polarity and amplitude ratio data already acquired.

[0026] Step 1: Use polarity data weighted fitting to initially screen for optional solutions The HCN algorithm takes into account the fact that the errors of factors such as earthquake location and regional velocity structure will bring uncertainty to the position of the observation station back to the source sphere (reflected in the off-source angle and azimuth of the input data). The perturbation errors of the off-source angle and azimuth are set to err ih =5° and err az =2°, the angle from the source and the azimuth of the original observation data are ±err ih ±err az Do within the scope N D times (default is 10 times) random perturbation inversion (default is no perturbation for the first time), generating N D This error range is set based on the previous research and understanding of predecessors and the applicant.

[0027] According to the spatial angle step of 5°, the entire space of focal mechanism parameters (fault plane strike 0°~355°, fault plane dip 0°~90°, fault slip angle 0°~355°) is divided into 72×19×72=98496 parameter combinations; this 5° step is the value that balances the 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 weighted polarity fit ( PF ) is not less than the maximum polarity fit among all parameter combinations ( PF M ) minus the downward amplitude threshold ( DPF ) parameter combination as the initial screening optional solution, that is PF ≥ PF M - RPF Due to the existence N D error-perturbed data samples, the initial screening solution is N D The union of the first sieve solutions.

[0028] The HCN algorithm will DPF The value is set to 10% because statistics show that the error rate in picking initial polarity reaches 10% to 20%. The polarity data weighting in the weighted polarity fit involves two levels of weighting. First, the user can pre-set the initial data weight based on the quality of polarity picking, assigning a value between -1 and 1, where ± indicates positive or negative polarity, and the absolute value is the set initial weight. Second, the HCN algorithm dynamically adjusts the weights based on the initial weights according to the radiation patterns of the P, SH, and SV waves. For example, near the nodal plane, where the theoretical amplitude of the P wave is very small, the P-wave polarity weight is reduced. This weighting method is used because polarity with small theoretical amplitude is more prone to misinterpretation in the presence of noise and because the position of station observation data backtracked on the focal sphere has certain errors. Using dynamic weighting can prevent polarity data from over-controlling the focal mechanism nodal plane and interfering with the proper evaluation of the polarity fit of the final result.

[0029] Figure 2 Figure ① shows the distribution of polarity data and preliminary screening results for focal mechanisms on the focal sphere: + and - indicate positive (upward) and negative (downward) vertical polarity for P waves, respectively; F and B indicate positive (away from the earthquake source) and negative (toward the earthquake source) radial polarity for SV waves, respectively; R and L indicate positive (rightward) and negative (leftward) tangential polarity for SH waves, respectively; the thin yellow line indicates the possible solutions obtained through preliminary screening using polarity-weighted fitting. It can be seen that due to the suboptimal amount of polarity data and its distribution on the focal sphere, the constraints on possible solutions are relatively poor, resulting in a highly discrete preliminary screening result across the entire focal mechanism space.

[0030] Step 2: Secondarily select the result with the highest correlation of amplitude ratio curve from the initial screening optional solutions Based on the polar data fitting, the optional solutions are preliminarily screened at each observation station (for example, N stations, STA1...,STA N )'s three theoretical amplitude ratios ( A P / A SV , A P / A SH , A SV / A SH ) sets the secondary screening criteria and selects the optional solutions that meet the amplitude ratio fitting criteria. In the past, the algorithms calculated the difference between the theoretical amplitude ratio and the actual observed amplitude ratio of the N stations during this process, and then used different statistical methods to calculate the statistical minimum of the N residuals for secondary screening; the HCN algorithm is essentially different. It combines the three theoretical amplitude ratios of each observation station into three overall curves (which can also be regarded as three arrays), and calculates the correlation with the three curves (or arrays) composed of actual observations, and selects the result with the highest comprehensive correlation as the secondary screening result. When calculating the comprehensive correlation, the HCN algorithm defaults to equal weights for the three curves, and also allows users to set different weights based on the actual amplitude reliability assessment. For N D error disturbance data samples, the HCN algorithm selects these N D The result of this stage is the union of the optional solutions obtained through the second independent calculation and the second screening.

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

[0032] Step 3: Calculate the average tensor final solution and evaluate its quality according to N DThe final focal mechanism inversion result is determined by calculating the average tensor of the union of the alternative solutions selected by the secondary perturbation. The polarity-weighted fit between this result and the actual observation data, the correlation coefficient of the amplitude ratio curve, and the spatial angular dispersion of the alternative solutions are then calculated to comprehensively assess the quality of the final focal mechanism inversion result. Based on extensive inversion applications and testing of the HCN algorithm, focal mechanism inversion results with a polarity-weighted fit ≥ 0.8, an amplitude ratio curve correlation coefficient ≥ 0.5, and an angular dispersion of the alternative solutions ≤ 25° are considered to be highly reliable and high-quality results.

[0033] If reference results are available from other methods, such as focal mechanisms and regional stress field results from waveform fitting, the HCN algorithm will also calculate the spatial angle difference between the inversion results and the reference results. The HCN algorithm also outputs parameter files for the selected focal mechanism solutions at different stages. The algorithm also uses a built-in result plotting module to display the input data used in the inversion, the selection of alternative solutions at different stages, and the data fit of the final inversion results, using a common seismological paradigm.

[0034] Figure 2 ③ This figure shows the HCN algorithm's inversion results for the Kamchatka earthquake focal mechanism, along with the final data fit. The left figure is a focal sphere projection, where the thick green nodal planes represent the final focal mechanism inversion results. The overall fit between the final results and the input data is numerically displayed on both sides of the focal sphere diagram, with the fit of each data point displayed using different symbols. In the left figure, the distribution of different symbols represents P, SV, and SH polarity, with red and black colors indicating whether the theoretical and actual polarities match, respectively. The right figure shows the fit of the theoretical and actual observed amplitude ratios at each station.

[0035] For the 8.7 magnitude Kamchatka earthquake, the Global Center of Moment Tensor (GCMT) also provided the results of the focal mechanism study (red thick line section in the left figure), which can be considered to be relatively reliable. It can be seen that the inversion results given by the HCN algorithm are quite close to the GCMT results. Figure 2 ③The right figure shows that the PS amplitude ratio ( A P / A SV , A P / A SH) exhibit significant systematic deviations from observed values ​​(most theoretical and actual logarithmic values ​​differ by more than 1, meaning the theoretical and actual PS amplitude ratios differ by more than a factor of 10). This is likely due to the long propagation distances and significant attenuation of teleseismic waveforms. Traditional amplitude ratio constraints require reasonable corrections through various technical means; otherwise, it is difficult to constrain a reliable focal mechanism solution that approximates the "true" value. The HCN algorithm, however, relies on correlation fitting between the theoretical and observed amplitude ratio curves. It exhibits exceptionally robust stability and accuracy despite this strong amplitude attenuation, even without correction.

[0036] Compared with similar methods and improved technologies for inverting focal mechanisms based on polarity and amplitude ratio data, the main advantages of the method of the present invention are reflected in the following aspects.

[0037] (1) Compared with the traditional fitting residual statistics, the HCN algorithm is more theoretically reasonable in fitting the relative amplitude ratios between different stations. In the past, algorithms that used amplitude ratios to constrain focal mechanisms were based on the residual statistics of the actual observed amplitude ratios and theoretical values ​​at different stations to screen optional solutions. For example, FOCMEC used a limited fitting residual threshold and the number of records exceeding the threshold, while HASH and SKHASH (Skoumal et al., 2024) used the minimum sum of the absolute values ​​of the residuals. Given that the theoretical basis of the amplitude ratio constrained focal mechanism is that the amplitude ratios of different stations are controlled by their positions traced back to the source sphere, there must be a correlation between the theoretical and actual observed amplitude ratios, and this correlation is also supported by a large amount of amplitude ratio measured data. Therefore, in principle, the correlation fitting of the theoretical and actual amplitude ratios used by the HCN algorithm can more accurately match the relative differences in amplitude ratios between different stations; in contrast, even for optional solutions that meet the residual statistical fitting criteria in previous algorithms, the relative differences in amplitude ratios between different stations may still have large uncertainties. For example, assuming that the logarithmic residuals of the theoretical and actual amplitude ratios of two stations are both 0.4, the relative difference between the two amplitude ratios is The maximum possible difference is 10 0.8 =6.3 times, which is enough to lead to the misjudgment of the "real" positions of the two on the source sphere, making it difficult to constrain a reasonable result that is closer to the truth.

[0038] (2) The amplitude ratio data is largely affected by factors such as regional medium differences and path attenuation. The most prominent theoretical advantage of the HCN algorithm is that the correlation fitting has linear transformation invariance, that is, This linear transformation invariance indicates that the HCN algorithm fully matches the technical improvement proposed by Trugman to use linear transformation to correct the amplitude ratio in principle, and at the same time fundamentally avoids the uncertainty of parameter transformation that may be required when applied in different regions. Returning to the basic principle, formulas (2) and (3) show that when the S / P amplitude ratio is used 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 about a fixed amplitude scaling of the amplitude ratio (or a fixed amplitude offset of the log value), which has no effect on the HCN algorithm that uses the correlation fitting of theoretical and actual amplitude ratios. The same is true for the path attenuation factor. The systematic overall attenuation of the 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 eliminate the systematic deviation of regional wave velocity ratio differences and path attenuation in principle, and only the local differences between different stations, which are often small and high-order relative to the systematic deviation, may remain.

[0039] Taking the Kamchatka earthquake data as an example, assuming a realistic focal mechanism, the theoretical PS amplitude ratios for the five observation stations are [22, 42, 62, 82, 102]. If factors such as significant differences in the attenuation of the P and S waves from the source to the observation stations and significant differences between the PS velocity ratios and theoretical values ​​exist, the observed amplitude ratios may be approximately [1, 2, 3, 4, 5]. Using conventional amplitude ratio fitting methods, the residuals between the theoretical and actual amplitude ratios for each station are large, and the overall similarity coefficient is also low. Without proper correction, it is difficult to screen out reasonable results. For this kind of systematic deviation that is sufficient to change the results, various existing improved algorithms and technologies, including relative focal mechanism inversion algorithms such as REFOC and RFMI, and amplitude ratio correction technology, have many limitations, such as the requirement of a lot of prior information, high computing power cost, and low universality. The HCN algorithm, based on the correlation fitting of the amplitude ratio change curve between theory and actual observations, fundamentally avoids the interference of this systematic deviation. This principle advantage is clearly demonstrated in the Kamchatka earthquake inversion.

[0040] (3) The HCN algorithm takes into account various data errors in a reasonable manner, which can reduce the probability of losing the "true" solution and improve the reliability of the final focal mechanism inversion results. For polarity picking errors, the HCN algorithm sets the polarity fitting lowering amplitude threshold ( DPF) is 10%, which is higher than the corresponding range of many similar algorithms, naturally increasing the dispersion of the initial screening results. This setting is based on the fact that, on the one hand, previous studies have shown that the initial motion polarity picking error rate reaches 10%-20%. Therefore, results within this drop range are likely to be true. If this setting is too low for similar algorithms, it may lead to "false" concentration of alternative solutions and loss of "true" solutions. On the other hand, it reflects confidence in the reliability of the HCN algorithm. Even if the initial screening solution is very discrete, the HCN algorithm can still reasonably and efficiently constrain reliable focal mechanism results through secondary screening based on amplitude ratio correlation fitting. Furthermore, based on user-defined initial weights, the HCN algorithm dynamically adjusts the polarity data weights according to the respective radiation patterns of P, SH, and SV. This can prevent the polarity data from over-controlling the focal mechanism nodes and interfering with the reasonable evaluation of the polarity fit of the final results. Furthermore, based on the accumulated research, the error perturbations of the ray path from the source angle and azimuth angle are reasonably set, and the union of optional solutions from multiple perturbation tests is selected, which further reduces the probability of losing the "true" solution and improves the reliability of the final result and the rationality of its evaluation.

[0041] Furthermore, compared to the step-by-step screening approach used by the HCN algorithm, some similar algorithms may consider a comprehensive screening approach, where they calculate the normalized sum of the fits of multiple data points (such as polarity and amplitude ratio) to obtain a comprehensive fit, and then select the result with the highest comprehensive fit. However, the sensitivity of different data fits to the difference between the alternative solution and the true solution varies greatly. Therefore, screening based on the comprehensive fit obtained by the normalized sum of the fits of different data points does not achieve the effect of joint constraint on multiple data points, and the stability of the resulting screening results is also difficult to guarantee. It is precisely because of the different sensitivities of polarity and amplitude ratio data to the difference between the alternative solution and the true solution that the HCN algorithm adopts a step-by-step screening approach, and the stability and rationality of this approach have been verified during our algorithm development and testing.

[0042] (4) The HCN algorithm does not require prior information such as reliable focal mechanisms in the vicinity or accurate regional velocity structure. It does not rely on conditions such as earthquake cluster distribution and intensive observation, nor does it involve large-scale calculations of earthquake clusters. The parameter settings have scientific basis and do not need to be changed in most cases. These all indicate that it has many technological advantages such as low requirements for prior information, low computing power cost, and strong universality. At the same time, the HCN algorithm is developed based on the MATLAB language and integrates a full chain of functional modules from data reading, inversion calculation, perturbation testing, to result drawing. It can achieve "one-click" from input data to result display, with excellent user-friendliness and ease of use. It uses matrix operations to replace a large number of loops, significantly improving the efficiency of focal mechanism inversion. It also supports parallel computing and can achieve task distribution through the multi-process module provided by MATLAB, quickly improving the computational tasks of a large number of focal mechanism inversions. These technological and practical advantages make the HCN algorithm a classic improved algorithm for focal mechanism, supporting the application and promotion of small earthquake focal mechanism research, and playing a greater role in other major earth science problems.

[0043] (5) The HCN algorithm is applied to different earthquake research scenarios such as teleseismic, regional earthquake, and small mine earthquake, and two different earthquake activities of strike-slip and dip-slip types are selected for each scenario; the HCN algorithm shows excellent accuracy and robustness in the inversion of actual seismic data in these different scenarios and different mechanisms. Figure 3 and Figure 4As shown, from top to bottom are the actual results of the inversion of two teleseismic event data (the Sagaing earthquake in Myanmar and the Kamchatka earthquake), two regional earthquake data inversion (the Gaoyou earthquake in Jiangsu and the Changning aftershock in Sichuan), and two small mine earthquake data inversion (the Yanzhou coal mine earthquake in Shandong and the Gonghe hot dry rock mine earthquake in Qinghai). The left figure shows the distribution of polarity data and focal mechanism inversion results on the focal sphere: + and - indicate positive (upward) and negative (downward) vertical polarity of P waves, respectively; F and B indicate positive (away from the earthquake source) and negative (toward the earthquake source) radial polarity of SV waves, respectively; R and L indicate positive (rightward) and negative (leftward) tangential polarity of SH waves, respectively. The red and black colors of these symbols indicate a match or mismatch with the final solution, respectively. The thin yellow line represents the alternative solution obtained through the initial screening based on polarity-weighted fitting, the thin blue line represents the alternative solution obtained through the secondary screening, the thick green line represents the final focal mechanism inversion result, and the thick red line represents the reference focal mechanism solution. The right figure shows the variation of the amplitude ratio between theoretical and observed amplitudes at different stations. All of the above results are based on 10 samples, accounting for random perturbations due to ray path errors. It can be seen that even though the initial polarity screening results are very weakly constrained due to limited data, the HCN algorithm can still accurately constrain the focal mechanism results that are close to the true one (for example, using the Global Center of Moment Tensor Research Center or the waveform fitting algorithm results as a reference) through amplitude ratio curve correlation fitting, and it is very robust in ray path error perturbation tests. In addition, compared with the limitation of waveform fitting algorithms that are difficult to apply to the inversion of the focal mechanism of small earthquakes, the HCN algorithm is clearly more universal. At the same time, the HCN algorithm requires less prior information about the study area and the earthquake itself, and requires less computer computing power. This provides a universal, simple, and efficient technical approach for accurately determining the focal mechanism of earthquakes in various scenarios.

[0044] Those skilled in the art will appreciate that all or part of the processes in the above-described method embodiments can be implemented by instructing related hardware through a computer program. The program can be stored in a computer-readable storage medium, and when executed, the program can include the processes in the above-described method embodiments. The storage medium can be a magnetic disk, an optical disk, a read-only memory (ROM), or a random access memory (RAM).

[0045] It should be understood that the detailed description of the technical solutions of the present invention using the preferred embodiments above is illustrative and not restrictive. A person skilled in the art, after reading the present specification, may modify the technical solutions described in the embodiments or replace some of the technical features therein with equivalents; such modifications or replacements do not deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A focal mechanism inversion method based on polarity and amplitude ratio correlation fitting, characterized by: The method includes: Based on quantitative seismological theory, the systematic deviations of the actual amplitude ratio data caused by regional medium differences and path attenuation 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, the original data are perturbed with random errors to generate multiple data sets. The optional solutions are preliminarily screened based on the polarity data of P waves, SV waves, and SH waves in the entire focal mechanism parameter space. Then, the best fitting result of the correlation between the theoretical and actual observed amplitude ratio change curves is selected as the optional solution for secondary screening. The secondary screening results of multiple sets of random error perturbations are unioned, and the final focal mechanism inversion result is determined by calculating the average tensor. The quality of the final focal mechanism inversion result is comprehensively evaluated based on the polarity-weighted fitting degree of the result with the actual observed data, the correlation coefficient of the amplitude ratio curve, and the spatial angular dispersion of the optional solution.

2. The focal mechanism inversion method based on polarity and amplitude ratio correlation fitting according to claim 1, characterized in that: The quantitative analysis of systematic deviation is specifically as follows: According to quantitative seismology theory, the principle of using amplitude ratio to invert focal mechanism is based on the fact that P-waves and S-waves radiate differently in the wave displacement field generated by a double-couple point source, which can be expressed as: （1） Considering that the seismic wave will attenuate from the source to the station in accordance with the propagation path, it can be expressed as ,in, 、 、 represent the frequency, travel time, and attenuation factor of the seismic wave, respectively. In addition, considering the possible regional differences in the actual PS wave velocity ratio, the systematic deviation caused by these two factors can be quantitatively expressed as: （2） If the logarithmic value (log value) commonly used in amplitude ratio fitting is used, equation (2) can be written as: （3） In the above formula, A P 、A SV 、A SH Represent the amplitudes of P wave, SV wave and SH respectively. α、β denote the propagation speeds of P and S waves, respectively. θ is the angle between the ray and the normal to the fault plane, is the angle between the projection of the ray in the fault plane and the fault slip direction, theo Indicates theoretical value, upper left prac represents the actual observed value; 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.

3. The focal mechanism inversion method based on polarity and amplitude ratio correlation fitting according to claim 1, characterized in that: The inversion process of the HCN algorithm includes: Step 1: Use polarity data weighted fitting to preliminarily screen for optional solutions; Step 2: Secondarily screen the results with the highest correlation of amplitude ratio curve from the initially screened optional solutions; Step 3: Calculate the average tensor final solution and evaluate its quality.

4. The focal mechanism inversion method based on polarity and amplitude ratio correlation fitting according to claim 3, characterized in that: The optional solutions for the initial screening using polarity data weighted fitting include: Set the disturbance errors of the source angle and azimuth angle to err ih =5° and err az =2°, the angle from the source and the azimuth of the original observation data are ±err ih ±err az Do within the scope N D times (default is 10 times) random perturbation inversion (default is no perturbation for the first time), generating N D disturbance data samples; with a spatial angle step of 5°, the entire focal mechanism parameter space (fault plane strike 0°~355°, fault plane dip 0°~90°, fault slip angle 0°~355°) is divided into 72×19×72=98496 parameter combinations; for each focal mechanism parameter combination, the theoretical polarity of each station is calculated and compared with the actual observed polarity data, and the weighted polarity fitting degree ( PF ) is not less than the maximum polarity fit among all parameter combinations ( PF M ) minus the downward amplitude threshold ( DPF ) parameter combination as the initial screening optional solution, that is PF ≥ PF M - RPF; Due to the existence N D error-perturbed data samples, the initial screening solution is N D The union of the first sieve solutions.

5. The focal mechanism inversion method based on polarity and amplitude ratio correlation fitting according to claim 3, characterized in that: The secondary screening of the results with the highest correlation of the amplitude ratio curve from the primary screening optional solutions specifically includes: Based on the polar data fitting, the optional solutions are preliminarily screened at each observation station (N stations, STA1...,STA N )'s three theoretical amplitude ratios ( A P / A SV , A P / A SH , A SV / A SH ) are correlated with the three curves composed of actual observations, and the result with the highest comprehensive correlation is selected as the secondary screening result; when calculating the comprehensive correlation, the weights of the three curves are equal, or different weights are set according to the actual amplitude reliability assessment; for N D error disturbance data samples, select this N D The result of this stage is the union of the optional solutions obtained through the second independent calculation and the second screening.

6. The focal mechanism inversion method based on polarity and amplitude ratio correlation fitting according to claim 3, characterized in that: The calculation of the average tensor final solution and evaluation of its quality include: right N D The union of the optional solutions screened out by the random error perturbation is taken, and the average tensor of these optional solutions is calculated to determine the final focal mechanism inversion result; then, the polarity-weighted fitting degree of the inversion result with the actual observation data, the correlation coefficient of the amplitude ratio curve, and the spatial angular discreteness of the optional solution are calculated to comprehensively evaluate the quality of the final inversion result of the focal mechanism.

Images