CN106646609A - Multi-scan microseism multi-parameter joint quick inversion method - Google Patents

Abstract

The invention discloses a multi-scan microseismic multi-parameter joint fast inversion method, the method steps are: (1) data acquisition; (2) preprocessing; (3) dynamic correction and stacking; (4) microseismic positioning and source Mechanism joint inversion; (5) Carry out the initial scan of space and angle with large steps; (6) Judge whether there is an event, if not, do not output; if so, perform secondary scanning with small space steps and large angle steps Scanning; (7) Preliminarily determine the spatial coordinates of the seismic source and reduce the scanning range; (8) Perform three scans with small steps in space and angle; (9) Output the results. The joint fast inversion method of the present invention can identify the shear source with polarity change, while the traditional seismic source scanning method cannot identify the change of seismic source polarity; the present invention needs to scan tens of millions of gathers, and the calculation amount is huge. Multiple scans with large and small steps are optimized for acceleration, which can effectively save calculation time while ensuring positioning accuracy.

Description

Multi-scan microseismic multi-parameter joint fast inversion method

technical field

The invention relates to the technical field of microseismic monitoring of hydraulic fracturing in shale gas development, in particular to the technology of microseismic positioning and focal mechanism inversion.

Background technique

Microseismic monitoring is currently one of the most accurate, timely, and informative monitoring methods in reservoir fracturing, and the primary task of microseismic monitoring is to determine the source location, time and intensity of microseismic events.

Microseismic monitoring mainly includes surface monitoring and borehole monitoring. Downhole monitoring and positioning technology is relatively mature, but the number of monitoring wells, spatial distribution and number of geophones are very limited, so its lateral resolution is low, and Monitoring costs are also high. On the other hand, surface microseismic monitoring technology has fewer restrictions on the lateral distribution range and density, can obtain higher lateral resolution, and has lower monitoring cost. Therefore, the present invention mainly studies the ground microseismic monitoring system.

With the maturity of surface microseismic monitoring technology and the improvement of positioning accuracy, the interpretation of the distribution of microseismic events combined with the focal mechanism will become the focus of future work. The focal mechanism of an earthquake can effectively explain the physical properties of the source and can indicate the regional or local stress and strain distribution. Combined with the surrounding geological conditions, the fracturing effect can be systematically and accurately evaluated.

Therefore, it is necessary to propose a multi-scan microseismic multi-parameter joint fast inversion method for the above problems.

Contents of the invention

Aiming at the deficiencies in the above-mentioned prior art, the present invention proposes a multi-scanning microseismic multi-parameter joint fast inversion method, which can quickly and accurately explain the orientation and shape of fractures and achieve the purpose of evaluating fracturing effects.

The multi-scan microseismic multi-parameter joint fast inversion method, the method steps are: (1) Input the monitoring system to collect seismic data; (2) Perform band-pass filtering and static correction on the original data collected in step (1) Preprocessing; (3) Perform dynamic correction on the preprocessed seismic records with a virtual source point as the source; (4) Use the virtual source as the source to add changes in direction, dip, and slip angle to perform polarity correction and stacking (5) Carry out joint inversion of microseismic location and focal mechanism; (6) Perform a large-step initial scan of the focal space and focal mechanism; (7) Determine whether there is Microseismic events, if not, then do not output; if so, perform a second scan with a small step in space and a large step in angle; (8) preliminarily determine the spatial coordinates and focal mechanism of the seismic source, and reduce the scanning range; (9) scan the source The spatial and focal mechanisms perform three scans with small steps; (10) output the results.

Preferably, wherein, the method for joint inversion of microseismic location and focal mechanism is an improved source scanning algorithm, which

where _uj is the observed seismic record, N is the total number of stations, τ is the time of the hypocenter i, t _ij is the travel time from the virtual hypocenter i to the station j, and E is the space-time function of the hypocenter mechanism;

m _ij ＝f(φ _i ,δ _i ,λ _i ,α _ij ,β _ij ) (2)

Among them, φ _i , δ _i and λ _i are the strike, dip and slip angle of virtual source point i respectively, α _ij is the azimuth of station j relative to virtual source point i, and β _ij is the azimuth of station j relative to virtual source point i The angle from the source, m _ij reflects the influence of the hypocenter mechanism of the virtual hypocenter point i on the P wave polarity of the seismic record of the station j;

Where ( _xi , y _i , z _i ) is the source coordinates, (x _j , y _j , z _j ) is the station coordinates, w _ij is the weight function of station j relative to the virtual source point i, the weight function and The virtual source point is related to the geometric position of the monitoring station. The closer the station is to the virtual source point, the higher the weight value.

The calculation process of the above formula is to convert the original seismic records into the stacked gathers of each virtual source, and then in the time-space domain of the source mechanism, the long-short time-window ratio is calculated for each sampling point one by one, and the "brightness" value R of each point is calculated;

R _k =STA _k /LTA _k (4)

In the formula, R _k is the "brightness" of the kth sampling point, STA _k is the root mean square amplitude of the short time window after the kth sampling point, N _sta is the number of sampling points in the short time window; LTA _k is the kth sampling point The root mean square amplitude of the long time window before the point, N _lta is the number of sampling points in the long time window, and A(k) is the amplitude of the kth sampling point.

In the time-space domain of the focal mechanism, when the "brightness" R of a certain sampling point is greater than a certain threshold, it is regarded as a seismic event. The spatial position and focal mechanism corresponding to the sampling point are the inversion results of the event, and the corresponding moment is the moment of the earthquake.

In the time-space domain of the focal mechanism, it takes a very long time to calculate and scan tens of millions of gathers point by point. Even a multi-core multi-threaded workstation cannot meet the needs of real-time positioning and inversion. Therefore, the main purpose of the present invention is to optimize the acceleration of this method. For the spatial location and focal mechanism of the seismic source, first scan the spatial location and focal mechanism of the seismic source with a larger step size; if there is a microseismic event, then scan the source location with a smaller spatial step The spatial position is finely positioned; then, according to the fine positioning results, the scanning range of the source space is narrowed, and the source space and the source mechanism are scanned three times with a smaller step size. In this way, the calculation time can be effectively saved, and the positioning accuracy can also be ensured.

Due to the adoption of the above technical scheme, the joint fast inversion method of the present invention can identify the shear source with polarity change, while the traditional source scanning method cannot identify the change of source polarity; the present invention needs to scan tens of millions of gathers , the amount of calculation is huge, and the optimization and acceleration of multiple scans with large and small steps can effectively save calculation time while ensuring positioning accuracy, and has strong practicability.

Description of drawings

Figure 1 is the basic flow chart of the combined rapid inversion of microseismic location and focal mechanism of multiple scans;

Figure 2 is a schematic diagram of the ground monitoring system and source space;

Figure 3 is a seismic record with a signal-to-noise ratio of 0.5;

Fig. 4 is the seismic record after preprocessing;

Figure 5 is the seismic record after motion correction;

Figure 6 is the virtual seismic source superposition gather;

Figure 7 is the polarity distribution of the ground monitoring system relative to the seismic source;

Figure 8 is the P-axis inversion results of data with different SNRs using small step size scanning in both space and angle;

Figure 9 is the P-axis inversion results of data with different signal-to-noise ratios for scanning with a small step size in space and a large step size in angle;

Figure 10 is the P-axis inversion results of data with different signal-to-noise ratios scanned with a large spatial step and a small angular step;

Figure 11 is the P-axis inversion results of data with different SNRs using multiple scans;

Figure 12 is a time-consuming histogram for the inversion of different scanning methods in Table 1;

Fig. 13 is the inversion results of different focal mechanisms when joint fast inversion is adopted;

Fig. 14 is the P-axis inversion result of the column source when the joint fast inversion is adopted;

Fig. 15 is the P-axis inversion result of the column source with a signal-to-noise ratio of 0.5 when the joint fast inversion is adopted.

detailed description

The embodiments of the present invention will be described in detail below with reference to the accompanying drawings, but the present invention can be implemented in many different ways defined and covered by the claims.

Implementation case one:

As shown in Fig. 1 and Fig. 2 to Fig. 15, the multi-scanning microseismic multi-parameter joint fast inversion method has the following steps: (1) input the seismic data into the monitoring system; (2) combine the steps (1 ) band-pass filtering and static correction preprocessing for the original data collected; (3) dynamic correction of the preprocessed seismic records with a virtual source point as the source; (4) adding trend, After changing the dip angle and slip angle, polarity correction and superposition processing are carried out to obtain a virtual focal mechanism stacked gather; (5) joint inversion of microseismic location and focal mechanism; (6) focal space and focal mechanism The first scan with a large step; (7) judge whether there is a microseismic event, if not, then do not output; if so, perform a second scan with a small step in space and a large step in angle; (8) initially determine the spatial coordinates of the source and focal mechanism, narrowing the scanning range; (9) performing three small-step scans on the focal space and focal mechanism; (10) outputting the result.

Further, among them, the joint inversion of microseismic location and focal mechanism adopts an improved source scanning algorithm, and its

where _uj is the observed seismic record, N is the total number of stations, τ is the time of the hypocenter i, t _ij is the travel time from the virtual hypocenter i to the station j, and E is the space-time function of the hypocenter mechanism;

m _ij ＝f(φ _i ,δ _i ,λ _i ,α _ij ,β _ij ) (2)

Among them, φ _i , δ _i and λ _i are the strike, dip and slip angle of virtual source point i respectively, α _ij is the azimuth of station j relative to virtual source point i, and β _ij is the azimuth of station j relative to virtual source point i The angle from the source, m _ij reflects the influence of the hypocenter mechanism of the virtual hypocenter point i on the P wave polarity of the seismic record of the station j;

Where ( _xi , y _i , z _i ) is the source coordinates, (x _j , y _j , z _j ) is the station coordinates, w _ij is the weight function of station j relative to the virtual source point i, the weight function and The virtual source point is related to the geometric position of the monitoring station. The closer the station is to the virtual source point, the higher the weight value.

The calculation process of the above formula is to convert the original seismic records into the stacked gathers of each virtual source, and then in the time-space domain of the source mechanism, the long-short time-window ratio is calculated for each sampling point one by one, and the "brightness" value R of each point is calculated;

R _k =STA _k /LTA _k (4)

In the formula, R _k is the "brightness" of the kth sampling point, STA _k is the root mean square amplitude of the short time window after the kth sampling point, N _sta is the number of sampling points in the short time window; LTA _k is the kth sampling point The root mean square amplitude of the long time window before the point, N _lta is the number of sampling points in the long time window, and A(k) is the amplitude of the kth sampling point.

In the time-space domain of the focal mechanism, when the "brightness" R of a certain sampling point is greater than a certain threshold, it is regarded as a seismic event. The spatial position and focal mechanism corresponding to the sampling point are the inversion results of the event, and the corresponding moment is the moment of the earthquake.

In the time-space domain of the focal mechanism, it takes a very long time to calculate and scan tens of millions of gathers point by point. Even a multi-core multi-threaded workstation cannot meet the needs of real-time positioning and inversion. Therefore, the main purpose of the present invention is to optimize the acceleration of this method. For the spatial location and focal mechanism of the seismic source, first scan the spatial location and focal mechanism of the seismic source with a larger step size; if there is a microseismic event, then scan the source location with a smaller spatial step The spatial position is finely positioned; then, according to the fine positioning results, the scanning range of the source space is narrowed, and the source space and the source mechanism are scanned three times with a smaller step size. In this way, the calculation time can be effectively saved, and the positioning accuracy can also be ensured.

Implementation case two:

Combined with Figure 1, the basic process of the joint rapid inversion of microseismic events is described, and the ground monitoring system shown in Figure 7 is adopted, which has a relatively large range of coverage, and the interfaces of polarity changes are basically distributed in the coverage of the geophones Within the range, the change of the polarity of the seismic source can be fully reflected.

Step 1: After the field data collection is completed, preprocessing such as band-pass filtering and static correction is performed on the original data. In the original seismic records with a signal-to-noise ratio of 0.5 in Fig. 3, almost no seismic events can be seen. After (1,5, Figure 4 is obtained after 75,85) Hz bandpass filtering and static correction at a longitudinal wave velocity of 2500m/s. It can be seen that there are longitudinal waves and shear waves near 700ms in Figure 4, but they are not clear enough. Since the microseismic signals monitored on the surface are often very weak and cannot even be identified from the seismic waveform, filtering and static correction are essential steps.

Step 2, the preprocessed seismic records are dynamically corrected with a virtual source point as the source to obtain Figure 5, and then the traces in Figure 5 are superimposed and averaged to obtain the one in Figure 6, which is the blasting source without the source mechanism The virtual source overlay gather of . In the next few tracks, polarity correction and superposition processing are performed with the virtual seismic source as the seismic source after adding changes in trend, dip angle, and slip angle, and a virtual seismic source stacking gather including the focal mechanism is obtained. This process is the process of converting Figure 4 into Figure 6 through formula (1) and formula (2).

Step 3, for the virtual seismic source superposition gather in Fig. 6, the long and short time window ratio is calculated point by point for each trace, R = STA/LTA, and the R value of each point is calculated. R reflects the microseismic signal in the time and space domain of the focal mechanism The ratio of the rms amplitude relative to the background noise, called "brightness". STA is the short-time window average value, which reflects the root-mean-square amplitude of the microseismic signal; LTA is the long-time window average value, which reflects the root-mean-square amplitude of the background noise. For each point in the virtual source stacked gather, when the "brightness" R is greater than a certain threshold, it is considered a microseismic event, and the spatial position and focal mechanism of the virtual source corresponding to this point are the inversion results of the event. Here we take the threshold value of R as 3.0, the time window length of STA as 0.05 seconds, and the time window length of LTA as 1.0 seconds.

Generally speaking, the scanning range of the source space is in a three-dimensional grid with a side length of 200 meters centered on the perforation, and the angle range of the source mechanism is strike 0-360°, dip angle 0-90°, and slip angle 0-180° . If the spatial step is 10 meters and the angular step is 15°, then the superimposed traces of the virtual source are 20*20*20*24*6*12=13824000 traces, and it takes a long time to scan and calculate trace by trace and point by point. Optimize for acceleration. Therefore, we adopted different scanning schemes and compared their inversion results: traditional scanning with small step size (Figure 8), scanning with small step size in space and large step size in angle (Figure 9), scanning with large step size in space and small step size in angle (Fig. 10), multiple scan scheme (Fig. 11).

For the inversion process in Figure 11, the first scan is performed with a larger step size (spatial step size 40 m, angle step size 60°), and the stacked gathers of the first scan are 5*5*5*6*2*3= 4500 roads. Then scan each channel point by point, in the virtual source grid where the event exists, and then perform a second scan with a smaller space step (10 meters in space and 60° in angle), and the second scan is 4 *4*4*6*2*3=2304 traces, determine the spatial coordinates of the seismic source and the preliminary focal mechanism. According to the focal coordinates and preliminary focal mechanism, the scanning range is appropriately reduced (spatial range 40 meters, strike range 0-120°, dip angle range 0-90°, slip angle range 0-120°), with a smaller space and mechanism step (The space step length is 10 meters, and the angle step length is 15°). Three scans are performed, 4*4*4*8*6*8=24576 traces. Therefore, a total of 4500+2304+24576=31380 traces need to be calculated for the first scan, the second scan and the third scan, which is 13824000-31380=13792620 scans less than before.

Comparing the scan with small step size (Figure 8, 13824000 tracks), the scan with small step size and angle with large step length in space (Figure 9, 288000 tracks) and the scan with small step size with large space step size and angle (Figure 10, 216000 tracks), multiple scans (Figure 11, track 31380) It will save several times or even hundreds of times of time, and the inversion result is also more reasonable. As shown in Figure 11, the positioning errors are all within 20 meters, and most of the included angles of the P axis are around 10°, and the lower the signal-to-noise ratio of the synthetic data, the greater the positioning errors and focal mechanism errors. Overall, the inversion results in Figure 11 are the closest to those in Figure 8. While the positioning error in Figure 10 is above 30 meters, most of the included angles of the P axis are 11°. The positioning error in Figure 9 is about 20 meters, but because the scanning angle step is just the same as the theoretical source mechanism, both are 60 °, so the included angle of the P axis is 0°, which does not have a reasonable reference value.

Therefore, we compared the inversion of the source with a depth of 1000 meters, a strike angle of 32°, a dip angle of 23°, and a slip angle of 72°. Such a source is more general. Synthetic SNR 0.5, sampling rate 4ms, length 2s, segy file size 2.13MB, theoretical source coordinates (0, 0, 1000), source mechanism (32, 23, 72). The specific inversion results of various scanning methods are as follows:

Table 1 Comparison of inversion results of different scanning methods

It can be seen from Table 1 that the joint fast inversion method of multiple scans can effectively save calculation time while ensuring the inversion accuracy for general seismic sources. Combined with Figure 12, it can be seen that the time-consuming of multiple scans is about 500 times less than that of small-step scans, and it is also 5 or 6 times less than the other two scan methods.

Implementation case three:

For this joint fast inversion method, we tested point sources with different source mechanisms and column sources with different signal-to-noise ratios, and the results are as follows:

In Fig. 13, the events of different focal mechanisms at a depth of 1000 m, with a strike of 60°, a dip of 60°, and a slip angle varying from 30° to 150°, combined with the rapid inversion focal mechanism (figure a) and the theoretical focal mechanism (figure b) Basically the same, indicating that this method can also obtain more accurate inversion results for different focal mechanisms.

In Fig. 14, the strike is 60°, the dip angle is 60°, the slip angle is 90°, and the depth is 1000 meters. There are 11 seismic sources arranged at intervals of 50 meters from -250 meters to 250 meters on the y-axis. The included angle between the theoretical P-axis and the inversion P-axis is 0°, the inversion result of the horizontal y-axis is completely consistent with the theoretical coordinates, and the vertical positioning error is also within 20 meters; Figure 15 adds random Noise synthetic data with a signal-to-noise ratio of 0.5, in which the angle of the P-axis is about 10° on average, the inversion results of the horizontal y-axis have an error of 10m at 100m and 150m, and the vertical positioning error is mostly within 30m. It shows that the lower the signal-to-noise ratio, the larger the inversion space error and focal mechanism error, but the error is within a reasonable error range. Therefore, this method can also obtain reasonable inversion results for seismic sources with different signal-to-noise ratios.

To sum up, this multi-scan microseismic location and focal mechanism joint fast inversion method can obtain more accurate inversion results for seismic sources with different signal-to-noise ratios and different focal mechanisms, and the calculation efficiency should also be greatly improved. .

The joint fast inversion method of the present invention can identify the shear source with polarity change, while the traditional seismic source scanning method cannot identify the change of seismic source polarity. Tens of millions of gathers need to be scanned, and the amount of calculation is huge. Optimizing and accelerating long and multiple scans can effectively save computing time while ensuring positioning accuracy, and has strong practicability.

The above descriptions are only preferred embodiments of the present invention, and are not intended to limit the patent scope of the present invention. Any equivalent structure or equivalent process transformation made by using the description of the present invention and the contents of the accompanying drawings, or directly or indirectly used in other related All technical fields are equally included in the scope of patent protection of the present invention.

Claims (2)

1. The multi-scan microseismic multi-parameter joint fast inversion method is characterized in that: the method steps are: (1) Input into the monitoring system to collect seismic data; (2) performing band-pass filtering and static correction preprocessing on the raw data collected in step (1); (3) The preprocessed seismic records are dynamically corrected with a virtual source point as the source; (4) Take the virtual seismic source as the seismic source to add changes in the trend, dip angle, and slip angle, and then perform polarity correction and stacking processing to obtain a virtual seismic source stacking gather including the focal mechanism; (5) Carry out joint inversion of microseismic location and focal mechanism; (6) Perform a large-step initial scan of the source space and focal mechanism; (7) Determine whether there is a microseismic event, if not, then do not output; if so, perform a secondary scan with a small step size in space and a large step size in angle; (8) Preliminarily determine the spatial coordinates and focal mechanism of the seismic source, and narrow the scanning range; (9) Carry out three scans with a small step size on the focal space and focal mechanism; (10) output the result. 2. the microseismic multi-parameter joint rapid inversion method of multiple scans according to claim 1 is characterized in that: wherein, the method for carrying out microseismic location and focal mechanism joint inversion is an improved focal scanning algorithm, which $\mathrm{<span\; class="notranslate">\; E.</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; \&tau;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; N</span>}}\underset{\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; N</span>}}{<span\; class="notranslate">\; \&Sigma;</span>}}{\mathrm{<span\; class="notranslate">\; w</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; j</span>}}{\mathrm{<span\; class="notranslate">\; m</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&tau;</span>}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>$ where _uj is the observed seismic record, N is the total number of stations, τ is the time of the hypocenter i, t _ij is the travel time from the virtual hypocenter i to the station j, and E is the space-time function of the hypocenter mechanism; m _ij ＝f(φ _i ,δ _i ,λ _i ,α _ij ,β _ij ) (2) Among them, φ _i , δ _i and λ _i are the strike, dip and slip angle of virtual source point i respectively, α _ij is the azimuth of station j relative to virtual source point i, and β _ij is the azimuth of station j relative to virtual source point i The angle from the source, m _ij reflects the influence of the hypocenter mechanism of the virtual hypocenter point i on the P wave polarity of the seismic record of the station j; $\begin{array}{c}{\mathrm{<span\; class="notranslate">\; w</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; )</span>\\ <span\; class="notranslate">\; =</span>{<span\; class="notranslate">\; (</span>\sqrt{{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; )</span>}^{<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}\end{array}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 3</span>}<span\; class="notranslate">\; )</span>$ Where ( _xi , y _i , z _i ) is the source coordinates, (x _j , y _j , z _j ) is the station coordinates, w _ij is the weight function of station j relative to the virtual source point i, the weight function and The virtual source point is related to the geometric position of the monitoring station. The closer the station is to the virtual source point, the higher the weight value. The calculation process of the above formula is to convert the original seismic records into the stacked gathers of each virtual source, and then in the time-space domain of the source mechanism, the long-short time-window ratio is calculated for each sampling point one by one, and the "brightness" value R of each point is calculated; R _k =STA _k /LTA _k (4) ${\mathrm{<span\; class="notranslate">\; STA</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; =</span>\sqrt{\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{{\mathrm{<span\; class="notranslate">\; N</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}\mathrm{<span\; class="notranslate">\; t</span>}\mathrm{<span\; class="notranslate">\; a</span>}}}\underset{\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{{\mathrm{<span\; class="notranslate">\; N</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}\mathrm{<span\; class="notranslate">\; t</span>}\mathrm{<span\; class="notranslate">\; a</span>}}}{<span\; class="notranslate">\; \&Sigma;</span>}}{\mathrm{<span\; class="notranslate">\; A</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; )</span>}$ ${\mathrm{<span\; class="notranslate">\; LTA</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; =</span>\sqrt{\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{{\mathrm{<span\; class="notranslate">\; N</span>}}_{\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; t</span>}\mathrm{<span\; class="notranslate">\; a</span>}}}\underset{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{{\mathrm{<span\; class="notranslate">\; N</span>}}_{\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; t</span>}\mathrm{<span\; class="notranslate">\; a</span>}}}{<span\; class="notranslate">\; \&Sigma;</span>}}{\mathrm{<span\; class="notranslate">\; A</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 5</span>}<span\; class="notranslate">\; )</span>$ In the formula, R _k is the "brightness" of the kth sampling point, STA _k is the root mean square amplitude of the short time window after the kth sampling point, N _sta is the number of sampling points in the short time window; LTA _k is the kth sampling point The root mean square amplitude of the long time window before the point, N _lta is the number of sampling points in the long time window, and A(k) is the amplitude of the kth sampling point.