Three-dimensional Rayleigh wave attenuation relative variation imaging method with stable plot as reference value

Abstract

The invention discloses a three-dimensional Rayleigh wave attenuation relative variation imaging method taking a stable plot as a reference value, relates to the technical field of Rayleigh wave attenuation relative variation measurement, and aims to solve the problems that in the prior art, a seismic source and a path effect are difficult to effectively strip, and a medium attribute change signal cannot be directly extracted. The invention aims to improve the capability of monitoring the change of the state of the earth crust medium and improve the recognition precision of the earthquake precursor abnormity. Comprising the following steps: respectively calculating Rayleigh wave amplitude ratios of seismic stations in a research area and a reference area to two seismic events; correcting the amplitude ratio of the seismic station in the research area by using the amplitude ratio of the seismic station in the reference area, and calculating to obtain the relative variation of the Rayleigh wave attenuation coefficient of the medium of the station in the research area between the two seismic events; on the basis of relative variation of attenuation coefficients calculated by multiple stations and multiple frequencies, spatial interpolation is performed on each single-frequency data to generate a two-dimensional plane distribution diagram, and then different-frequency two-dimensional diagrams are combined according to corresponding depths to generate a three-dimensional space distribution image.

Description

Technical Field

[0001] This invention relates to the technical field of measuring the relative change in Rayleigh wave attenuation during earthquakes, specifically a three-dimensional Rayleigh wave attenuation relative change imaging method using stable landmasses as reference values. Background Technology

[0002] Earthquakes are frequent in my country, and earthquake prediction is a crucial aspect of earthquake prevention and mitigation. The key to earthquake prediction lies in accurately identifying precursory anomalies. Before an earthquake occurs, the physical properties of the Earth's crustal medium in the source region (such as stress state and fluid content) change, thereby altering the attenuation characteristics of seismic waves. Therefore, monitoring the relative changes in the attenuation coefficient of the crustal medium holds promise for providing important information for identifying earthquake precursors.

[0003] Traditional methods for measuring seismic wave attenuation changes typically involve directly comparing the attenuation coefficients calculated after two earthquake events. However, this method has inherent limitations: the calculation of the attenuation coefficient heavily relies on the accuracy of source parameters (such as source spectrum and radiation pattern) and propagation path effects (such as geometric diffusion). In practical applications, even when two earthquakes are located close to each other, differences in their source mechanisms and precise paths can be introduced into the attenuation coefficient calculation, creating a significant source of error. This results in a low signal-to-noise ratio in the final attenuation change measurement, making it difficult to effectively identify subtle precursor signals truly caused by changes in the physical properties of the medium.

[0004] Furthermore, in order to achieve spatial monitoring, traditional methods require tedious calculations or complex tomographic inversions of a large amount of station and seismic event data, resulting in low computational efficiency and making it difficult to achieve rapid and dynamic monitoring applications.

[0005] Therefore, there is an urgent need in this field for a seismic wave attenuation change imaging method that can effectively separate the source and path effects and directly extract the medium property change signal, so as to improve the monitoring capability of crustal medium state changes and the identification accuracy of earthquake precursor anomalies. Summary of the Invention

[0006] To address the aforementioned problems, specifically those raised in the background section, this invention proposes a three-dimensional Rayleigh wave attenuation relative change imaging method using stable land parcels as reference values, comprising the following steps: S10. Seismic data preprocessing: Obtain the Rayleigh wave amplitude spectra of the two seismic events recorded by seismic stations in the study area and the reference area, respectively; S20. Amplitude ratio calculation: Calculate the Rayleigh wave amplitude ratio for the two earthquake events at the seismic stations in the study area and the reference area, respectively. S30. Calculation of relative attenuation change: The amplitude ratio of the seismic stations in the study area is corrected using the amplitude ratio of the seismic stations in the reference area to eliminate source and common path effects, and the relative change of Rayleigh wave attenuation coefficient of the medium in the study area between the two seismic events is calculated. ; S40, 3D Imaging: Based on the relative changes in attenuation coefficients calculated from multiple stations at multiple frequencies, spatial interpolation is first performed on the data for each single frequency to generate a two-dimensional planar distribution map; then, the two-dimensional maps of different frequencies are combined according to their corresponding depths to finally generate a three-dimensional spatial distribution image.

[0007] A further provision of the present invention is that the S10 seismic data preprocessing includes: S11. Determine the study area and reference area: Delineate the study area according to the monitoring target; Based on the geological structural stability of the surrounding area of ​​the study area, select a stable geological block as the reference area; S12. Data preprocessing and instrument response removal: Select two earthquake events with similar epicentral distances, and download the waveform data of the earthquake events recorded by stations in the study area and reference area respectively; preprocess the waveform data to remove instrument response; S13. Waveform Extraction: For the waveform data of each earthquake event, a time window containing the main earthquake signal is determined, and the waveform data within the time window is extracted, starting from the initial arrival time of the earthquake wave. S14. Calculate the amplitude spectrum: Perform a fast Fourier transform on the truncated waveform data to obtain the corresponding amplitude spectrum; extract the amplitude value of the spectrum of each earthquake event within the predetermined Rayleigh wave dominant frequency range, thereby obtaining the amplitude dataset of each station corresponding to each earthquake event.

[0008] A further provision of the present invention is as follows: In step S20, the Rayleigh wave amplitude ratios of the two earthquake events observed by seismic stations in the study area and the reference area are calculated respectively, wherein the expression for the amplitude ratio of the two earthquakes observed by seismic stations in the study area is:

[0009] Among them, the seismic stations in the study area used The two earthquakes are indicated by subscripts 1 and 2, where S represents the source spectrum, D represents the source radiation pattern, and G represents the geometric diffusion. For frequency, The attenuation coefficient is... The average great circle ray path length in the study area; The expression for the amplitude ratio of two earthquakes observed by seismic stations in the reference area is:

[0010] Among them, the reference area seismic stations used The two earthquakes are indicated by subscripts 1 and 2, where S represents the source spectrum, D represents the source radiation pattern, and G represents the geometric diffusion. For frequency, The attenuation coefficient is... The average great circle ray path length is the reference area; for the reference area, when it contains multiple stations, Take the average value; when the earthquake is far from the study area and the stable block is close to the study area, it is considered that... .

[0011] A further provision of the present invention is that: S30, the correction of the amplitude ratio of the seismic stations in the study area using the amplitude ratio of the seismic stations in the reference area, specifically includes: adjusting the amplitude ratio of the seismic stations in the reference area... Substitute the amplitude ratio of seismic stations in the study area The expression, through formula

[0012] The relative change in Rayleigh wave attenuation coefficient in the study area was calculated. .

[0013] A further provision of the present invention is that the reference area is a stable land parcel, and the relative change in its attenuation coefficient is considered to be zero, i.e. .

[0014] A further provision of the present invention is that during the three-dimensional imaging process in step S40, two-dimensional planar distribution maps of different frequencies are combined according to the relationship between the wavelength of their corresponding Rayleigh waves and the depth of the medium.

[0015] A further provision of the present invention is that, during the three-dimensional imaging process in step S40, the spatially interpolated two-dimensional distribution map is subjected to spatial filtering to smooth out anomalies and enhance meaningful variation features.

[0016] The beneficial technical effects of this invention are as follows: 1. Effectively eliminates source and path effects, improving the signal-to-noise ratio. Traditional methods heavily rely on the accuracy of source parameters and path effects when calculating attenuation coefficient changes, which easily introduces errors. This invention introduces a stable block as a reference area, utilizing its stable physical properties. In step S30, the amplitude ratio of the reference area is used as a benchmark and substituted into the amplitude ratio expression of the study area through the amplitude ratio correction method. This effectively removes common path effects such as source spectrum, radiation pattern, and geometric diffusion, significantly improving the signal-to-noise ratio of attenuation changes and making minute medium change signals more prominent.

[0017] 2. Achieving high-precision, dynamic monitoring of medium changes: Through the reference area correction mechanism introduced in step S30, this invention can directly extract the relative change in attenuation coefficient caused by changes in the physical properties of the medium. This method avoids the systematic errors caused by differences in seismic source and path in traditional methods. It is particularly suitable for long-term, dynamic monitoring of crustal medium conditions, providing a more reliable data foundation for earthquake precursor identification.

[0018] 3. Simplified calculation process and improved efficiency: Traditional methods often require complex dual-station cross-correlation or tomographic inversion, resulting in high computational load and low efficiency. This invention significantly simplifies the data processing flow and reduces computational complexity by directly comparing the amplitude ratio with a reference zone benchmark, making the method more suitable for rapid, large-scale real-time or near-real-time monitoring applications. By spatially interpolating and combining attenuation changes at multiple frequencies, a three-dimensional attenuation change distribution image is generated, which can intuitively display the changes in medium properties in three-dimensional space, helping to identify anomalous areas and their trends with depth, providing richer spatial information for earthquake precursor analysis. Attached Figure Description

[0019] Figure 1 A schematic diagram of the method of this application is shown.

[0020] Figure 2 A schematic diagram showing the tectonic environment and distribution of observation stations in Yangbi is provided.

[0021] Figure 3 The study area is shown in a bar chart of daily precipitation and a curve showing the variation of well water temperature.

[0022] Figure 4 The Rayleigh wave attenuation curve near the epicenter of the Yangbi earthquake is shown.

[0023] Figure 5 The diagram shows the identification of abnormal water temperature in the observation wells before the Yangbi earthquake.

[0024] Figure 6 Three-dimensional imaging of Rayleigh wave attenuation in the study area is shown. Detailed Implementation

[0025] Preferred embodiments of the present invention will now be described with reference to the accompanying drawings. Those skilled in the art should understand that these embodiments are merely illustrative of the technical principles of the present invention and are not intended to limit the scope of protection of the present invention.

[0026] This invention proposes a three-dimensional Rayleigh wave attenuation relative change imaging method using stable land parcels as reference values, comprising the following steps: S10. Seismic data preprocessing: Obtain the Rayleigh wave amplitude spectra of the two seismic events recorded by seismic stations in the study area and the reference area, respectively; Specifically, this includes: S11, determining the study area and reference area: delineating the study area based on the monitoring target; and selecting stable geological blocks as reference areas based on the geological structural stability of the surrounding area. S12. Data preprocessing and instrument response removal: Select two earthquake events with similar epicentral distances, and download the waveform data of the earthquake events recorded by stations in the study area and reference area respectively; preprocess the waveform data to remove instrument response and eliminate possible differences in gain between different models of seismographs; S13. Waveform Extraction: For the waveform data of each earthquake event, determine the time window containing the main earthquake signal, and extract the waveform data within the time window starting from the initial arrival time of the earthquake wave. S14. Calculate the amplitude spectrum: Perform a fast Fourier transform on the truncated waveform data to obtain the corresponding amplitude spectrum; extract the amplitude value of the spectrum of each earthquake event within the predetermined Rayleigh wave dominant frequency range, thereby obtaining the amplitude dataset of each station corresponding to each earthquake event.

[0027] S20. Amplitude ratio calculation: Calculate the Rayleigh wave amplitude ratio for the two earthquake events at the seismic stations in the study area and the reference area, respectively. Seismic waves travel a certain distance in a medium, and the attenuation performance of the waves is characterized by the change in their amplitude. Let be the attenuation coefficient, when The medium is perfectly elastic. A larger value indicates greater internal friction within the medium. If the origin of the coordinate system is the excitation point, the generated amplitude is... The seismic waves propagate to a distance The amplitude decays at time t= ,but After rearranging the terms and taking the natural logarithm, it can be transformed into: Therefore, it can be seen that the traditional calculation method is to calculate the average attenuation coefficient of a certain section of the medium. The method for obtaining the change in attenuation is to take the difference between two calculation results. The subscripts 1 and 2 represent two attenuation coefficient measurement experiments, and the relative change in attenuation coefficient is... .

[0028] In the Earth model, Figure 1 In the diagram, E (a pentagram) represents an earthquake; A, B, C, D, and F (triangles) represent seismic stations; M (shaded area) represents a seismic station in the reference zone of a solid block; subscripts 1 and 2 indicate the observation results of two consecutive earthquakes; and the straight lines represent the propagation paths of seismic waves. The laboratory scheme is as follows: Figure 1As shown in Figure a, two earthquakes located close to each other are used as hypocenters. Measurement stations A, B, and C are approximately aligned with the two earthquakes E1 and E2 (generally, the angle should be within 5°), allowing for the calculation of attenuation changes between stations. That is, through... , , , and Approximate estimate of the decay change at point A Earthquakes located close to each other have a low probability of occurrence; most earthquakes occur randomly in Taiwan and neighboring regions such as Japan, the Philippines, and Indonesia. When two earthquakes (E1 and E2) are significantly apart, such as... Figure 1 As shown in b. Then through , , , , and Approximate estimate of the decay change at point A This traditional method is computationally intensive and requires significant calculations. First, it necessitates locating stations along the earthquake's propagation path based on its azimuth. Then, it calculates the attenuation changes between multiple stations and approximates the attenuation change at a specific observation point. This invention proposes a method that uses stable ground features to correct the attenuation coefficient, enabling rapid calculation of attenuation changes at observation points (e.g.,...). Figure 1 As shown in c), this reduces the amount of computation.

[0029] In the Earth model, the spectral function A of the seismic wave received by a station at an epicentral distance of R ( ) expressed as Where S is the source spectrum, D is the source radiation pattern, G is the geometric diffusion, and L is the site response. For frequency, The attenuation coefficient is... The path length of the great circle ray; This application selects two earthquake events with epicenter locations as close as possible to minimize changes in epicentral distance R and geometric spread G. The two earthquakes are distinguished by subscripts 1 and 2. , , , The site response parameters for fixed stations remain unchanged. .

[0030] The expression for the amplitude ratio of two earthquakes observed by seismic stations in the study area is as follows:

[0031] Among them, the seismic stations in the study area used The two earthquakes are indicated by subscripts 1 and 2, where S represents the source spectrum, D represents the source radiation pattern, and G represents the geometric diffusion. For frequency, The attenuation coefficient is... The average great circle ray path length in the study area; The expression for the amplitude ratio of two earthquakes observed by seismic stations in the reference area is:

[0032] Among them, the reference area seismic stations used The two earthquakes are indicated by subscripts 1 and 2, where S represents the source spectrum, D represents the source radiation pattern, and G represents the geometric diffusion. For frequency, The attenuation coefficient is... The average great circle ray path length is the reference area; for the reference area, when it contains multiple stations, Take the average value; when the earthquake is far from the study area and the stable block is close to the study area, it is considered that... .

[0033] S30. Calculation of Attenuation Change: The amplitude ratio of the seismic stations in the study area is corrected using the amplitude ratio of the seismic stations in the reference area to eliminate source and common path effects. The relative change of Rayleigh wave attenuation coefficient of the medium where the stations in the study area are located between the two seismic events is calculated. ; Specifically, this includes: stable land parcels have relatively high rigidity; when the tectonic stress field changes, the physical properties of stable land parcels change little, and can be considered as... The variation in the attenuation coefficient of stable blocks can be used as a reference benchmark. The expression for the amplitude ratio of two earthquakes observed by seismic stations in the reference area can be simplified to:

[0034] With similar epicentral distances, the geometric spread term of seismic waves can be approximated. ,Right now Substituting the amplitude ratio expression of two earthquakes observed by seismic stations in the study area, and rearranging the terms and taking the natural logarithm of both sides, it can be transformed into...

[0035] The relative change in Rayleigh wave attenuation coefficient in the study area was calculated. Its physical meaning is the average attenuation change per unit length along the great circle propagation path from the earthquake source to the observation station. S40, 3D Imaging: Based on the relative changes in attenuation coefficients calculated from multiple stations at multiple frequencies, spatial interpolation is first performed on the data for each single frequency to generate a two-dimensional planar distribution map; then, the two-dimensional maps of different frequencies are combined according to their corresponding depths to finally generate a three-dimensional spatial distribution image.

[0036] In the three-dimensional imaging process of step S40, two-dimensional planar distribution maps of different frequencies are combined according to the relationship between the wavelength of their corresponding Rayleigh waves and the depth of the medium. The process also includes spatial filtering of the spatially interpolated two-dimensional distribution maps to smooth out anomalies and enhance meaningful variation features.

[0037] Application examples: Groundwater anomalies include rapid changes in water level and temperature in observation wells, as well as variations in the levels of radon, mercury, carbon dioxide, and sulfides in the well water. Prior to the 1975 Haicheng earthquake, numerous groundwater anomalies were detected near the epicenter, providing a scientific basis for the accurate prediction of the Haicheng earthquake. However, many factors can cause groundwater anomalies, and numerous earthquake cases demonstrate that identifying the causes of groundwater anomalies, especially those caused by increased tectonic stress, is extremely difficult.

[0038] Taking the 2021 Yangbi earthquake in Yunnan as an example, this paper introduces the difficulties in identifying the causes of groundwater anomalies and the role of this invention patent in promoting this work.

[0039] On May 21, 2021, an earthquake occurred in Yangbi County, Yunnan Province. A magnitude 6.4 earthquake occurred in the study area (95°–107°E, 21°–31°N), which is divided into the South China Block, Sichuan Basin Block, Songpan-Ganzi Block, Western Sichuan Block, Central Yunnan Block, Lanping-Simao Block, Baoshan Block, and the eastern part of the Qinghai-Tibet Plateau. The distribution of seismic stations within the study area is shown in [reference needed]. Figure 2 The diagram shows the following fault zones: F1 Honghe Fault Zone, F2 Jinshajiang Fault, F3 Xianshuihe Fault, F4 Anninghe-Zemuhe Fault, F5 Xiaojiang Fault, F6 Lijiang-Xiaojinhe Fault, F7 Longmenshan Fault Zone, and F8 Lancangjiang Fault.

[0040] The following is a bar chart showing the daily precipitation values ​​from January 1, 2020 to December 31, 2021 at the seismic monitoring stations in Qujiang and Yaoan within the study area: Figure 3 As shown in Figure a, precipitation in the study area exhibits strong regularity, with the maximum precipitation concentrated in summer and relatively less precipitation in winter. The well water temperature curves of three wells (numbered Kunming, Jiangchuan, and Mile) in the study area during the same period are shown in Figure a. Figure 3As shown in b, c, and d, the annual variation pattern of the water temperature curve in the Kunming well is obvious. The high water temperature value corresponds to the peak precipitation value with a lag of half a month, reflecting that surface water gradually infiltrates into the Kunming well through multiple strata (which may include high-temperature strata), causing its water temperature to show seasonal increases and decreases. The water temperature in the Jiangchuan well shows a steady and slow increase, with no obvious annual variation pattern. The Mile well has no obvious correlation with precipitation, but experienced several large fluctuations before the Yangbi earthquake. The significantly different annual variation patterns indicate that the three observation wells are located in different hydrogeological environments.

[0041] Well water temperature is affected not only by precipitation but also by solid tides, water level, human groundwater extraction activities, and the effects of underground tectonic stress. When the changes in well water temperature are complex, it is difficult to identify which stage is caused by enhanced tectonic stress leading to earthquake precursor anomalies. Without accurately identifying earthquake precursor anomalies, it is difficult to make accurate earthquake predictions.

[0042] Using the earthquakes in Table 1 as hypocenters, and based on the method of this patent, the hypocenters near the epicenter of the Yangbi earthquake were calculated, with the Sichuan Basin (a stable block) as the reference area. Figure 2 The relative change curve of Rayleigh wave attenuation coefficient of seismic stations within the area indicated by the middle circle. Figure 4 ).from Figure 4 As can be seen, a positive anomaly in Rayleigh wave attenuation occurred in early February 2021 prior to the Yangbi earthquake. Due to the slow thermal conductivity of crustal rocks and the linear decrease in the arrival time of seismic waves near the epicenter, it is believed that the main reason for the dramatic fluctuations in Rayleigh wave attenuation in early February 2021 was the upwelling of fluids within the crust, leading to an increase in fluid content in the middle and upper crust.

[0043] Table 1 List of earthquake events that served as epicenters № Date (UTC) Latitude (°) Longitude (°) Magnitude Depth (km) Region 1 March 18, 2020 -11.04 115.23 6.3 10 BALI, INDONESIA 2 August 5, 2020 -9.74 119.08 5.5 10 SUMBA, INDONESIA 3 January 14, 2021 -2.94 118.93 6.2 20 SULAWESI, INDONESIA 4 February 10, 2021 -5.62 101.73 6.2 10 SUMATRA, INDONESIA 5 March 8, 2021 -4.34 101.1 5.5 10 SUMATRA, INDONESIA 6 April 10, 2021 -8.54 112.51 6.0 89 JAVA, INDONESIA 7 21 May 2021 -8.35 112.37 5.7 117 JAVA, INDONESIA 8 July 26, 2021 -0.65 121.94 6.2 10 SULAWESI, INDONESIA 9 December 6, 2022 -10.72 113.39 5.9 18 JAVA, INDONESIA Data source: European-Mediterranean Seismological Centre (www.emsc-csem.org).

[0044] According to the source information in Table 1, the attenuation changes occurred between January 14, 2021, and March 8, 2021. (Excerpt) Figure 3 Water temperature variation curves of the three wells in the middle section from January 1, 2021 to March 31, 2021 (as shown in the figure). Figure 5 It was discovered that all three wells exhibited abnormally high water temperatures from late January to mid-February 2021. Therefore, it is believed that... Figure 5 The rapid rise in water temperature simultaneously observed in three wells located in different hydrogeological environments within a medium-gray columnar background period is a typical precursory anomaly to an earthquake.

[0045] To further investigate the extent and scale of crustal fluid upwelling before the Yangbi earthquake, the earthquakes in Table 2 were used as hypocenters to calculate the relative changes in Rayleigh wave attenuation coefficients at seismic stations in the study area.

[0046] Table 2 List of earthquake events with epicenters № Date (UTC) Latitude (°) Longitude (°) Magnitude Depth (km) Region 1 January 14, 2021 -2.94 118.93 6.2 20 SULAWESI, INDONESIA 2 February 10, 2021 -5.62 101.73 6.2 10 SUMATRA, INDONESIA 3 January 6, 2021 -29.11 -177.06 6.2 40 KERMADECISLANDS, NEWZEALAND 4 February 10, 2021 -23.2 171.59 6.4 10 LOYALTY ISLANDS Data source: European-Mediterranean Seismological Centre (www.emsc-csem.org).

[0047] Interpolation is performed between observation points to generate a spatial distribution map of the relative attenuation of Rayleigh waves at different dominant frequencies. Figure 1-6 ).in Figure 6 The dominant frequency of wave a is 0.05 Hz. Based on the correlation between Rayleigh wave dominant frequency and depth, 0.05 Hz corresponds to a depth of approximately 18-23 km, which in the study area corresponds to the depth range from the bottom of the upper crust to the top of the middle crust. Figure 6 The frequency of b is 0.043Hz, corresponding to the middle crustal depth range. Figure 6 The darker the color, the greater the increase in Rayleigh wave attenuation, indicating a greater increase in fluid content. As deep hydrothermal fluids rise, the well water temperature in surface observation wells may rise rapidly, forming typical precursory anomalies of earthquakes. Therefore, this patent is of great significance for studying changes in the physical properties of the Earth's crust.

[0048] Furthermore, it should be noted that, in the description of this invention, unless otherwise explicitly specified and limited, the terms "installation," "connection," and "linking" should be interpreted broadly. For example, they can refer to a fixed connection, a detachable connection, or an integral connection; they can refer to a mechanical connection or an electrical connection; they can refer to a direct connection or an indirect connection through an intermediate medium; and they can refer to the internal connection of two components. Those skilled in the art can understand the specific meaning of the above terms in this invention according to the specific circumstances.

[0049] The term "comprising" or any other similar term is intended to cover non-exclusive inclusion, such that a process, article, or apparatus / device that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to those processes, articles, or apparatus / devices.

[0050] The technical solution of the present invention has been described above with reference to the preferred embodiments shown in the accompanying drawings. However, it will be readily understood by those skilled in the art that the scope of protection of the present invention is obviously not limited to these specific embodiments. Without departing from the principles of the present invention, those skilled in the art can make equivalent changes or substitutions to the relevant technical features, and the technical solutions after these changes or substitutions will all fall within the scope of protection of the present invention.

Claims

1. A method for imaging relative changes in attenuation of three-dimensional Rayleigh waves using stable block as a reference value, characterized in that: The method comprises the following steps: S10, seismic data preprocessing: obtaining Rayleigh wave amplitude spectrum recorded by seismic stations in the study area and reference area for two earthquake events respectively; S20, amplitude ratio calculation: calculating Rayleigh wave amplitude ratio of two earthquake events for seismic stations in the study area and reference area respectively; S30, attenuation change calculation: the amplitude ratio of the reference area seismic station is used to correct the amplitude ratio of the study area seismic station to eliminate the source and common path effects, and the relative change of the Rayleigh wave attenuation coefficient of the medium where the study area station is located between two seismic events is calculated ; S40, three-dimensional imaging: based on the relative change amount of attenuation coefficient calculated by multiple stations at multiple frequencies, firstly, spatial interpolation is carried out for each single frequency data to generate a two-dimensional plane distribution map; then, two-dimensional maps of different frequencies are combined according to their corresponding depths to finally generate a three-dimensional spatial distribution image.

2. The method of claim 1, wherein the method is characterized by: The S10 seismic data preprocessing comprises: S11, determining the study area and the reference area: determining the study area according to the monitoring target; selecting a stable geological block as the reference area based on the stability of the geological structure around the study area; S12, data preprocessing and instrument response removal: selecting two earthquake events with similar epicentral distance, downloading waveform data recorded by stations in the study area and reference area for the earthquake events respectively; preprocessing the waveform data and removing instrument response; S13, waveform interception: for the waveform data of each earthquake event, determining a time window containing the main signal of the earthquake, and taking the arrival time of the seismic wave as the starting point to intercept the waveform data in the time window; S14, amplitude spectrum calculation: performing fast Fourier transform on the intercepted waveform data to obtain the corresponding amplitude spectrum; in the predetermined Rayleigh wave main frequency range, extracting the amplitude value of the spectrum of each earthquake event, thereby obtaining the amplitude data set corresponding to each earthquake event for each station.

3. The method of claim 1, wherein the method is characterized by: The S20, respectively, calculates the Rayleigh wave amplitude ratio of the two earthquake events of the seismic station in the study area and the reference area, wherein the amplitude ratio expression of the two earthquakes observed by the seismic station in the study area is: Wherein the seismic station in the study area is represented by , the two earthquakes are represented by subscripts 1 and 2, S is the source spectrum, D is the source radiation pattern, G is the geometric diffusion, is the frequency, is the attenuation coefficient, is the average great circle ray path length in the study area; the amplitude ratio expression of the two earthquakes observed by the seismic station in the reference area is: Wherein the seismic station in the reference area is represented by , the two earthquakes are represented by subscripts 1 and 2, S is the source spectrum, D is the source radiation pattern, G is the geometric diffusion, is the frequency, is the attenuation coefficient, is the average great circle ray path length in the reference area; for the reference area, when multiple stations are included, Take the average value; when the earthquake distance is far from the study area and the stable block is close to the study area, it is considered that .

4. The method of claim 1, wherein the method is characterized by: S30, the amplitude ratio of the research area seismic station is corrected by using the amplitude ratio of the reference area seismic station, specifically including: substituting the expression of the reference area seismic station amplitude ratio into the expression of the research area seismic station amplitude ratio , and calculating the relative variation of the Rayleigh wave attenuation coefficient of the research area by the formula .

5. The method of claim 4, wherein the method is characterized by: The reference region is a stable block, whose relative variation of attenuation coefficient is considered as zero, i.e. .

6. The method of claim 1, wherein the method is characterized by: In the three-dimensional imaging process of step S40, the two-dimensional plane distribution maps of different frequencies are combined according to the relationship between the wavelength of the corresponding Rayleigh wave and the medium depth.

7. The method of claim 1, wherein the method is characterized by: In the three-dimensional imaging process of step S40, the two-dimensional distribution map generated by spatial interpolation is also subjected to spatial filtering processing to smooth abnormal points and enhance meaningful change characteristics.

Images