A method for tracking earthquake source location based on coda wave interferometry

Abstract

The purpose of the present invention is to provide a method for tracking the location of an earthquake source based on coda wave interference, comprising the following steps: obtaining time domain signals before and after the change in the location of the earthquake source; signal segmentation; time delay estimation of the change in the location of the earthquake source; distance difference calculation; location estimation of the earthquake source; and tracking the location of the earthquake source. From a theoretical perspective, the present invention combines the characteristics of the coda wave, utilizing the characteristics of the coda wave being sensitive to medium changes and the strong energy of the coda wave signal under long-distance reception, and analyzes the change in the coda wave phase to distinguish the time delays caused by medium changes and source location changes. In particular, when the medium changes in an inhomogeneous medium, the time delay caused by the medium change is removed, and the time delay caused by the change in the location of the earthquake source is accurately and reliably calculated, thereby achieving high-precision tracking of the location of the earthquake source. From an engineering application perspective, the present invention has lower requirements for station distribution, is less dependent on the medium velocity structure, and is more practical.

Description

Technical Field

[0001] The present invention relates to an earthquake detection method, in particular to an earthquake source location tracking method. Background Art

[0002] The focal locations and their variations of small and moderate earthquakes are crucial for studying seismic tectonics, seismicity, fault properties, and the evolution of seismic stress fields. Commonly used methods for tracing earthquake focal locations are primarily based on the processing of direct waves. While direct waves can estimate the time delay between earthquake events, they cannot determine whether the delay is caused by changes in the focal location or by changes in the medium. Furthermore, the accuracy of the tracking results depends on the station distribution and the medium velocity structure. When the detailed medium velocity structure is poorly understood, the paths of received signals differ between stations, and the direct wave may not be the fastest, resulting in suboptimal results for direct-wave-based focal location tracing. In complex geological environments, even with minor perturbations in the propagation path or medium velocity, the direct wave portion remains highly consistent between earthquake events, making them difficult to distinguish. Furthermore, the direct wave signal is transformed into a coda signal as it scatters, rapidly decaying in amplitude while the coda energy gradually increases. In long-distance propagation scenarios, the direct wave energy may even be weaker than the coda energy. For example, the papers "Time-lapse Monitoring with Coda Wave Interferometry" and "Monitoring Velocity Variations in the Crust Using Earthquake Doublets' An Application to the Calaveras Fault, California" show that the coda wave energy of the received signal at the station is stronger than the direct wave. Therefore, existing earthquake source location tracking methods have difficulty distinguishing between earthquake source location changes and changes in the medium structure, and also have difficulty achieving high-precision tracking of earthquake source location changes at long distances.

[0003] In summary, the existing source location tracking methods are difficult to achieve long-distance reception, source location disturbances, and high-precision tracking of the source location when there is an inhomogeneous velocity variation field, and cannot meet the high-precision requirements of earthquake tectonic research. Summary of the Invention

[0004] The purpose of the present invention is to provide a source location tracking method based on coda wave interferometry that can meet the research needs of seismic activity, fault properties, etc.

[0005] The object of the present invention is achieved like this:

[0006] The present invention provides a method for tracking earthquake source position based on coda wave interferometry, which is characterized by:

[0007] (1) Obtaining time domain signals before and after the earthquake source position changes;

[0008] (2) Signal segmentation;

[0009] (3) Time delay estimation of earthquake source location change;

[0010] (4) Distance difference calculation;

[0011] (5) Estimation of earthquake source location;

[0012] (6) Tracking of earthquake source location.

[0013] The present invention may also include:

[0014] 1. Step (1) obtains the time domain signals before and after the earthquake source position changes as follows:

[0015] Determine the initial source position S1(x1,y1) and the receiving point position R i =(X i ,Y i ), 1≤i≤n and medium propagation speed v, obtain the receiving point R i Time domain signal A before and after the change of the earthquake source position i (t) and B i (t).

[0016] 2. The specific method of step (2) signal segmentation is:

[0017] Based on the total time window t all , the window length t of the moving time window mw , moving step t step , get the segmented signal A i,j (t),1≤j≤k and B i,j (t), where k refers to the number of segments.

[0018] 3. Step (3) The time delay estimation of the earthquake source position change is specifically as follows:

[0019] Obtaining the offset α between segmented signals based on coda wave interferometry i,j , by the offset α i,j and the elapsed time t j Fitting, determining the polynomial equation, performing differential and inverse differential operations on the polynomial equation, and calculating the time delay β caused by the speed change i,j , remove the time delay and obtain the time delay τ caused by the change of the source position i :

[0020] a. Calculate the offset: Calculate the receiving point R based on the coda wave interferometry method i Segmented signal A before and after the source position changes i,j (t) and B i,j(t) The offset α between i,j ;

[0021] b. Delay estimation of speed change: for the offset α i,j and the elapsed time t j Perform polynomial fitting, determine the polynomial order n by the curve change trend, and determine the polynomial coefficient a based on the minimum residual principle n ,a n-1 ,...,a0, to determine the polynomial equation; perform differential and inverse differential operations on the polynomial equation to obtain the time delay β caused by the speed change i,j , β i,j With t j Do a linear fit, the slope of which is the receiving point R i The estimated velocity disturbance Δv i ;

[0022] c. Calculate the time delay of the earthquake source position change: offset α i,j Subtract the delay β i,j , denoted as τ i,j , and then τ i,j Take the average, which is the receiving point R i Signal A i (t) and B i (t) Time delay τ caused by the change of earthquake source position i .

[0023] 4. Step (4) distance difference calculation is specifically as follows:

[0024] Based on the delay τ obtained in step (3) i Multiplying by the propagation velocity v, we can get the distance from the source point S1 to the receiving point R. i Distance difference Δr i .

[0025] 5. Step (5) of earthquake source location estimation is as follows:

[0026] Based on the distance difference Δr obtained in step (4) i , get the unknown coordinates of the source point S2 to the receiving point R i The distance is combined with the coordinate information of the receiving point and hyperbola positioning is used to estimate the coordinate position of the source point S2;

[0027] A. Estimation of the distance between the source position and the receiving point after the change: First calculate the distance from the initial source position S1 (x1, y1) to the receiving point R i =(X i ,Y i ), the distance d of i≥1 i , combined with the distance difference Δr i , get the source point S2 to the receiving point R iDistance D i ;

[0028] B. Hyperbola positioning: take receiving point R i and R j ,j≠i is the center of the circle, D i and D j As the radius, draw a hyperbola and find the focus of the hyperbola, select the focus close to the initial source point S1 as the source point S2 and the receiving point R i and R j An estimated coordinate position S i,j ,i≠j;

[0029] C. Position averaging: any two receiving points estimate a source position S i,j ,i≠j, combined with S i,j and the initial earthquake source point S1, will be very far away from S1 or from the rest of S i,j The coordinate values ​​of are eliminated and the average of the remaining estimated source positions is the coordinate position of source S2.

[0030] 6. Step (6) of earthquake source location tracking is as follows:

[0031] Based on the velocity disturbance Δv i ,1≤i≤n, take the average value to obtain the average velocity disturbance Δv, and then use Δv to compensate the medium propagation velocity v; then use the position of the source S2 located in step 5 as the initial source point, repeat steps (1) to (5) to locate the next source movement position; and so on, to achieve source position tracking.

[0032] The advantages of the present invention are:

[0033] 1. From a theoretical perspective, the present invention combines the characteristics of coda waves, taking advantage of the fact that coda waves are sensitive to medium changes and have strong coda wave signal energy when received at long distances. It analyzes the coda wave phase changes and distinguishes the time delays caused by medium changes and source position changes. Especially when the medium changes are inhomogeneous, the time delay caused by medium changes is removed, and the time delay caused by source position changes is accurately and reliably calculated, achieving high-precision source position tracking.

[0034] 2. From the perspective of engineering application, the present invention has lower requirements for station distribution, is less dependent on the medium velocity structure, and is more practical. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 is a flow chart of the present invention;

[0036] Figure 2 This is a velocity model diagram of a strong scattering medium according to the present invention;

[0037] Figure 3It is a velocity model diagram of a weak scattering medium of the present invention;

[0038] Figure 4 is a diagram of the movement of the earthquake source position of the present invention;

[0039] Figure 5 This is a signal diagram at the strong scattering medium receiving point R1 of the present invention;

[0040] Figure 6 This is a signal diagram at the weak scattering medium receiving point R1 of the present invention;

[0041] Figure 7 This is a graph showing the nonlinear time delay estimation result of the present invention;

[0042] Figure 8 is the error bar graph of the tracking results of the present invention;

[0043] Figure 9 It is the cumulative distribution function CDF graph of the present invention. DETAILED DESCRIPTION

[0044] The present invention will be described in more detail below with reference to the accompanying drawings:

[0045] Combine Figure 1-9 , the present invention specifically includes the following steps:

[0046] The detailed process of step 1 to obtain the time domain signal before and after the change of the earthquake source position is to construct the simulation signal in the present invention. First, the scattering medium velocity model is constructed, as shown in the attached figure. Figure 2 and attached Figure 3 As shown, determine the medium propagation velocity v, set the initial source point position S1 (x1, y1) and the receiving point R i =(X i ,Y i ), i≥1 coordinate position. Then the source is excited and the receiving point R i Receive and record the time domain signal A at the earthquake source point S1 i (t). Then, the source point is changed to S2(x,y), the source is excited, and the receiving point R i Receive and record the time domain signal B with the earthquake source at S2 i (t);

[0047] Step 2: Signal segmentation: Signal segmentation is used to estimate local delay and observe the change trend of local delay with the elapsed time, which is used to distinguish the cause of delay change later. First, determine the total time window t all , ensuring that the coda signal is dominant in the total time window and the signal-to-noise ratio is high, and the receiving point R i Time domain signal A i (t) and B i (t) Reassignment, i.e. Ai (t) = A i (t all ) and B i (t) = B i (t all ) so that it only contains the signal within the total time window. Set the window length t of the moving time window mw , the moving time window refers to a time window with a specified length t mw The time window is included in the total time window and moves according to the set step size t step Slide the moving window along the time axis to achieve signal segmentation. It should be noted that the moving time window needs to contain enough information, that is, the waveform contained must be complex enough. However, if the window length of the moving time window is too large, the calculation speed will be reduced, so it is necessary to set the window length reasonably. At the same time, in order to ensure the continuity of the analysis results, the moving step size needs to be small so that the signals in the moving time window overlap. Assume that the total time window t all =[t1-t w ,t k +t w ], the window length of the moving time window is t mw =2t w , where the center position of the sliding time window is t j ,1≤j≤k, where k is the number of segmented signals. The present invention makes the moving time window contain five cycles of waveform changes, and the moving step is set to 1 / 15 of the window length, that is, t step =|t j -t j-1 |=1 / 15×t mw , obtain segmented signal A i,j (t),1≤j≤k and B i,j (t),1≤j≤k;

[0048] Step 3: The time delay estimation process of the earthquake source position change includes:

[0049] 3.1: Calculate the offset. When the source position or medium changes, the offset between signals will show a long-term trend of change. Therefore, the offset between signals can be used to detect medium changes, source position changes, or the presence of both medium changes and source position changes. Taking the Windowed Cross-Correlation (WCC) method as an example, it is a coda wave interferometry method to calculate the receiving point R i Segmented signal A before and after the change of earthquake source position i,j (t) and B i,j Assume that the time delay variable t s In [-2t w ,2t w ] range, find the receiving point Ri The jth segment signal A i,j (t) and B i,j (t), t∈[t j -t w ,t j +t w ]'s mutual correlation coefficient R(t s ), the formula is as follows, the time delay variable t corresponding to the maximum value of the mutual correlation coefficient s , that is, obtain the segmented signal A i,j (t) and B i,j The offset α of (t) i,j ,1≤j≤k, where k refers to the number of segmented signals;

[0050]

[0051] 3.2: Delay estimation for speed changes. Uniform speed changes result in linear delays between signals, while non-uniform speed changes result in non-linear delays. Plot the offset α calculated in step 3.1. i,j As time t passes j The polynomial order is determined by the trend of the curve change. It is recommended to set it to 2nd order or above to avoid the situation where the delay may show linear changes when the total time window is later, but may actually be more inclined to nonlinear changes. i,j and the elapsed time t j Perform polynomial fitting and determine the appropriate polynomial coefficient a based on the principle of minimum residual n To a0, that is, to determine the polynomial equation y = a n x n +a n-1 x n-1 +...+a1x+a0, where x∈[1,t k ], t k represents the elapsed time t j The last value of y is obtained to obtain the change value of y. The difference operation based on the variable y is to subtract the previous number from the next number, that is, X(m) = y(m+1)-y(m), and m∈[1,t k -1], the X variable is the result of the differential operation. Then do the reverse differential operation on X, Z(n)=X(n)+Z(n-1), n∈[1,t k -1] and set Z(0) = 0, the Z variable represents the result of the reverse difference operation. Finally, the delay β caused by the non-uniform speed change is calculated i,j , β i,j =Z(t j -1),1≤j≤k; i,j With t j Do a linear fit, the slope of which is the receiving point R iThe estimated velocity disturbance Δv i ;

[0052] 3.3: Find the time delay of the change of the earthquake source position. i,j Subtract the delay β i,j , denoted as τ i,j , that is, remove the time delay caused by speed change. Then τ i,j Take the average, that is, get the receiving point R i Signal A i (t) and B i (t) Time delay τ caused by the change of earthquake source position i , τ i =(τ i,1 +...+τ i,k ) / k;

[0053] Step 4 The distance difference calculation process is for each receiving point R i The calculated delay τ i Converted to the receiving point R before and after the source position changes i The distance difference Δr i , Δr i =τ i ×v, v refers to the medium propagation speed;

[0054] The detailed process of the earthquake source location estimation in step 5 is:

[0055] 5.1: Estimation of the distance between the source and the receiving point after the source position changes. Initial source point S1 and receiving point R i The coordinate position of is known, according to the formula Get the initial source point S1 to each receiving point R i The distance d i Then calculate the distance difference Δr according to step 4 i , get the source point S2 to the receiving point R i Distance D i , D i =d i -Δr i ;

[0056] 5.2: Hyperbolic positioning. From step 5.1, we get the source point S2 to the receiving point R i Distance D i , although S2 is obtained to the receiving point R i The distance, but the coordinates of S2 can be R i is the center of the circle, D i Any point on the arc with radius R. i is the center of the circle, D i As the radius; then take the receiving point R j,j≠i is the center of the circle, D j = radius, find the focus through the hyperbola, that is, solve the system of simultaneous equations to find the coordinates S2(x,y). The system of equations is shown below. Theoretically, there are two solutions, that is, there are two coordinates of S2. However, when calculating the relative source position by coda wave interferometry, the distance between the sources must be less than one wavelength, otherwise there will be period jumps. Therefore, the focus close to the initial source point S1 is selected as the source point S2 by the receiving point R i and R j An estimated coordinate position S i,j ,i≠j;

[0057]

[0058] 5.3: Position averaging. Every two receiving points can estimate a source position S i,j ,i≠j, combined with S i,j and the initial earthquake source point S1, will be very far away from S1 or from the rest of S i,j The coordinate values ​​of are eliminated and the average of the remaining estimated source positions is the coordinate position of source S2;

[0059] Step 6: The velocity disturbance Δv obtained at all receiving points i ,1≤i≤n take the average value and obtain the average velocity disturbance Δv=(Δv1+,...,+Δv n ) / n, compensating for the medium propagation velocity v, that is, v = v + Δv. Using the location of source S2 located in step 5 as the initial source point, repeat steps 1 to 5 to locate the next source movement. This process is repeated to achieve source location tracking.

[0060] The following is a simulation example to further illustrate it:

[0061] First, construct the scattering medium velocity model diagram, such as Figure 2 and Figure 3 As shown, Figure 2 is the velocity model of strong scattering medium, Figure 3 The velocity model of the weak scattering medium is used. The coordinates of the initial earthquake source point S1 are (1606m, 1600m), and the coordinates of S2 are (1606m, 1602m). They move in an approximately counterclockwise circular motion. Figure 4As shown, three receiving points R1 (3958m, 3900m), R2 (3958m, 3958m), and R3 (1978m, 3958m) were selected to receive the recorded signals. Subsequently, the SUFDMOD2 module in the Seismic Unix software was used to obtain the simulation signal. The velocity field boundary of this simulation was set to the absorbing boundary, and a Ricker wavelet with a main frequency of 200Hz was used as the source. The propagation time was 5s, the time sampling rate was 4000 points, and the medium propagation velocity v was the velocity after taking the average of the velocity model, for example, Figure 2 The velocity v obtained by the medium-strong scattering medium model is 3322m / s. Figure 5 A1(t) and Figure 6 are the time domain signals at the initial source point S1 received by the receiving point R1 in strong scattering media and weak scattering media, respectively. Figure 6 The direct wave has strong energy and short duration. The signal energy is between 10 -2 level, and Figure 5 The direct wave is not obvious, and the signal energy is 10 -3 The level indicates that the direct wave is not obvious in the strong scattering medium and the overall signal energy is weak, which is similar to the energy of the coda wave in the weak scattering medium. It also reflects that the coda wave signal may dominate in the strong scattering medium, and the traditional positioning method based on the direct wave will not be applicable. In the same way, the time domain signal of the source point S2 received at the receiving point R1 is nonlinearly stretched to construct the simulation signal B1(t) with both non-uniform velocity changes and source position changes, as shown in Figure 2. Figure 5 As shown in B1(t), controllable medium changes can be obtained.

[0062] In step 2, take the signal A1(t) received by R1 as an example, the total time window t all Set it to [3.75s, 4.25s], and reassign A1(t), that is, A1(t)=A1(t all ). For ease of processing, time is expressed in time points, and the conversion relationship is the time in seconds multiplied by the time sampling rate, that is, t all is [15001, 17000], the window length of the moving time window is t mw is 150, the moving step is t step Set to 10, then the number of segmented signals k = 184, and the segmented signal A is obtained. 1,1 ,A 1,2 ,...,A 1,184 , where A 1,1 =A1(15001:15150), A 1,2 =A1(15151:15300), and so on. Similarly, when the parameters are set the same, the segmented signal B will also be obtained from the B1(t) signal. 1,1 ,B 1,2,...,B 1,184 The time domain signals A2(t) and B2(t) as well as A3(t) and B3(t) at the receiving points R2 and R3 can be processed in the same manner to obtain segmented signals.

[0063] Signal A processed by step 3.1 1,1 (t) and B 1,1 (t), get the delay α 1,1 =-2.62, and the same applies to A. 1,2 (t) and B 1,2 (t) Until A 1,184 (t) and B 1,184 (t) Obtain the horizontal offset α at the receiving point R1 1,2 ,α 1,2 ,...,α 1,184 The segmented signals of receiving points R2 and R3 are processed in the same way to obtain the horizontal offset α at different receiving points. 2,1 ,α 2,2 ,... and α 3,1 ,α 3,2 ,.... Draw α 1,1 ,α 1,2 ,...,α 1,184 As time t passes j ,1≤j≤k, where t j The time corresponding to the center position of the moving time window, for example, t1 = (15001 + 15150) / 2 = 15076, as shown in the attached figure. Figure 7 midpoint dash line α 1,j As shown, a second-order polynomial is selected to fit the curve and the polynomial equation y = -1.695e -08 x 2 +1.077e -05 x + 0.7598, where x∈[1,16906]. Perform a differential operation on y, for example, X(1) = y(2) - y(1) = 1.0722e -05 , and so on, we get X(2), X(3), ..., X(16095). Then we get the Z value based on X, for example, Z(1) = X(1) + Z(0) = 1.0722e-05, Z(2) = Z(1) + X(2) = 2.1411e -05 , where Z(0)=0, and so on, until Z(16095) is obtained. Then, Z is taken together with the elapsed time t j The corresponding value is the delay β caused by non-uniform speed change 1,j , β 1,j =Z(t j -1),1≤j≤k, such as Figure 7As shown by the black solid line in the middle, it is found that the trend of the non-uniform velocity change is similar to that of the theory, and the nonlinear time delay can be estimated well, but there is a certain error. 1,j With the elapsed time t j Do a linear fit, and its slope is the velocity disturbance Δv1 = -0.0005312. 1,j and β 1,j Subtract and take the average, which is the time delay caused by the change of the source position estimated by the receiving point R1, that is, τ1=(α 1,1 -β 1,1 )+,...+(α 1,k -β 1,k ) / k. Similarly, estimate the time delays τ2 and τ3 caused by the change in the source position of the receiving points R2 and R3, as well as the velocity disturbances Δv2 and Δv3.

[0064] According to the formula Δr i =τ i ×v to obtain the distance differences Δr1, Δr2, and Δr3. Calculate the distances d1, d2 and d3 from the initial earthquake source point S1 to the receiving points R1, R2 and R3. Then use formula D i =d i -Δr i , estimate the distance from the source point S2 to each receiving point R i The specific values ​​are shown in Table 1. Combined with the position of the receiving point, the simultaneous equations are used to estimate the coordinate position of the source point S2. The position of the hyperbolic positioning of the receiving points R1 and R2 is S 1,2 =(3490.2402,645.1071), similarly, we get S 1,3 =(1605.7274,1602.1878) and S 2,3 =(1605.7791,1601.8645), where S 1,3 and S 2,3 The positioning result of S is close to that of S1, while S 1,2 The result is far from the S1 position, so S is not considered. 1,2 Position, for S 1,3 and S 2,3 The average of the results is taken as the estimated value of S2, that is, S2 = (1605.7533, 1602.0261), which is basically consistent with the theoretical value S2 = (1606, 1602).

[0065] All receiving points R i The estimated velocity disturbance Δv i,1≤i≤3, take the average, and obtain the average velocity disturbance Δv=(Δv1+Δv2+Δv3) / 3, and compensate the medium propagation velocity, that is, v=v+Δv, and use the estimated coordinate information of S2 as the initial source point to estimate the coordinate position of S3, and so on, to achieve source point tracking. The positioning results are plotted as error bars, as shown in Figure 8 As shown in the figure, the length of the end line indicates the size of the error. The up, down, left, and right directions indicate that the Z coordinate estimate is too large, the Z coordinate estimate is too small, the X coordinate estimate is too small, and the X coordinate estimate is too large. Overall, the X coordinate error is small and the Z coordinate error is obvious, but the trajectory of the source position change can also be estimated more accurately. In order to more intuitively describe the positioning error, the cumulative distribution function (CDF) evaluation index is selected to measure it. The cumulative distribution function CDF is as follows: Figure 9 The CDF represents the positioning measurement error on the horizontal axis, and the CDF on the vertical axis represents the percentage of positioning times at a certain positioning accuracy relative to the total number of positioning times. The maximum positioning error of the WCC method, corresponding to a positioning error of 1, is 1.58 m. Overall, the positioning results are not significantly affected, demonstrating that coda wave interferometry can accurately track earthquake source locations.

[0066] Table 1 Estimated parameters

[0067]

Claims

1. A method for tracking earthquake source location based on coda wave interferometry, characterized by: (1) Set up multiple groups of receiving points R i , determine the initial source position and medium propagation velocity v, at each receiving point R i Get the time domain signal A before the source position changes i (t); After the source position changes, at each receiving point R i Obtain the time domain signal B after the source position changes i (t); 1≤i≤n; (2) The time domain signal A i (t) Segmentation is performed to obtain segmented signal A i,j (t); the time domain signal B i (t) is segmented to obtain segmented signal B i,j (t); 1≤j≤k; k is the number of segments; (3) Obtaining segmented signal A based on coda wave interferometry i,j (t) and B i,j (t) The offset α between i,j ; For offset α i,j and the elapsed time t j Perform polynomial fitting, perform differential and inverse differential operations on the fitted polynomial equation, and calculate the time delay β caused by non-uniform speed changes i,j ; Delay β i,j With the elapsed time t j Do a linear fit and use the slope of the linear fit as the receiving point R i The velocity disturbance Δv at i ; The offset α i,j Subtract the delay β i,j Denoted as τ i,j , and then τ i,j Take the average and use it as the receiving point R i Signal A i (t) and B i (t) Time delay τ caused by the change of earthquake source position i ; (4) The delay τ i Multiplying by the medium propagation velocity v, we can get the value of the source position before and after the change to the receiving point R. i The distance difference Δr i ; (5) Receiving point R i The coordinate position of the earthquake source before the position change is known, and the distance difference Δr i , using hyperbola positioning, estimate the coordinates of the earthquake source after the change; (6) The velocity disturbance Δv obtained at all receiving points i Take the average value to obtain the average velocity disturbance Δv, and compensate the medium propagation velocity v, that is, v = v + Δv; use the coordinates of the source position after the change in step (5) as the initial source point, repeat steps (1) to (5) to locate the next source movement position, and realize source position tracking.

2. The method for tracking earthquake source position based on coda wave interferometry according to claim 1, wherein: Step (1) is specifically as follows: Determine the initial source position S1(x1,y1) and the receiving point position R i =(X i ,Y i ) and the medium propagation speed v, at each receiving point R i Get the time domain signal A before the source position changes i (t); After the source position changes to S2(x,y), at each receiving point R i Obtain the time domain signal B after the source position changes i (t).

3. The method for tracking earthquake source position based on coda wave interferometry according to claim 1, wherein: The specific method of step (2) signal segmentation is: Based on the total time window t all , the window length t of the moving time window mw , moving step t step , obtain segmented signal A i,j (t) and B i,j (t).

4. The method for tracking earthquake source position based on coda wave interferometry according to claim 1, wherein: In step (3), the offset α i,j and the elapsed time t j Perform polynomial fitting, determine the polynomial order N by the curve change trend, and determine the polynomial coefficient a based on the minimum residual principle N ,a N-1 ,...,a0, to determine the polynomial equation.

5. The method for tracking earthquake source position based on coda wave interferometry according to claim 2, wherein: Step (5) is specifically as follows: A. Estimation of the distance between the source position and the receiving point after the change: First calculate the distance from the initial source position S1 (x1, y1) to the receiving point R i =(X i ,Y i ), the distance d i , combined with the distance difference Δr i , get the source point S2 to the receiving point R i Distance D i ; B. Hyperbola positioning: take receiving point R i and R j ,j≠i is the center of the circle, D i and D j As the radius, draw a hyperbola and find the focus of the hyperbola, select the focus close to the initial source point S1 as the source point S2 and the receiving point R i and R j An estimated coordinate position S i,j ,i≠j; C. Position averaging: any two receiving points estimate a source position S i,j ,i≠j, combined with S i,j and the initial earthquake source point S1, will be very far away from S1 or from the rest of S i,j The coordinate values ​​of are eliminated and the average of the remaining estimated source positions is the coordinate position of source S2.

Images