KR101614138B1 - Method for 3-dimensional moment tensor inversion by dividing modeling region - Google Patents

Abstract

The present invention relates to a method for moment tensor inversion used for identifying the direction of a fault when an earthquake occurs. According to the present invention, in a method for moment tensor inversion used for copying a waveform in a mid- or large-scale earthquake and the epicenter occurring in case of hydraulic fracturing in the field of seismology studies, to settle the issues of the methods of moment tensor inversion of the conventional technologies where there was a limit, wherein it was necessary to suppose a first-dimensional model and obtain an interpretational answer to reduce the calculation time for inversion of moment tensor waveform and a great calculation time and economic costs were required when a three-dimensional model is applied to overcome those issues, without needing to perform an inversion for the entire three-dimensional model, calculation is conducted by dividing the entire area of the three-dimensional model between a source and geophone, and a perfectly matched layer (PML) is applied to a time region three-dimensional viscoelastic finite-difference method to remarkably reduce the calculation time and costs.

Description

A three-dimensional moment tensor inversion method using a modeling domain decomposition method

The present invention relates to a moment tensor inversion technique used for grasping the direction of a fault in an earthquake. More particularly, the present invention relates to a moment tensor inversion technique used in earthquake research, The conventional moment tensor inversion method used for the simulation of the moment tensor waveform inversion has a disadvantage in that it takes much time to inverse the moment tensor waveform. In order to solve the problem that the calculation time is reduced, the three-dimensional moment tensor inversion ≪ / RTI >

In order to overcome the limitations of the prior art, in which the analytical solution has to be obtained assuming a one-dimensional model in order to reduce the calculation time, the present invention adopts a moment inversion technique using waveform simulation for a three- In order to overcome the disadvantages of the conventional moment tensor inversion techniques which require a lot of computation time and economic cost to simulate the model waveform, the entire area is decomposed into a model between the source and the receiver, and time domain three-dimensional viscoelastic The present invention relates to a three-dimensional moment tensor inversion method using modeling domain decomposition configured to significantly reduce computation time and cost by applying a PML (perfectly matched layer) to a finite difference method.

Generally, a wave inversion technique is performed by constructing the same environment as a survey environment using a computer simulation, measuring the error with the actually measured data, and measuring the underground medium velocity It is a technique that obtains the propagation velocity of elastic wave by updating the distribution, and it has been attracting public attention in recent years due to the advantage of completely disregarding the subjective judgment of the constructor.

In addition, the elastic-wave imaging technology is the most important part of the exploration technology. It uses real seismic data to collect diffraction and scattered data at one point, corrects the reflection surface to the correct position, In recent years, the reverse imaging technique using wave equation has been actively studied worldwide due to the development of computer performance.

That is, as an example of the conventional technique for the seismic wave exploration method using the waveform inversion as described above, for example, a method and an apparatus for estimating the underground medium structure through waveform inversion stepwise as disclosed in Korean Patent Registration No. 10-1262990, And a method and apparatus for imaging an underground structure using waveform inversion, which are disclosed in Korean Patent Publication No. 10-1092668, there have been proposed methods for imaging an underground structure using waveform inversion.

As another example of the conventional art for the seismic wave detection method using the waveform inversion as described above, Korean Patent Registration No. 10-1182838, for example, discloses a method for calibrating a frequency domain inverse time structure through transmission source estimation, &Quot;, and Korean Patent Publication No. 10-2013-0054236 discloses "total wave length inversion using time-varying filters ".

As described above, conventionally, various methods using waveform inversion have been proposed for estimating and imaging the underground structure through seismic wave exploration. However, the contents of the above-mentioned prior arts include such waveform inversion methods as applied to seismic exploration There is no description of the contents of the technique.

That is, in the earthquake research in the related art, as a method used for grasping the direction of an oscillation mechanism, that is, a fault layer in the event of an earthquake, for example, a FOCMEC (see Reference 1) method using a vertical component polarity and an amplitude ratio , And the moment tensor inversion technique (see Reference 2) using waveform matching is the most widely used.

Recently, in the field of seismology, moment tensor inversion is frequently used to simulate waveforms for an earthquake equal to or more than a medium scale, and in recent years, moment tensor inversion techniques have been widely used for the origin occurring during hydraulic fracturing.

More specifically, in the past, a method of using an analytic solution assuming a one-dimensional model due to the time-consuming time of the moment tensor waveform inversion has been widely used. In recent years, however, Moment inversion techniques using waveform simulations for 3D model are applied.

However, the conventional moment inversion technique for a three-dimensional model has a disadvantage in that calculation time is long due to a large number of calculation areas since both the source for the acoustic wave simulation and the area for the water depth for the 3D model are required to be calculated.

In addition, the conventional moment tensor inversion method for a three-dimensional model has a problem that it is necessary to perform a complicated calculation process repeatedly in order to simulate the waveform of the three-dimensional model. That is, when the moment tensor inversion is performed, There is a disadvantage in that a lot of calculation time and economic cost are required because the elastic wave simulation must be performed.

Therefore, in order to solve the problem of the conventional moment tensor inversion techniques, it is required to reduce the processing time required for the simulation of the elastic wave for the 3D model. For this purpose, for example, It is preferable to decompose the three-dimensional model and calculate only the model part between the transmission source and the water depth, so that the calculation time is shortened. However, there is not provided a device or a method satisfying all of such requirements yet In fact.

[references]

1. Snoke, JA, JW Munsey, AC Teague, and GA Bollinger, 1984, "A Program for Focal Mechanism Determination of Polarity and SV-P Amplitude Ratio Data", Earthquake Notes, Vol. 3. pp. 15

2. Minson, S. E. and Dreger, D. S., 2008, "Stable inversions for complete moment tensors", GJI, Vol. 174, pp. 585-592.

3. JH Chang, SJ Min, 2012, "Application of ADE-PML Boundary Conditions for Numerical Simulation of 3-Dimensional SEM Seismic Wave using Velocity-Stress Variation Method", Geophysical and Geophysical Exploration, Vol 15, No.2 pp.57- 65.

4. Komatitsch, D., and Martin, R., 2007, "An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation", Geophysics, Vol. 72, pp. SM155-SM167.

5. Robertsson, J., Blanch, J., and Symes, W., 1994, "Viscoelastic finite-difference modeling", Geophysics, Vol. 59, pp. 1444-1456

6. Carcione, J. M., 2001, Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic and Porous Media, Elsevier, New York.

7. Ben-Menahem and Singh, 1981, Ben-Menaehe and Singh S J 1981 Seismic waves and sources, Springer-Verlag, New York.

8. Carcione, J. M., Kosloff, D. and Kosloff, R., 1988, "Wave propagation simulation in a linear viscoelastic medium", Geophys. J. R. Astr. Soc. vol 95, 597-611.

9. Diaz, J., Ezziani, A., and Le Goff, N., 2011, "Version 2.0 Gar6more3D", http://web.univ-pau.fr/~jdiaz1/gar63DCecill.html.

10. Kang, TS and Shin, JS, 2006, "Surface-wave tomography from ambient seismic noise of accelerograph networks in southern korea", GEOPHYSICAL RESEARCH LETTERS, VOL. 33, L17303, doi: 10.1029 / 2006GL027044

[Prior Art Literature]

1. Korean Patent Registration No. 10-1262990 (Feb.

2. Korean Patent Registration No. 10-1092668 (Dec. 2011)

3. Korean Patent Registration No. 10-1182838 (2012.09.07.)

4. Korean Patent Publication No. 10-2013-0054236 (May 31, 2014).

SUMMARY OF THE INVENTION The present invention has been made to solve the problems of the prior art as described above, and it is therefore an object of the present invention to provide a method and apparatus for grasping the direction of a fault in an earthquake, In order to solve the problems of the conventional moment tensor inversion techniques, which take a long time in the moment tensor waveform inversion in the moment tensor inversion technique, there is no need to perform inversion on the entire three-dimensional model Dimensional moment tensor inversion using a modeling domain decomposition configured to decompose the entire area of the 3D model into a model between the source and the receiver, thereby reducing the computation time significantly.

It is another object of the present invention to provide a method for estimating the moment tensor inversion by applying a moment inversion technique to a three-dimensional model in order to overcome the problem, In order to solve the problems of the conventional moment tensor inversion techniques which take a long calculation time and increase the economic cost, there is no need to perform the inverse calculation for the entire 3D model, Dimensional moment tensor using modeling domain decomposition, which is constructed to decompose into a model between a water surface and a water depth, and to apply a PML (perfectly matched layer) to a time domain three-dimensional viscoelastic finite difference method, And to provide a method of inversion.

In order to achieve the above object, according to the present invention, in the moment tensor inversion used for grasping the fault direction and waveform simulation for an earthquake when an earthquake occurs, Dimensional moment tensor inversion method using modeling domain decomposition to solve the problem of the conventional moment tensor inversion methods which requires calculation time and cost because of necessity of calculation, A modeling step of constructing a 3D velocity structure model; A modeling domain decomposition step of decomposing the 3D velocity structure model constructed in the modeling step into a model for a region between a source and a receiver; Mapping the coordinates of the moment tensor with respect to the model between the transmission source and the water depth decomposed in the modeling region decomposition step to coordinates on the X-axis; A simulation step of simulating a model between the transmission source and the water depth by using a moment tensor mapped to coordinates on the X-axis in the mapping step; A reverse mapping step of returning the result calculated in the simulation step back to the original coordinates; And a moment tensor inverse calculation step of performing moment tensor inverse calculation on a model between the transmission source and the receiver obtained in the inverse mapping step, by executing the process on the computer or dedicated hardware, Wherein the model region is calculated by performing only the inverse calculation on the model between the source and the receiver without performing the inverse calculation on the modeling domain decomposition, Dimensional moment tensor inversion method using the three-dimensional moment tensor inversion method.

Here, the mapping step may be configured to perform a process of rotating the axis of the moment tensor by an angle between the transmission source and the water depth by using the rotation operator (R) shown in the following equation.

Further, in the simulation step, when the moment tensor (m) is represented by the following equation,

A new moment tensor (m ') for the coordinate system mapped in the mapping step is calculated using the following equation, and a process of performing a simulation using the new moment tensor (m') is performed do.

(Where R ^t is a transpose matrix of R)

In addition, the reverse mapping step may be configured to perform a process of rotating the components calculated in the simulation step by an angle between the transmission source and the water depth calculator using the rotation operator (R ') shown in the following equation do.

Further, the inverse mapping step is characterized in that the processing for obtaining the original composite wave S is performed by rotating the component of the calculated composite wave S 'using the following equation.

Also, the moment tensor inversion step may be performed using a three-dimensional viscoelastic finite difference method to perform three-dimensional moment tensor inversion, and to reduce an artificial reflection wave generated at a finite boundary by a calculation area reduction, (Auxiliary Differential Equations) -PML (Perfectly Mached Layer).

Here, the inversion method may include a process of performing three-dimensional moment tensor inversion using an SEM (Spectral Element Method) or a Finite Element Method (FEM) instead of the finite difference method in the moment tensor inversion step To be performed.

According to the present invention, there is also provided a recording medium on which a program configured to cause a computer to execute a process according to the three-dimensional moment tensor inversion method using the modeling domain decomposition described above.

Furthermore, according to the present invention, a moment tensor inversion system is provided that is configured to perform moment tensor inversion using the three-dimensional moment tensor inversion method using the modeling domain decomposition described above.

As described above, according to the present invention, in the moment tensor inversion technique used for grasping the direction of a fault in an earthquake occurrence, waveform simulation for a medium or large scale earthquake, and exact simulation occurring during hydraulic fracturing, The 3D model is decomposed into a model between the source and the receiver so as to calculate the entire area of the 3D model without performing the inversion of the entire model, By providing the moment tensor inversion method, it is possible to solve the problem of the conventional moment tensor inversion methods which take a long time in inversion of the moment tensor waveform.

According to the present invention, as described above, the entire area of the three-dimensional model is decomposed into a model between the source and the receiver, and the time domain three-dimensional viscoelasticity The present invention provides a three-dimensional moment tensor inversion method using modeling domain decomposition configured to significantly reduce computation time and cost by applying a perfectly matched layer (PML) to a finite difference method. In order to overcome this problem, there is a problem in that the calculation time is long and the economic cost is increased by applying the moment inversion technique to the three-dimensional model in order to overcome this problem. Can solve the problems of moment tensor inversion techniques.

1 is a view showing a velocity-stress staggered finite difference lattice.
2 is a diagram schematically showing a configuration of a moment tensor transmission source in a finite difference method.
FIG. 3 is a diagram for explaining rotation of domain decomposition and moment tensor components for a three-dimensional model; FIG.
FIG. 4 is a diagram for explaining a concept of mapping to the X-axis in order to facilitate the calculation of the decomposed model area.
FIG. 5 is a view schematically showing a velocity structure model applied to a verification test for verification of a three-dimensional finite difference algorithm applied to the present invention, and a layout of a transmission source and a receiver.
6 is a diagram showing the results of a verification test performed using the velocity structure model as shown in Fig.
FIG. 7 is a diagram showing a result of verifying the safety in calculating a long-term numerical model using the method according to the embodiment of the present invention.
FIG. 8 is a graph showing the relationship between the composite wave for the entire three-dimensional model and the transmittance of the transmitted wave in the three-dimensional model shown in FIG. 3, when a monopod strand 45.degree., Oblique 45.degree. And a composite wave obtained by decomposition into a model.
FIG. 9 is a graph showing the relationship between the composite wave for the entire three-dimensional model and the transmittance of the transmitted wave in the three-dimensional model shown in FIG. 3, when 45 degrees, 45 degrees, And a composite wave obtained by decomposition into a model.
10 is a diagram illustrating a three-dimensional velocity structure model applied to verify a moment tensor inversion method according to an embodiment of the present invention.
FIG. 11 is a graph showing a result of synthesized waves calculated by SEM using the three-dimensional velocity structure model shown in FIG. 9 and a result of synthesized waves calculated by model decomposition regions according to an embodiment of the present invention, respectively.
FIG. 12 is a diagram showing a result of moment inversion using a conventional three-dimensional SEM.
13 is a diagram showing a result of moment inversion using the decomposition area three-dimensional FDM according to an embodiment of the present invention.
FIG. 14 is a flowchart schematically showing the overall configuration of a three-dimensional moment tensor inversion method using modeling domain decomposition according to an embodiment of the present invention.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS Hereinafter, a three-dimensional moment tensor inversion method using modeling region decomposition according to the present invention will be described with reference to the accompanying drawings.

Hereinafter, it is to be noted that the following description is only an embodiment for carrying out the present invention, and the present invention is not limited to the contents of the embodiments described below.

In the following description of the embodiments of the present invention, parts that are the same as or similar to those of the prior art, or which can be easily understood and practiced by a person skilled in the art, It is important to bear in mind that we omit.

That is, as described later, in the moment tensor inversion technique used for grasping the tomographic direction in case of an earthquake, waveform simulation for a medium or larger scale earthquake, and exact simulation occurring during hydraulic crushing, In order to solve the problem of the conventional moment tensor inversion techniques which have a disadvantage that it takes a long time to invert the tensor waveform, it is necessary to perform the inverse calculation on the whole three-dimensional model, The present invention relates to a three-dimensional moment tensor inversion method using modeling domain decomposition configured to be decomposed into a model so as to significantly reduce calculation time.

In order to reduce computation time in the moment tensor inversion, the present invention imposes a one-dimensional model on the assumption that an analytical solution has to be obtained. To overcome this problem, a moment inversion technique for a three- In order to solve the problem of the conventional moment tensor inversion techniques, which takes a long calculation time and increases the economic cost, there is no need to carry out the inverse calculation for the entire 3D model, The modeling decomposition is constructed by decomposing and calculating the model between the source and the receiver, and applying the PML (perfectly matched layer) to the time domain three-dimensional viscoelastic finite difference method to reduce the computation time and cost remarkably. Dimensional moment tensor inversion method.

Next, a specific embodiment of the three-dimensional moment tensor inversion method using the modeling domain decomposition according to the present invention constructed as described above will be described with reference to the drawings.

Before describing the details of the three-dimensional moment tensor inversion method using the modeling domain decomposition according to the present invention, a viscoelastic wave equation for applying moment tensor inversion and an algorithm for obtaining the solution will be described. First, The equation can be expressed by the following equation (1).

[Equation 1]

Where ρ is the density, v is the velocity component, σ is the stress, f is the force, 1 is the x-axis, 2 is the y-axis, and 3 is the z-axis.

In addition, the particle velocity-stress relation can be expressed as shown in [Equation 2] (see Reference 5).

&Quot; (2) "

Where lambda is the Lame constant, and mu is the shear modulus.

In addition, the memory variable can be expressed by the following equation (3).

&Quot; (3) "

및

는 각각 이하에 제시되는

및

과 같고,

는

=

과 같다. here,

And

Lt; RTI ID = 0.0 >

And

Lt; / RTI >

The

=

Respectively.

In order to adjust to a constant Q value in a specific frequency domain, the relaxation frequency should be firstly parameterized. In this case, the fitting should be performed on a logarithmic frequency domain scale , Which provides an optimal fit for the viscoelastic model.

In other words, the P wave and the S wave in the elastic medium are characterized by attenuation and dispersion due to the difference in Q values (refer to Reference 6), and the unrelaxed shear wave modulus for the S wave is expressed by [Mathematical As shown in Fig.

&Quot; (4) "

In addition, the unrelaxed modulus is expressed by Equation (5) below.

&Quot; (5) "

Where v _s is the S wave velocity and ρ is the density.

In addition, the Q value for the S wave is expressed by the following equation (6).

및

은 전단응력(sheer stress)과 전단변형(shear strain)에 대한 완화(relaxation) 시간을 나타내는 것으로 할 때, 완화 주파수는 이하의 [수학식 7]과 같이 나타낼 수 있다.
Furthermore, L _s is the maximum number for the mitigation mechanism,

And

Represents a shear stress and a relaxation time for shear strain, the relaxation frequency can be expressed by the following equation (7). &Quot; (7) "

&Quot; (7) "

Here, the weighting factor is expressed by the following equation (8).

&Quot; (8) "

Therefore, the relaxation frequency for defining the Q value for the p-wave can be expressed by the following equation (9).

&Quot; (9) "

The weighting factor is expressed by Equation (10) below.

&Quot; (10) "

과

은 팽창응력(dilatation stress) 및 팽창변형(dilatation strain)에 대한 완화(relaxation) 시간을 나타낸다.
Where L _p is the number of P-wave damping mechanisms

and

Represents the relaxation time for the dilatation stress and the dilatation strain.

Also, for ease of calculation, L = Ls = L _p and _{ω l = ω sl = ω pl} , l = 1, ..., if as L, L increases if the number of Q _s (ω) and Q _p ([omega]) is improved, which represents the usable frequency range of the source.

Therefore, if Q is expressed with respect to the P wave, the following Equation (11) is obtained (refer to Reference 7).

&Quot; (11) "

Where v _p is the relaxed P wave velocity, and v _s is the relaxed S wave velocity.

Next, the free surface boundary condition will be described. The expression for the free interface of the three-dimensional viscoelastic medium is expressed by the following equation (12).

&Quot; (12) "

Here, using the equation of stress, it can be expressed by the following equation (13).

&Quot; (13) "

Also, using the imaging method of stress, the surface boundary is set to antisymmetry.

The PML boundary condition is applied to remove the artificial reflection wave. In the conventional PML, it is necessary to separate the vertical component and the horizontal component of the attenuation region, but the CFS-PML (Complex Frequency Shifted PML) Auxiliary Differential Equations (ADE) using -PML has the advantage that it can reduce the storage capacity by using it without the need for component separation like classical PML and it does not need to modify the wave equation itself.

That is, the PML is intended to attenuate the wave according to space and time by taking a complex domain in a spatial variable, and the derivative with respect to the direction can be expressed by the following equation (14).

&Quot; (14) "

이고, d_x는 댐핑(damping) 축 단면, α_x는 감쇠상수이며, k ≥ 1는 ω는 주파수를 나타낸다. here,

D _x is the damping axis section, α _x is the attenuation constant, and k ≥ 1 represents the frequency.

In addition, if k _x = 1 and? _X = 0, this becomes the same as the conventional PML as shown in the following Equation (15).

&Quot; (15) "

Here, the expression (15) can be rewritten as the following expression (16).

&Quot; (16) "

At this time,? _X can be rewritten as shown in the following Equation (17).

&Quot; (17) "

Further, by expanding the equation (17) again and converting it into the time domain, it becomes equal to the following equation (18).

&Quot; (18) "

Further, Equation (18) is referred to as Auxiliary Differential Equation (ADE) for PML application, and it can be expressed as a following equation (19).

&Quot; (19) "

Here, the variable of PML is _{d x (0) = - (} N POWER + 1) * V p (max) * log 10 (R coef) / (2 * L PML), α x = α max * (1 - ( x / L _PML ) ² ), and α _max = π * f _m were used.

At this time, N _POWER is 2, f _m is the main frequency of the transmission source, and R _coef = 0.001 reflection coefficient.

First, referring to FIG. 1, FIG. 1 is a diagram illustrating a velocity-stress staggered finite difference grating.

That is, the present inventors used a staggered grid finite difference as shown in FIG. 1 for the above-described elastic wave simulation.

Although the present invention has been described with reference to the case where the finite difference method is used for the elastic wave simulation as described above in the embodiment of the present invention described below, the present invention is not necessarily limited to this case, It should be noted that the present invention can be applied to an SEM (Spectral Element Method) method or a FEM (Finite Element Method) in addition to the above-described finite difference method.

That is, referring to FIG. 2, FIG. 2 is a diagram schematically illustrating a configuration of a moment tensor transmission source in a finite difference method.

More specifically, in the finite difference method, the moment tensor transmission source generates nine components (M _{xx (t)} , M _{xy (t)} , M _{xz (t)} , M _{yx (t), M yy (} t), M yz (t), M zx (t), M zy (t), M zz (t)) is to be given.

3 is a view for explaining the 3D model region decomposition and the rotation of moment tensor components. FIG. 3A is a three-dimensional model constructed in the same manner as the conventional method, and FIG. 3B is a three- 5 is a view showing an exploded view of a model region between a fifth receiver and a transmission source.

More specifically, as shown in FIG. 3A, when there is a source and a receiver for the elastic wave simulation for a three-dimensional model, the conventional method requires a long calculation time because these areas must be calculated.

Particularly, in performing moment tensor inversion, it is necessary to repeatedly perform the elastic wave simulation from 7th to 9th times. Therefore, in order to reduce the elastic wave simulation time, the present inventors have proposed a three- The new inversion method is implemented by considering that the computation time can be reduced considerably by decomposing only the model part between the source and the receiver without the need to calculate the whole.

As shown in FIG. 3B, even if the 3D model is decomposed into a model between the source and the receiver, there is a problem that it is very difficult to calculate the 3D model directly due to the problem of the coordinate system. Accordingly, And the coordinate transformation is performed on the x-axis of the decomposed model so that the calculation can be performed very easily.

That is, referring to FIG. 4, FIG. 4 is a diagram for explaining a concept of mapping to the X-axis to facilitate calculation of the decomposed model area.

At this time, not only the model region but also the moment transmission source must be rotated. To do so, the model and moment tensor may be rotated by -θ angle.

The rotation operator at this time is as shown in the following equation (20).

&Quot; (20) "

In addition, the moment tensor can be expressed as shown in the following equation (20).

&Quot; (21) "

Here, the moment tensor m 'in the new coordinate system must be calculated, and the newly calculated moment tensor can be expressed by the following equation (22).

&Quot; (22) "

Here, R ^t is a transpose matrix of R, that is, according to the embodiment of the present invention, the simulation may be performed using the newly calculated moment tensor.

In addition, it is necessary to rotate the calculated component through the simulation. The rotation conversion equation used here is as shown in the following equation (23).

&Quot; (23) "

In addition, when the composite wave calculated using the model decomposition is referred to as S 'component and the desired composite wave is referred to as S, this can be obtained by performing rotation as shown in the following equation (24).

&Quot; (24) "

In the above-described embodiments of the present invention, the present invention has been described by taking the case of using the model domain decomposition and the component rotation of the moment tensor as the finite difference method. However, the present invention is not limited to this case, The concept of model domain decomposition and moment tensor component rotation proposed in the present invention can be applied to, for example, a spectral element method (SEM) method and a finite element method (FEM) as it is, It should be noted that this is applicable to

Next, the actual performance of the moment tensor inversion method according to the embodiment of the present invention constructed as described above will be described by simulation.

That is, in order to verify whether the algorithm of the three-dimensional finite difference method according to the embodiment of the present invention as described above is actually good, the inventors of the present invention performed a calculation algorithm verification by the Lamb's Problem.

More specifically, referring to FIG. 5, FIG. 5 is a diagram schematically illustrating a velocity structure model applied to a verification test for verification of a three-dimensional finite difference algorithm applied to the present invention, and a transmission source and a receiver arrangement of the velocity structure.

Referring to FIG. 6, FIG. 6 is a view illustrating a result of a verification test performed using the velocity structure model shown in FIG. 5, wherein the interpretation of the cumulative data for the Lamb's problem The result obtained by Gar6More3D (refer to Reference 8) as the final solution is compared with the SEM data and the results calculated by the method according to the embodiment of the present invention, respectively.

That is, as shown in FIG. 5, the present inventors arranged a receiver at a distance of 500 m and a distance of 200 m from the transmission source, applied a waveform of a main frequency of 8 Hz, The result of the spectral element method (SEM) according to the conventional method and the analytical solution by GarMore3D 2.0 (refer to Reference 9) are simultaneously shown for the validity comparison of the obtained results As a result of the comparison, it can be confirmed that the respective results are in good agreement as shown in FIG.

Referring to FIG. 7, FIG. 7 is a diagram illustrating a result of verifying safety in calculating a long-term numerical model using a method according to an embodiment of the present invention.

That is, in order to verify the instability at the artificial interface occurring in the calculation of the long-term numerical model using the method according to the embodiment of the present invention, the present inventors have found that the Lamb's problem (PML boundary condition) Respectively.

At this time, the PML coefficient is applied to the parameter of the PML as d _x (0) = - (N _POWER + 1) * V _p (max) * log ₁₀ (R _coef ) / (2 * L _PML ) Step was used.

As a result of the above-described stability verification, as shown in FIG. 7, it can be confirmed that the method according to the embodiment of the present invention is not divergent even during long-time simulation, and is stable.

Here, in the present invention, the computation using the three-dimensional finite difference method takes a long calculation time, so that the calculation is performed in parallel through the MPI. Also, by applying the PML as the boundary condition, the artificial finite boundary reflection wave can be efficiently reduced Respectively.

That is, in general, when the PML is applied, it tends to diverge in long time simulation, or when the transmission source is close to the PML boundary, a sudden amplitude may increase and disappear at the PML interface. However, , This problem can be solved by selecting an appropriate PML coefficient.

Next, the verification result of the algorithm using the model decomposition between the transmission source and the receiver will be described with reference to FIGS. 8 and 9. FIG.

8 and 9, FIGS. 8 and 9 are graphs showing the relationship between the three-dimensional model shown in FIG. 3 and the total three-dimensional model when a monotonic 45.degree. Slope 45.degree. And the synthesized wave in the case of decomposition into a model between the transmission source and the receiver in accordance with the embodiment of the present invention.

More specifically, for verification of a three-dimensional finite difference method using model decomposition between a source and a receiver according to an embodiment of the present invention as described above, , 45 ° of slope 45 ° and 45 ° of moving direction. The main frequency of the transmission source was a 100 Hz Riker waveform. The receiver was placed at a position of 45m from the transmission source. At each 45 ° position, And the positions of the respective water depths are indicated by numerals 1 to 8 in FIG.

In addition, the composite waves of the three-dimensional model as a whole and the composite waves of the result of decomposition into the model between the transmission source receivers according to the embodiment of the present invention are shown in FIG. 8 according to the conventional method. As a result, , It can be confirmed that the composite waves of the two cases agree well.

In other words, as described above, in the conventional art, there are many calculation areas for performing three-dimensional modent tensor inversion and six model calculations are required in one iteration calculation. Therefore, in the present invention, In order to shorten the computation time, a three-dimensional velocity model, a three-dimensional velocity model, and a space between the transmission source and the water depth calculator, as shown in FIG. 3B, The model area of the model is limited.

FIGS. 8 and 9 show the result of calculation based on the decomposition region with limited models and the result based on the original three-dimensional model. In this case, the transmission source used for the waveform simulation calculation A strike of 45 °, a dip of 45 °, and a rake of 45 ° were used for the moment tensor. The Ricker transmission waveform with a main frequency of 100 Hz was used. The model area is calculated by mapping in the x-axis direction.

Here, the moment tensor should be rotated, and each calculated component should be rotated to calculate the components of the original data direction. As shown in FIGS. 8 and 9, the results show that the two results are in good agreement have.

Therefore, from the above, the method according to the present invention is configured to limit the 3D model region in order to reduce the three-dimensional waveform simulation calculation time for the inversion of the moments using the 3D finite difference modeling method, It can be calculated very efficiently in the case of a three-dimensional velocity structure with few spawning sites, and is particularly useful for inversion of seismic moment using low frequency region.

Next, the verification with the result of the moment tensor inversion will be described.

In the three-dimensional moment tensor algorithm, the variable vector for the moment tensor and the earthquake location information can be expressed by the following equation (25).

&Quot; (25) "

In addition, the change amount with respect to the variables of the moment tensor variable and the seismic source position can be expressed by the following equation (26).

&Quot; (26) "

At this time, the observation data D obtained at each observing station is expressed by the following equation (27).

&Quot; (27) "

In addition, when the synthetic wave data obtained at the respective stations are S as shown in the following Equation (28), the synthetic wave data and the observation data difference r can be expressed by the following Equation (29).

&Quot; (28) "

&Quot; (29) "

Therefore, when the moment tensor inversion is such that the change amount of the synthetic wave for each variable is equal to the following equation (30), the relationship between the variable change amount and the data change amount can be expressed by the following equation (31).

 &Quot; (30) "

 &Quot; (31) "

Since inequality (31) is not a square matrix, it is impossible to carry out inverse calculation immediately, so that inverse calculation can be performed by transforming as in the following expression (32).

(32)

Next, the inversion method according to the present invention as described above is applied to practical examples, and the results compared with the conventional method will be described.

10 and 11, FIG. 10 is a diagram illustrating a three-dimensional velocity structure model applied to verify a moment tensor inversion method according to an embodiment of the present invention, FIG. 11 is a three- And comparing the result of the synthetic wave calculated by SEM with the model and the result of synthesized wave calculated by the model decomposition region according to the embodiment of the present invention.

As shown in FIGS. 10 and 11, the inventors of the present invention have found that the results obtained by inversion of moments using the decomposition domain finite difference three-dimensional simulation according to the embodiment of the present invention as described above, And the results were compared with those obtained by the above.

That is, the present inventors performed verification by comparing the inversion results using the existing SEM three-dimensional simulations and the inversion using the three-dimensional finite difference model applied in the present invention for the Odaja earthquake of January 20, 2007.

More specifically, for the moment inversion comparison of the Mt. Odae earthquake on January 20, 2007, using the three-dimensional velocity structure model of the Korean peninsula as shown in Fig. 10 for the calculation of moment tensor (see Reference 9) 11 are compared with the composite wave results calculated by the SEM and the composite wave results calculated by the model decomposition region according to the embodiment of the present invention, respectively.

In this case, 30 grids were used in the y direction for model decomposition, and one grid interval was 500 m.

As a result, as shown in FIG. 11, the results of comparing the green function synthetic waves of 0.04 to 0.1 Hz at the CHC station position with respect to the three-dimensional velocity structure shown in FIG. 10 are almost identical .

12 and 13, FIG. 12 is a view showing a result of moment inversion using a conventional three-dimensional SEM, and FIG. 13 is a diagram showing a result of moment inversion using the decomposition region three-dimensional FDM according to an embodiment of the present invention. Fig.

12 and 13, it can be seen that there is no significant difference between the two methods. As shown in FIG. 13, the case of using the decomposition area according to the embodiment of the present invention is the same as that of the conventional SEM method I can see that it shows slightly better results.

Therefore, the three-dimensional moment tensor inversion method using the modeling domain decomposition according to the embodiment of the present invention can be implemented as described above.

That is, referring to FIG. 14, FIG. 14 is a flowchart schematically showing the overall configuration of a three-dimensional moment tensor inversion method using modeling region decomposition according to an embodiment of the present invention.

More specifically, as shown in FIG. 14, a three-dimensional moment tensor inversion method using modeling region decomposition according to an embodiment of the present invention includes a 3D velocity structure model for an area to be analyzed, (S141) for constructing a three-dimensional velocity structure model, a modeling domain decomposition step (S142) for decomposing the 3D velocity structure model constructed in the step into a model for a region between a source and a receiver, A mapping step (S143) of transforming the coordinate of the moment tensor with respect to the model between the transmission source and the receiver to the coordinates on the x-axis decomposed in the mapping step, and a simulation step A reverse mapping step (S145) of returning the result calculated in the step back to the original coordinates, and a step And a moment tensor inversion step (S146) for performing moment tensor inversion with respect to the model between the first and second models.

Here, the modeling step S141 and the modeling area decomposition step S142 may be performed as described above with reference to FIGS. 1 to 3, and the mapping step S143, S144) and the inverse mapping step S145 may be performed as described above with reference to FIGS. 3 and 4 and [Equation 20] to [Equation 24].

In addition, the moment tensor inversion step S146 may be performed using a three-dimensional viscoelastic finite difference method as described above with reference to Equations 1 to 19 and Equations 25 to 32. [ To apply a PML (Perfectly Mached Layer) boundary condition to the three-dimensional moment tensor inversion.

Furthermore, the three-dimensional moment tensor inversion method using the modeling domain decomposition according to the embodiment of the present invention as described above can be applied not only to the finite difference method but also to the SEM (Spectral Element Method) or the Finite Element Method (FEM) It should be noted that it is applicable.

Further, the present invention performs a series of processes including processing steps as shown in Fig. 14 by using the three-dimensional moment tensor inversion method using the modeling domain decomposition according to the embodiment of the present invention configured as described above Or may be implemented in the form of a program configured to execute the above steps on a computer or dedicated hardware so that the computation time compared to the conventional moment tensor inversion method The moment tensor inversion apparatus and inversion system that can be shortened can be implemented easily and inexpensively.

As described above, according to the present invention, a three-dimensional point acoustic wave simulation is required at the time of three-dimensional moment inversion. However, for example, in the case of a Korean peninsula earthquake, In order to solve the problem that a lot of calculation time is required to calculate, the calculation time can be drastically reduced by dividing the entire model into a model between the source and the destination and having the minimum calculation area.

In addition, according to the present invention, ADE-PML can be applied to reduce the artificial reflection wave generated at the finite boundary that can be generated by reducing the model area as described above, In order to solve the problem that the moment tensor and the component are different from each other by decomposing the model into a model between the source and the receiver, the axis of the moment sensor is rotated by the angle between the source and the receiver, The computed component is rotated by the angle between the source and the receiver, and the decomposed model area can be decomposed and computed in multiple CPUs or cluster machines using multiple MPIs.

That is, according to the present invention, the entire model area is decomposed into a small area between the transmission source and the water depth calculator without calculating the entire area in the conventional three-dimensional moment inversion, and ADE (Auxiliary Differential Equations) -PML Since the frequency band to be calculated at the time of inversion of the moments is low frequency, it is possible to reduce the time required for the calculation of the moment inversion, There is almost no difference between the results of the three-dimensional simulation and those methods.

Therefore, the three-dimensional moment tensor inversion method using the modeling domain decomposition according to the present invention can be implemented as described above.

Further, by implementing the three-dimensional moment tensor inversion method using the modeling domain decomposition according to the present invention as described above, according to the present invention, when an earthquake occurs, grasping the tomographic direction, waveform simulation for medium- In the moment tensor inversion technique used to simulate the origin occurring at a given time, the whole area of the 3D model is decomposed into a model between the source and the receiver, without performing the inversion of the entire 3D model A moment tensor inversion method using a modeling domain decomposition configured to remarkably reduce a calculation time by using a moment tensor inversion method of the prior art which has a disadvantage that it takes a long time to inverse the moment tensor waveform, Can be solved.

In addition, according to the present invention, as described above, the entire area of the 3D model is decomposed into a model between the source and the receiver without performing inversion on the 3D model as a whole, and at the same time, The present invention provides a three-dimensional moment tensor inversion method using modeling domain decomposition configured to significantly reduce computation time and cost by applying a perfectly matched layer (PML) to a finite difference method. In order to overcome this problem, there is a problem in that the calculation time is long and the economic cost is increased by applying the moment inversion technique to the three-dimensional model in order to overcome this problem. Can solve the problems of moment tensor inversion techniques.

While the present invention has been particularly shown and described with reference to exemplary embodiments thereof, it is to be understood that the invention is not limited to the disclosed exemplary embodiments. It will be understood by those skilled in the art that various changes, modifications, combinations and substitutions may be made without departing from the scope of the present invention. .

Claims (9)

In case of moment tensor inversion, which is used for understanding the fault direction in case of an earthquake and simulating a waveform for an earthquake, it is necessary to simulate the entire waveform in the three-dimensional model, so that calculation time is long and cost is required A three-dimensional moment tensor inversion method using modeling domain decomposition to solve the problems of the conventional moment tensor inversion methods,
A modeling step of constructing a 3D velocity structure model for an area to be analyzed;
A modeling domain decomposition step of decomposing the 3D velocity structure model constructed in the modeling step into a model for a region between a source and a receiver;
Mapping the coordinates of the moment tensor with respect to the model between the transmission source and the water depth decomposed in the modeling region decomposition step to coordinates on the X-axis;
A simulation step of simulating a model between the transmission source and the water depth by using a moment tensor mapped to coordinates on the X-axis in the mapping step;
A reverse mapping step of returning the result calculated in the simulation step back to the original coordinates; And
And a moment tensor inverse calculation step of performing a moment tensor inverse calculation on the model between the transmission source and the receiver obtained in the inverse mapping step,
It is possible to reduce the calculation time and cost by decreasing the model area to be calculated by performing only the inversion of the model between the source and the receiver without performing the inverse calculation for the whole three-dimensional velocity structure model Dimensional moment tensor inversion method using modeling domain decomposition.
The method according to claim 1,
Wherein the mapping step comprises:
And a process of rotating the axis of the moment tensor by an angle between the transmission source and the water depth by using the rotation operator (R) shown in the following equation is performed: < EMI ID = .

3. The method of claim 2,
Wherein the simulating step comprises:
When the moment tensor (m) is expressed by the following equation,

A new moment tensor (m ') for the coordinate system mapped in the mapping step is calculated using the following equation, and a process of performing a simulation using the new moment tensor (m') is performed Three - dimensional moment tensor inversion method using modeling domain decomposition.

(Where R ^t is a transpose matrix of R)
The method of claim 3,
Wherein the inverse mapping step comprises:
Wherein processing for rotating components calculated in the simulation step by an angle between the transmission source and the water depth calculator is performed using a rotation operator (R ') shown in the following equation: < EMI ID = Moment tensor inversion method.

5. The method of claim 4,
Wherein the inverse mapping step comprises:
A three-dimensional moment tensor inversion method using modeling domain decomposition, which is characterized in that processing of obtaining the original composite wave S is performed by rotating the component of the calculated composite wave S 'using the following equation .

6. The method of claim 5,
Wherein the moment tensor inversion step comprises:
Dimensional moment tensor inversion is performed using a three-dimensional viscoelastic finite difference method,
To reduce the artificial reflection wave generated at the finite boundary by the reduction of the calculation domain, a modeling feature that is configured to apply a boundary condition according to ADE (Auxiliary Differential Equations) -PML (Perfectly Mached Layer) to the three-dimensional viscoelastic finite difference method Three - Dimensional Moment Tensor Inversion Method Using Domain Decomposition.
The method according to claim 6,
The inversion method includes:
In the moment tensor inversion step, a process of performing a three-dimensional moment tensor inversion process is performed using an SEM (Spectral Element Method) or a Finite Element Method (FEM) instead of the finite difference method Three - dimensional moment tensor inversion method using modeling domain decomposition.
A computer-readable recording medium on which is recorded a program for causing a computer to execute a process according to the three-dimensional moment tensor inversion method using the modeling domain decomposition described in any one of claims 1 to 7.
The moment tensor inversion system is configured to perform moment tensor inversion using the three-dimensional moment tensor inversion method using the modeling domain decomposition described in any one of claims 1 to 7.

Images