Forward modeling method for seismic point source excited atmospheric acoustic gravimetric waves

Abstract

The invention discloses a forward modeling method for atmospheric sound weight waves excited by a seismic point source, and relates to the field of geophysics. Comprising the following steps that a surface harmonic coordinate system is introduced, a frequency domain atmospheric control equation and a stratum elastodynamic equation are expanded in the surface harmonic coordinate system, an atmospheric control equation set of a frequency-wavenumber domain and a PSV equation set and an SH equation set which are decoupled in the stratum are obtained, and a frequency-wavenumber domain wave field in a horizontal layering model is solved through a global matrix method. Under the condition that a source and medium parameters are given, all continuous boundary conditions and contribution of the source are utilized, a linear equation set with wave amplitude as an unknown number in each layer of medium can be obtained, the equation set is solved, and displacement and stress in the stratum and speed and pressure in the atmosphere are obtained. A frequency-wavenumber domain can be converted into a time-space domain wave field by utilizing Hankel transform on wavenumber and fast Fourier transform on frequency.

Description

Technical Field

[0001] The present invention relates to the field of geophysical technologies, and particularly to a forward modeling method for seismic point source excitation of atmospheric acoustic gravity waves. Background Art

[0002] Many observations show that natural disasters such as earthquakes, volcanoes, and tsunamis are accompanied by atmospheric and ionospheric disturbances, indicating the existence of the coupling phenomenon among the lithosphere - atmosphere - ionosphere spheres. Atmospheric acoustic gravity waves are one of the important mechanisms causing the sphere coupling. When an earthquake occurs, seismic waves will cause surface disturbances and then cause atmospheric disturbances. Propagating to the ionosphere will disturb the plasma and cause ionospheric disturbances. Since the atmospheric density decreases exponentially with height, the velocity of atmospheric particles will increase as they propagate upward. Even a very small surface disturbance will cause large disturbances when propagating to the atmosphere and ionosphere.

[0003] In the current forward modeling studies of seismic excitation of acoustic gravity waves, most use the finite difference method, which requires processing a huge amount of data, resulting in slow calculation and low efficiency. In addition, many studies simplify the seismic source and cannot consider the fault mechanism.

[0004] It should be noted that the above content belongs to the technical cognition scope of the inventor and does not necessarily constitute the prior art. Summary of the Invention

[0005] Based on this, it is necessary to provide a forward modeling method for earthquake-induced atmospheric acoustic gravity waves with high efficiency and saving computing resources in view of the problems mentioned in the above background art, and to investigate the influence of different fault mechanisms on acoustic gravity waves.

[0006] The object of the present invention can be achieved through the following technical solutions: A forward modeling method for seismic point source excitation of atmospheric acoustic gravity wave field based on a horizontally stratified atmosphere - stratum model, the forward modeling method comprising:

[0007] S1: Establish a horizontally stratified model;

[0008] The horizontally stratified model includes N layers of atmosphere and L layers of stratum, where the topmost layer and the bottommost layer are uniform semi-infinite spaces. The seismic source is placed on the z-axis at a depth of z = z s with its coordinate being (0, 0, z s ). In the horizontally stratified model, the acoustic gravity wave field in each layer of medium above the surface obeys the atmospheric control equation, and the seismic waves in each layer of medium below the surface obey the elastic dynamics equation;

[0009] S2: Introduce the surface harmonic coordinate system, expand the frequency-domain atmospheric control equation and the stratum elastic dynamics equation in the surface harmonic coordinate, obtain the atmospheric control equation set and the decoupled PSV equation set and SH equation set in the frequency - wavenumber domain, and derive the variables with the amplitude coefficients as unknowns;

[0010] S3: Solve the frequency - wavenumber domain wave field in the horizontally layered model by the global matrix method. Given the source and medium parameters, using all the continuous boundary conditions and the contribution of the source, a linear system of equations with wave amplitude as the unknown in each layer of the medium is obtained. Solving this system of equations can yield the displacement and stress in the formation and the velocity and pressure in the atmosphere;

[0011] S4: The frequency - wavenumber domain results can be transformed into the time - space domain wave field by using the Hankel transform of the wavenumber and the fast Fourier transform of the frequency.

[0012] Furthermore, the introduction of the surface harmonic coordinate system in S2 includes:

[0013] Introduce a set of surface harmonic coordinate vectors:

[0014]

[0015] where e r , e θ and e z represent the unit vectors in the radial, tangential, and vertical directions in the cylindrical coordinate system respectively. Y m (r,θ) = J m (kr)e imθ (m = 0, ±1, ±2,...), and J m (kr) is the Bessel function of the first kind of order m, where k is the horizontal wavenumber. In the frequency domain, any vector in the cylindrical coordinate system is expanded in the following form:

[0016]

[0017] where the vector ξ is V, P1 in the atmosphere, where V is the atmospheric particle velocity and P1 is the atmospheric pressure perturbation, and is also the displacement u and the stress vector τ = τ rz e r +τ θz e θ +τ zz e z in any position on the horizontal plane in the formation. The value of m depends on the azimuthal symmetry of the source. For an explosive source, m = 0; for a seismic moment tensor source, m ≤ 2, are the components of the vector ξ in the cylindrical harmonic coordinate system.

[0018] Furthermore, the variable expression with the amplitude coefficient as the unknown in S2 is specifically operated as follows: Use variable substitution to find the expression of the pressure P in the atmospheric control equations, and then find the expression of the velocity v z . At any point in the atmosphere, there is

[0019]

[0020] Among them

[0021]

[0022] Among them, the subscript -n represents the parameters of the -nth layer in the atmosphere, and k z represents the vertical wave number of the acoustic gravity wave, g is the acceleration due to gravity, and c s is the speed of sound, and are the amplitude coefficients of the downward wave and the upward wave in the -nth layer respectively, and are variables introduced for convenience in solving.

[0023] Furthermore, the specific operations of S3 are as follows:

[0024] For the PSV mode wave, jointly solve the atmospheric wave equation and the PSV mode equation set, consider the continuous boundary conditions and the contribution of the source, and obtain the global matrix one;

[0025] For the SH mode wave, solve the SH mode equation set separately to obtain the global matrix two;

[0026] Solve the global matrix one and the global matrix two, and then obtain the amplitude expressions of the PSV mode wave and the SH mode wave at any point in the formation, and then calculate the frequency - wavenumber domain field at any point in the horizontal layered model;

[0027] The amplitude expression of the PSV mode wave at any point in the formation is:

[0028]

[0029] The SH mode amplitude expression is:

[0030]

[0031] Among them and are the amplitude coefficients of the downward wave and the upward wave of the P wave in the lth layer underground. and are the amplitude coefficients of the downward wave and the upward wave of the SV wave in the lth layer underground. and are the amplitude coefficients of the downward wave and the upward wave of the SH wave in the lth layer underground.

[0032] Furthermore, the global matrix one is:

[0033]

[0034] Among them is the eigenvector of the up-going wave on the Nth above-ground floor, and M -n are the eigenvectors of the down-going and up-going waves on the (n + 1)th and nth above-ground floors respectively. N -1 and N 1 are the eigenvectors of the down-going and up-going waves on the 1st above-ground floor and the 1st underground floor respectively, G l and are the eigenvectors of the down-going and up-going waves on the lth and (l + 1)th underground floors respectively. S V represents the source discontinuity vector in the PSV mode wave;

[0035] The global matrix two is:

[0036]

[0037] where R1 is the eigenvector of the down-going and up-going waves on the 1st underground floor, H l and are the eigenvectors of the down-going and up-going waves on the lth and (l + 1)th underground floors respectively, J L is the eigenvector of the down-going wave on the bottommost Lth underground floor, S H represents the source discontinuity vector in the SH mode wave.

[0038] Furthermore, the specific process of S4 is:

[0039] Using the formula

[0040]

[0041] transform the velocity V and pressure P1 in the atmosphere and the displacement u and stress τ in the formation in the wavenumber domain to the spatial domain, i.e., cylindrical coordinates;

[0042] Then use the inverse Fourier transform formula to transform the atmospheric particle velocity V, atmospheric pressure P1, formation displacement u and formation stress τ in the frequency domain to the time domain, and finally transform from cylindrical coordinates to Cartesian coordinates to display the results;

[0043] From the formula

[0044]

[0045] transform the frequency - spatial domain vector to the time - spatial domain, where the vector ξ is V and P1 in the atmosphere, or the displacement u and stress vector τ = τ rz e r +τ θz e θ +τ zz e z .

[0046] The present invention has the following advantages compared with the prior art:

[0047] 1. The present invention provides a forward modeling method for seismic-induced atmospheric acoustic gravity waves in a horizontally stratified atmosphere-stratosphere model. The present invention introduces a set of surface harmonic coordinate basis vectors, expands the frequency-domain control equations in the surface harmonic coordinates, and decouples the seismic wave equations into PSV mode and SH mode; then, the frequency-wavenumber domain acoustic gravity wave field in the horizontally stratified model is solved by the global matrix method; finally, the Hankel transform and the fast Fourier transform are used to integrate the frequency and wavenumber to obtain the time-space domain wave field. Compared with the forward modeling method based on the simplified source-induced acoustic gravity waves, the calculation results of the present invention are closer to the observed values.

[0048] 2. The present invention obtains the three-dimensional atmospheric acoustic gravity wave field in the horizontally stratified model through a semi-analytical algorithm, and examines the influence of different fault mechanisms on the acoustic gravity waves, with high calculation efficiency, providing a theoretical basis for studying the influence of earthquakes on atmospheric waves and providing a forward operator for inverting seismic information from atmospheric waves. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] The drawings described herein are used to provide a further understanding of the present invention, and constitute a part of the present invention. The schematic embodiments and descriptions thereof of the present invention are used to explain the present invention, and do not constitute an improper limitation of the present invention. In the drawings:

[0050] Figure 1 is a schematic flow chart of a forward modeling method for acoustic gravity waves excited by a seismic point source in a horizontally stratified earth-atmosphere model provided by the present invention;

[0051] Figure 2 is a schematic diagram of the three-dimensional horizontally stratified model of the present invention;

[0052] Figure 3 is a diagram of the atmospheric model parameters used in the embodiment of the present invention;

[0053] Figure 4 is a schematic diagram of the response of seismic waves and acoustic gravity waves induced by an earthquake in a vertical plane in the horizontally stratified model of the present invention;

[0054] Figure 5 For the present invention Figure 4 is a schematic diagram of the velocity and pressure responses at the positions of the four black triangles;

[0055] Figure 6 is a schematic diagram of acoustic gravity waves excited by a seismic source at different depths at the same observation point;

[0056] Figure 7 is a schematic diagram of the influence of different strike angles of a vertical reverse fault on acoustic gravity waves;

[0057] Figure 8Schematic diagram of acoustic gravity wave waveforms excited by different fault mechanisms at the same observation point;

[0058] Figure 9 Schematic diagram of acoustic gravity wave signals excited by the source under different formation velocity structures. Detailed implementation manners

[0059] For a clearer explanation of the overall concept of the present invention, the following will be described in detail by way of examples in conjunction with the accompanying drawings of the specification.

[0060] A forward modeling method for acoustic gravity waves excited by a seismic point source in a horizontally stratified earth-atmosphere model, as Figure 1 shown, the forward modeling method includes:

[0061] S1: Establish a horizontally stratified model;

[0062] The structure inside the earth is inhomogeneous. Generally speaking, the inhomogeneity in the vertical direction is greater than that in the horizontal direction. The actual formation structure often presents a structure pattern similar to a horizontally stratified one. Thus, the earth medium can be regarded as a horizontally stratified medium model composed of multiple parallel medium layers. For the atmosphere, the inhomogeneity in the vertical direction is also greater than that in the horizontal direction. Therefore, a horizontally stratified model is established.

[0063] The horizontally stratified model includes N layers of the atmosphere and L layers of the formation. The topmost layer and the bottommost layer are uniform semi-infinite spaces. The source is placed on the z-axis at a depth of z = z s with its coordinates being (0, 0, z s ). In the horizontally stratified model, the acoustic gravity wave field in each layer of the medium above the ground surface obeys the atmospheric control equation, and the seismic wave in each layer of the medium below the ground surface obeys the elastic dynamics equation.

[0064] S2: In order to quickly solve three-dimensional problems using the method of coordinate transformation, a set of orthogonal and complete coordinate systems, namely the spherical harmonic coordinate system, is introduced. The atmospheric control equation and the formation elastic dynamics equation in the frequency domain are expanded in the spherical harmonic coordinates to obtain the atmospheric control equation set and the decoupled PSV mode equation set and SH mode equation set in the frequency-wavenumber domain, and the variable expressions with the amplitude coefficients as unknowns are derived. This step also transforms the equations in the frequency-space domain to the frequency-wavenumber domain.

[0065] Introducing the spherical harmonic coordinate system includes:

[0066]

[0067] where e r 、e θ and e z respectively represent the unit vectors in the radial, tangential, and vertical directions in the cylindrical coordinate system. Y m(r,θ) = J m (kr)e imθ (m = 0, ±1, ±2,...), and J m (kr) is the Bessel function of the first kind of order m, where k is the horizontal wavenumber. In the frequency domain, any vector in the cylindrical coordinate system is expanded in the following form:

[0068]

[0069] Transform the frequency - space domain vector to the time - space domain, where the vector ξ is V and P1 in the atmosphere, or the displacement u and the stress vector τ = τ rz e r + τ θz e θ + τ zz e z .

[0070] The value of m depends on the azimuthal symmetry of the source. For an explosive source, m = 0; for a seismic moment tensor source, m ≤ 2. are the components of the vector ξ in the cylindrical harmonic coordinate system. The expanded vector includes:

[0071]

[0072] Transform the frequency - wavenumber domain vector to the frequency - space domain.

[0073] The atmospheric governing equations expanded in the spherical harmonic coordinates are

[0074]

[0075] where ρ b and P b are the background atmospheric density and pressure distributions, is a variable introduced for convenience in solving, g is the acceleration due to gravity, and are the components of the velocity spherical harmonic vector respectively, and are the density and pressure perturbation spherical harmonic scalars respectively. Derive the variable expressions with the amplitude coefficients as unknowns. At any point in the atmosphere, there is

[0076]

[0077] where

[0078]

[0079] where the subscript - n represents the parameters of the - nth layer in the atmosphere, k zrepresents the vertical wavenumber of the acoustic gravity wave, g is the acceleration due to gravity, and c s is the speed of sound, and are the amplitude coefficients of the downgoing wave and the upgoing wave in the -nth layer respectively, and are variables introduced for the convenience of solution.

[0080] The decoupled PSV mode equation in the formation is:

[0081]

[0082] where and are the components of the displacement surface harmonic vector, and are the components of the stress surface harmonic vector, and are the elastic moduli, V p and V s are the velocities of the P-wave and S-wave, and ρ s is the density of the solid medium. and are the components of the body force surface harmonic vector.

[0083] The SH mode equation is:

[0084]

[0085] where and are the components of the displacement and stress surface harmonic vectors, is the component of the body force surface harmonic vector.

[0086] S3: Solve the frequency-wavenumber domain wave field in the horizontally layered model by the global matrix method. Given the source and medium parameters, using all the continuous boundary conditions and the contribution of the source, a linear system of equations with wave amplitudes as unknowns in each layer of the medium is obtained. Solving this system of equations can obtain the displacements and stresses in the formation and the velocities and pressures in the atmosphere. The specific operations are as follows:

[0087] For the PSV mode wave, by jointly solving the atmospheric wave equation and the PSV mode equation system, considering the continuous boundary conditions and the contribution of the source, the global matrix one can be obtained:

[0088]

[0089] where is the eigenvector of the upgoing wave in the N layers above the ground, and M -n are the eigenvectors of the down- and upgoing waves in the (n + 1)th and nth layers above the ground respectively.N -1 and N 1 are the downward and upward wave eigenvectors of the first above - ground floor and the first underground floor respectively. G l and are the downward and upward wave eigenvectors of the l - th underground floor and the (l + 1)-th underground floor respectively. S V represents the discontinuity vector of the source in the PSV - mode wave;

[0090] For the SH - mode wave, considering the continuous boundary condition and the contribution of the source, the global matrix two is obtained as:

[0091]

[0092] where R1 is the eigenvector of the downward and upward waves of the first underground floor, H l and are the eigenvectors of the downward and upward waves of the l - th and (l + 1)-th underground floors respectively, J L is the eigenvector of the downward wave of the bottom - most L - th underground floor, S H represents the discontinuity vector of the source in the SH - mode wave.

[0093] Solve the global matrix one and the global matrix two, and then obtain the amplitude expressions of the PSV - mode wave and the SH - mode wave at any point in the formation, and then find the frequency - wavenumber domain wave field at any point in the horizontal stratified model.

[0094] The amplitude expression of the PSV - mode wave at any point in the formation is:

[0095]

[0096] The SH - mode amplitude expression is:

[0097]

[0098] where and are the amplitude coefficients of the downward and upward waves of the P - wave in the l - th underground floor. and are the amplitude coefficients of the downward and upward waves of the SV - wave in the l - th underground floor. and are the amplitude coefficients of the downward and upward waves of the SH - wave in the l - th underground floor.

[0099] Among them, the PSV - mode wave is coupled with the atmospheric acoustic gravity wave, while the SH - mode wave is not coupled with the atmospheric acoustic gravity wave, and the two need to be solved separately.

[0100] S4: Use the inverse Fourier transform and Hankel transform to obtain the results in the time - space domain, and the specific process is as follows:

[0101] Using the formula

[0102]

[0103]

[0104] Transform the velocity V and pressure P1 in the atmosphere and the displacement u and stress τ in the formation in the wavenumber domain to the spatial domain, i.e., cylindrical coordinates;

[0105] Then use the inverse Fourier transform formula to transform the atmospheric particle velocity V, atmospheric pressure P1, formation displacement u, and formation stress τ in the frequency domain to the time domain, and finally transform from cylindrical coordinates to Cartesian coordinates to display the results;

[0106] From the formula

[0107]

[0108] Transform the frequency - spatial domain vector to the time - spatial domain, where the vector ξ is V and P1 in the atmosphere, or the displacement u and stress vector τ = τ rz e r +τ θz e θ +τ zz e z .

[0109] Table 1 and Figure 3 respectively give the formation and atmospheric model parameters used in the present invention. Considering a 42 - layer atmosphere and a 5 - layer formation model, simulating seismic waves and atmospheric acoustic gravity waves excited by an underground point source, the adopted fault mechanism is 0, 90, 90 (default unit is degree), the focal depth is 12 km, and the seismic moment component is M0 = 6.3×10 18 Nm, equivalent to an Mw 6.5 earthquake.

[0110] Figure 4 Gives Figure 2 The wave - field snapshot of the vertical plane shown in b. It can be seen that the earthquake excites P - waves, S - waves, and Rayleigh waves underground. After transmission through the ground, P - wave head waves (acoustic gravity waves caused by P - waves), R - waves (head waves of acoustic gravity waves caused by Rayleigh waves), and A - waves (direct waves generated at the epicenter) are generated in the atmosphere.

[0111] Among them, the P - wave head wave is relatively weak and not very obvious. It can be seen that due to the atmospheric stratification structure, there are obvious reflection and transmission phenomena, especially at an altitude of 100 km.

[0112] Figure 5 Gives Figure 4The velocity and pressure responses at the locations of the four black triangles shown are for an epicentral distance of 500 km and heights of 1 m, 100 km, 200 km, and 300 km respectively. It can be seen that the arrival time of the R-wave gradually delays with the increase in height. After the R-wave propagates upward, it reflects from a height of 100 km to the Earth's surface and then reflects and rises to a height of 300 km (R1). Compared with the R-wave, the R1-wave is delayed by approximately 800 s at each height. The arrival time of the A-wave at a height of 300 km is earlier than that at lower heights, which can also be seen from Figure 4 this as well.

[0113] Figure 6 The acoustic gravity waves excited by earthquake sources at different depths at the same observation point are given. It can be seen that with the increase in the depth of the earthquake source, the change in the R-wave is not significant, but the amplitude of the A-wave decreases significantly.

[0114] Figure 7 The influence of different strike angles of a vertical reverse fault on the acoustic gravity waves is given. It can be seen that both the R-wave and the A-wave have energy related to the strike angle, that is, they are weak in the strike direction.

[0115] Figure 8 The acoustic gravity waves excited by different fault mechanisms at the same observation point are given. The arrival times of the R-wave and the A-wave remain unchanged, but their waveforms and amplitudes vary greatly with the strike angle, dip angle, and slip angle. The above results show that the propagation of atmospheric acoustic gravity waves depends on the source mechanism and has potential application prospects in earthquake source location and earthquake source mechanism inversion.

[0116] Figure 9 The acoustic gravity wave signals excited by earthquake sources under different formation velocity structures are given. The blue line is for the half-space model, and the red line is for the 5-layer horizontal stratified model. It can be seen that the waveforms of the P-AGW (acoustic gravity wave caused by P-wave) and the R-wave in the 5-layer model are different from those in the half-space model, but the A-wave hardly changes. This indicates that the formation structure has little influence on the A-wave, but has a greater influence on the acoustic gravity waves caused by P-waves and Rayleigh waves.

[0117] Table 1 Formation velocity structure model used in the present invention

[0118]

[0119] Each embodiment in this specification is described in a progressive manner. For the same or similar parts among the embodiments, reference can be made to each other. Each embodiment focuses on the differences from other embodiments. In particular, for the system embodiment, since it is basically similar to the method embodiment, the description is relatively simple, and for the relevant parts, reference can be made to the partial description of the method embodiment.

[0120] The above are only embodiments of the present invention and are not intended to limit the present invention. For those skilled in the art, various modifications and variations can be made to the present invention. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included within the scope of the claims of the present invention.

Claims

1. A forward modeling method for atmospheric acoustic gravity waves excited by earthquake point sources, characterized in that: The forward modeling method comprises: S1: Establish a horizontal hierarchical model; The horizontal layered model contains N layers of atmosphere and L layers of strata, where the top and bottom layers are uniform semi-infinite spaces and the earthquake source is placed at a depth of z = z s On the z-axis, its coordinates are (0, 0, z s ), in the horizontal layered model, the acoustic-gravity wave field in each layer of medium above the surface obeys the atmospheric control equation, and the seismic wave in each layer of medium below the surface obeys the elastic-dynamic equation; S2: Introduce the surface harmonic coordinate system, expand the frequency domain atmospheric control equations and the formation elastic dynamics equations in the surface harmonic coordinates, obtain the atmospheric control equations in the frequency-wavenumber domain and the decoupled PSV equations and SH equations in the formation, and derive the variables with the amplitude coefficient as the unknown number; S3: The frequency-wavenumber domain wave field in the horizontal layered model is solved by the global matrix method. Given the source and medium parameters, all continuous boundary conditions and source contributions are used to obtain a linear equation system with wave amplitude as the unknown in each layer of the medium. Solving this equation system can obtain the displacement and stress in the formation and the velocity and pressure in the atmosphere. S4: The frequency-wavenumber domain results can be converted into the time-space domain wave field using the Hankel transform of wavenumber and the fast Fourier transform of frequency.

2. The forward modeling method of atmospheric acoustic gravity waves excited by earthquake point sources according to claim 1, characterized in that: The surface harmonic coordinate system is introduced in S2, including: Introduce a set of surface harmonic coordinate vectors: where e r 、e θ and e z Respectively represent the unit vectors in the radial, tangential and vertical directions in the cylindrical coordinate system, Y m (r,θ)=J m (kr)e imθ (m=0,±1,±2,...), and J m (kr) is the first kind m-order Bessel function, where k is the horizontal wave number. In the frequency domain, any vector in the cylindrical coordinate system Expands into the following form: The vector ξ is V in the atmosphere, P1, where V is the atmospheric velocity, P1 is the atmospheric pressure disturbance, and is also the displacement u and stress vector τ=τ on the horizontal plane at any position in the formation. rz e r +τ θz e θ +τ zz e z , the value of m depends on the annular symmetry of the source. For explosive sources, m = 0; for seismic moment tensor sources, m ≤ 2, is the component of the vector ξ in the cylindrical harmonic coordinate system.

3. The forward modeling method of atmospheric acoustic gravity waves excited by earthquake point sources according to claim 1, characterized in that: The variable expression with the amplitude coefficient as the unknown in S2 is specifically operated as follows: using variable substitution to find the expression of the pressure P in the atmospheric control equations, and then find the velocity v z The expression for any point in the atmosphere is in The subscript -n represents the parameters of the nth layer in the atmosphere, k z represents the vertical wave number of the acoustic-gravity wave, g is the gravitational acceleration, c is s is the speed of sound, and are the amplitude coefficients of the downgoing and upgoing waves in the -nth layer, and It is a variable introduced for the convenience of solving the problem.

4. The forward modeling method of acoustic and gravity waves excited by earthquake point sources based on a horizontally layered atmosphere-stratum model according to claim 1 is characterized in that: The specific operation of S3 is as follows: For PSV mode waves, the atmospheric wave equations and the PSV mode equations are jointly solved, taking into account the continuous boundary conditions and the contribution of the sources, and the global matrix 1 is obtained; For SH mode waves, the SH mode equations are solved separately to obtain the global matrix 2; Solve the global matrix 1 and global matrix 2, and then obtain the amplitude expressions of PSV mode waves and SH mode waves at any point in the formation, and then calculate the frequency-wavenumber domain field at any point in the horizontal layered model; The amplitude expression of the PSV mode wave at any point in the formation is: The expression of SH mode amplitude is: in and is the amplitude coefficient of the downgoing and upgoing P waves in the underground layer l; and is the amplitude coefficient of the downgoing and upgoing SV waves in the underground layer l; and It is the amplitude coefficient of the downgoing and upgoing SH waves in the underground layer 1.

5. The forward modeling method of atmospheric acoustic gravity waves excited by earthquake point sources according to claim 4, characterized in that: The global matrix 1 is: in is the characteristic vector of the upward wave in N layers above the ground, and M -n are the characteristic vectors of the downward and upward waves of the n+1th and nth layers on the ground, respectively. N -1 and N 1 are the down-going and up-going wave eigenvectors of the first floor above the ground and the first floor underground, respectively. G l and are the downward and upward wave eigenvectors of the underground layer l and layer l+1 respectively, S V represents the discontinuity vector of the source in the PSV mode wave; The global matrix 2 is: in R 1 is the characteristic vector of the downward and upward waves of the first underground layer, H l and are the characteristic vectors of the downward and upward waves of the lth and l+1th underground layers, respectively, J L is the characteristic vector of the downward wave of the lowest underground layer L, S H Represents the discontinuity vector of the source in the SH mode wave.

6. The forward modeling method of atmospheric acoustic gravity waves excited by earthquake point sources according to claim 1, characterized in that: The specific process of S4 is as follows: Using the formula Transform the velocity V and pressure P1 in the atmosphere, displacement u and stress τ in the formation in the wave number domain to the spatial domain, i.e., cylindrical coordinates; Then use the inverse Fourier transform formula to transform the atmospheric gas point velocity V, atmospheric pressure P1, formation displacement u and formation stress τ in the frequency domain into the time domain, and finally convert from cylindrical coordinates to rectangular coordinates to display the results; By formula Transform the frequency-space domain vector to the time-space domain, where the vector ξ is V and P1 in the atmosphere, or the displacement u and stress vector τ = τ on the horizontal plane at any position in the formation rz e r +τ θz e θ +τ zz e z .

Images