CN111797552A - Numerical data simulation method for undulating sea surface seismic wave field based on sea wave spectrum - Google Patents

Numerical data simulation method for undulating sea surface seismic wave field based on sea wave spectrum Download PDF

Info

Publication number
CN111797552A
CN111797552A CN202010535639.5A CN202010535639A CN111797552A CN 111797552 A CN111797552 A CN 111797552A CN 202010535639 A CN202010535639 A CN 202010535639A CN 111797552 A CN111797552 A CN 111797552A
Authority
CN
China
Prior art keywords
wave
sea surface
spectrum
formula
equation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010535639.5A
Other languages
Chinese (zh)
Other versions
CN111797552B (en
Inventor
张进
刘怀山
邢磊
王尧
申鹏
王林飞
尹燕欣
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Ocean University of China
Original Assignee
Ocean University of China
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Ocean University of China filed Critical Ocean University of China
Priority to CN202010535639.5A priority Critical patent/CN111797552B/en
Publication of CN111797552A publication Critical patent/CN111797552A/en
Application granted granted Critical
Publication of CN111797552B publication Critical patent/CN111797552B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Abstract

The invention belongs to the technical field of numerical simulation, and discloses a numerical data simulation method of a fluctuating sea surface seismic wave field based on a sea wave spectrum, which obtains a first-order speed-stress acoustic wave equation set expression through derivation and order reduction processing; respectively carrying out finite difference calculation on each acoustic wave equation of each grid point and each half grid point in time and space according to the calculation rule of the staggered grids; forward modeling of marine data; obtaining a space domain model of a surge background fluctuating sea surface under the P-M spectrum condition; and applying the obtained fluctuating sea surface numerical value to finite difference numerical simulation, and replacing the upper boundary numerical value of the model calculation area with the sea surface numerical value to obtain the simulation data under the surging and fluctuating sea surface. The numerical simulation method provided by the invention has the advantages that the equation is simple to use, the efficiency is high, and the precision can be improved by a certain means; the three-dimensional seismic data acquisition construction cost is low, the sinking depths of the seismic source and the detector are moderate, and the phenomenon of false in-phase axis can be effectively avoided.

Description

Numerical data simulation method for undulating sea surface seismic wave field based on sea wave spectrum
Technical Field
The invention belongs to the technical field of numerical simulation, and particularly relates to a numerical data simulation method of an undulating sea surface seismic wave field based on a sea wave spectrum.
Background
Currently, the closest prior art: in offshore oil and gas exploration, calm sea surfaces are rare due to the influence of factors such as wind, gravity and tide. The sea surface always shows irregular undulations and such undulations are time and space varying. The movement forms of sea water are various, wherein sea waves are the most common and common movement forms of sea water, and sea waves comprise storms and swell waves. During ghost studies, the sea surface is generally assumed to be horizontal, but is actually undulating. Sea surface undulations have a large influence on the delay time of the ghost, so the horizontal sea surface assumption that is commonly used affects the effect of ghost squashing. How to describe the undulating sea surface is an important problem, so the chapter introduces a wave spectrum method which is commonly used in ocean science research, researches the undulating sea surface by using the wave spectrum and models the forward modeling of the undulating sea surface. In physical oceanography, ocean waves are usually described and simulated by ocean wave spectrums, and the ocean waves simulated to a certain degree are considered to be relatively practical.
Various numerical simulation methods are developed, and commonly used methods include wave equation, ray tracing and the like. Among these methods, the ray tracing method is accurate when facing an ideal model, but the treatment effect on the more complex stratum is not ideal and the theoretical requirement is higher. In general, it is practical to simulate with elastic waves, but in indoor simulation, it is more complicated to simulate with elastic waves, and in actual seismic exploration, it is usually the longitudinal wave information that is most utilized.
In the process of acquiring seismic data by the marine streamer, in order to avoid sea surface wind, surge and other various sea surface environmental noises, the seismic source and the geophone are usually sunk to a certain position below the sea surface, so as to obtain data with high signal-to-noise ratio. However, under such construction conditions, the detectors receive not only primary reflected waves but also ghost waves, i.e., ghost waves. Generally, the depth of the seismic source and the detector is shallow in consideration of construction cost, which causes a ghost to appear behind the primary reflection, interfere with the primary wave and change the form of the primary wave, and even generates a false in-phase axis on the seismic record. In addition, due to the frequency trap characteristic of the ghost, the frequency band of the seismic data is narrowed, the resolution of the data is reduced, and the work band of interpretation, inversion and the like of the subsequent seismic data is difficult.
In summary, the problems of the prior art are as follows:
(1) in the existing various numerical simulation methods, the ray tracing method is accurate when facing an ideal model, but the treatment effect on a more complex stratum is not ideal and the theoretical requirement is high.
(2) It is practical to use elastic wave simulation, but in indoor simulation, it is more complicated to use elastic wave simulation, and in actual seismic exploration, it is usually the most utilized information of longitudinal waves.
(3) The existing seismic data acquisition construction cost is high, the sinking depths of a seismic source and a detector are shallow, and a false in-phase axis can be generated on seismic records.
(4) Due to the frequency trap characteristic of the ghost, the frequency band of the seismic data is narrowed, the resolution of the data is reduced, and the work band of interpretation, inversion and the like of the follow-up seismic data is difficult.
The technical problems are solved by the following difficulties:
(1) and (5) derivation of an acoustic wave equation.
(2) Seismic records combining the acoustic wave equation with the wave spectrum are being evolved.
The significance of solving the technical problems is as follows:
if the exact influence of the undulating sea surface and the sea waves on the seismic record can be obtained through numerical simulation calculation, the formation mechanisms of various noises in marine seismic exploration can be explored, and the noises can be suppressed, so that the signal-to-noise ratio of the seismic record is improved.
Disclosure of Invention
Aiming at the problems in the prior art, the invention provides a numerical data simulation method of an undulating sea surface seismic wave field based on a sea wave spectrum.
The invention is realized in such a way, the numerical data simulation method of the wave field of the undulating sea surface based on the wave spectrum explores the influence of the undulating sea surface with the waves on the seismic record by combining the wave spectrum with the sound wave equation, and further acquires the noise forming mechanism in the marine seismic exploration, thereby achieving the purpose of improving the signal-to-noise ratio by suppressing the noise. The method comprises the following steps:
step one, a second-order acoustic wave equation is obtained through derivation, and the second-order acoustic wave equation is subjected to order reduction processing in a mode of integrating time t at two sides of the equation at the same time, so that an expression of a first-order speed-stress acoustic wave equation set is obtained.
And step two, performing finite difference calculation on the acoustic wave equations of each grid point and each half grid point on time and space respectively according to the calculation rule of the staggered grids on the first-order acoustic wave equation obtained in the step one.
Step three, forward modeling of marine data: the PML layer above the model calculation region is not provided with an absorption attenuation coefficient, and the difference calculation of the upper boundary is not performed.
Step four, modeling the surge sea surface based on the wave spectrum: and obtaining a space domain model of the surge background undulating sea surface under the P-M spectrum condition by utilizing the P-M spectrum and combining with the inverse Fourier transform.
And fifthly, applying the obtained numerical value of the undulating sea surface to finite difference numerical simulation, and replacing the numerical value of the upper boundary of the calculation area of the model in the step four with the numerical value of the undulating sea surface to obtain the simulation data under the surging undulating sea surface.
Further, in the first step, the elastic wave equation in the two-dimensional homogeneous medium is derived according to a stress equation, a strain equation and a relationship between stress and strain, and can be expressed by the following formula:
Figure BDA0002536953790000031
equation (1-1) is an elastic wave equation of a two-dimensional homogeneous medium model, expressed in vector form as:
Figure BDA0002536953790000032
ideally, the shear stress τ is absent in the fluid, and has a value of 0, i.e.:
σxy＝σyz＝σzx＝0 (1-3)
the normal stress at this time is:
σxx＝σyy＝σzz＝-p (1-4)
in the above formula, p represents the average pressure of the fluid medium.
The following relationship exists in an ideal fluid:
σxx＝kθ＝λθ＝-p (1-5)
thus, in an ideal fluid, the elastic wave equation can be simplified as follows:
Figure BDA0002536953790000041
in the formula (1-6), where ρ represents density, divergence is taken simultaneously on both sides of the formula (1-6) and the order of differentiation is exchanged, and it can be obtained:
Figure BDA0002536953790000042
after simplification, the following can be obtained:
Figure BDA0002536953790000043
according to the relational expression
Figure BDA0002536953790000044
The formula (1-8) can be rewritten as
Figure BDA0002536953790000045
Wherein the content of the first and second substances,
Figure BDA0002536953790000046
when this is brought into the above formula (1-9), the acoustic wave equation required by the present invention is obtained, which is expressed as follows:
Figure BDA0002536953790000047
in the formula (1-10), V represents the acoustic velocity, ρ represents the density of the formation medium, and p represents the acoustic pressure.
In the actual forward modeling process, a seismic source function is added to the acoustic wave equation of the above formula to realize the true forward modeling of the acoustic wave equation. It is generally assumed that the density of the formation is constant, so that the acoustic equations (1-10) can be expressed as follows:
Figure BDA0002536953790000048
before solving the high-order difference equation of the acoustic wave equation, the formula (1-11) needs to be reduced, and the two sides of the equation are integrated for the time t at the same time to obtain an expression of a first-order velocity-stress acoustic wave equation set, which is as follows:
Figure BDA0002536953790000049
in the formula (1-12), u represents a stress component, vx、vzThe velocity components in the x and z directions are represented, respectively, with ρ representing the density of the formation and v representing the velocity.
Further, in step two, the time second-order precision staggered grid finite difference includes:
assuming a velocity component v in the x-directionxIs derivable at an m-th order time, then is paired at some time t
Figure BDA0002536953790000051
And
Figure BDA0002536953790000052
taylor-series expansion can be carried out to obtain the following formulas (1-13) and (1-14):
Figure BDA0002536953790000053
Figure BDA0002536953790000054
where Δ t is the time grid step, O (Δ t)m+1) Is an infinitesimal quantity of order m +1 with respect to Δ t.
Subtracting the equations (1-13) and (1-14) can obtain:
Figure BDA0002536953790000055
the invention adopts the time second order difference precision, and the formula (1-15) can be simplified as follows:
Figure BDA0002536953790000056
for the same reason, for the velocity component v in the Z directionzBy performing the calculation, we can obtain:
Figure BDA0002536953790000057
the stress u is calculated at the point of the time half-grid, i.e. at
Figure BDA0002536953790000058
The Taylor-series expansion is carried out on t and t + delta t, and the same approximation is carried out, so that:
Figure BDA0002536953790000059
the formulae (1-16), the formulae (1-17) and the formulae (1-18) are each vx、vzAnd u is a second order difference formula over time.
The spatial finite difference format of order 2N comprises:
assuming that any function f (x) has a partial derivative of order 2N +1, respectively
Figure BDA00025369537900000510
And
Figure BDA00025369537900000511
performing Taylor-series expansion:
Figure BDA00025369537900000512
Figure BDA00025369537900000513
Figure BDA0002536953790000061
Figure BDA0002536953790000062
Figure BDA00025369537900000612
Figure BDA0002536953790000063
Figure BDA0002536953790000064
where Δ x represents a spatial grid step length, and the difference between two formulas in the formula (2-21) can be obtained:
Figure BDA0002536953790000065
first partial derivative of equation (2-22)
Figure BDA0002536953790000066
Let its coefficient be 1 and the coefficients of the remaining terms be zero, in which case equation (1-22) can be written as follows:
Figure BDA0002536953790000067
at this time, the coefficient terms of the equations (1-23) can be written in the form of the following equation set
Figure BDA0002536953790000068
Solving equations (1-24) can be calculated using the Cramer's Law in linear algebra: when N is equal to 1, the compound is,
Figure BDA0002536953790000069
when the N is equal to 2, the N is not more than 2,
Figure BDA00025369537900000610
and so on.
In the seismic forward modeling process, in order to ensure the calculation accuracy and control the computation amount, the difference accuracy of the time 2 order, the space 10 order or the time 2 order, the space 12 order is generally used. The first six even-order difference coefficients, as shown in table 1:
TABLE 1 spatial difference coefficient tables of different orders
Figure BDA00025369537900000611
Figure BDA0002536953790000071
Dispersing partial derivatives in the formula (1-12) by using the difference method to obtain finite difference formulas of 2 orders of time and 2 Nth orders of space, wherein the finite difference formulas are shown in the formula (1-25):
Figure BDA0002536953790000072
wherein, Δ X and Δ Z respectively represent space grid step lengths in the X and Z directions, and i and j respectively represent corresponding space grid points; Δ t represents a time grid step, k represents a corresponding time grid point; u, V, W show stress u and velocity component vx、vzIn discrete form.
Further, in the third step, when the variable depth cable is simulated, if the detector is not on a grid point, the wave field value at the detector is calculated by utilizing the wave field interpolation of the two grid points above and below the detector.
The interpolation mode adopted by the inclined cable simulation in the invention is as follows: and calculating weights according to the distances between the ith channel detector (ith column) and two adjacent grid points (i, j) and (i, j +1), and performing interpolation according to the weights to obtain the wave field value at the detector.
Further, in step four, the mathematical formula of the P-M spectrum is:
Figure BDA0002536953790000073
where α is 0.0081, β is 0.74, both of which are constants, g is the gravitational acceleration, ω is the circular frequency, and U is the wind speed above the sea surface, which is an expression of the wave spectrum in the frequency domain. In water, the following relationship exists between frequency and wavenumber:
ω2(k)＝gk (1-37)
the expression of the P-M spectrum in the wavenumber domain is:
Figure BDA0002536953790000074
after obtaining the formula (1-38), multiplying the formula by a group of Gaussian distributed random complex numbers, and performing inverse Fourier transform to obtain a spatial domain model of the swell background undulating sea surface under the P-M spectrum condition.
And after the fluctuating sea surface numerical value is obtained, applying the numerical value to finite difference numerical simulation, and replacing the numerical value of the upper boundary of the model calculation area with the series numerical value to obtain the simulation data under the surging and fluctuating sea surface. When numerical simulation is performed, the size of the grid points cannot be decimal, and the wave spectrum obtained by the method can be decimal, so that rounding processing needs to be performed on the sea surface fluctuation numerical value.
Further, the finite difference conditions include source conditions, stability conditions, dispersion, boundary conditions.
1) Seismic source condition analysis
The seismic source wavelet of the analog seismic data selects a 30Hz zero-phase Rake wavelet, and the expression is as follows:
Figure BDA0002536953790000087
in the formula (1-26), fmIs the dominant frequency of the wavelet, t0Is the delay time of the seismic source.
2) Analysis of stability conditions
The stability condition of the interleaved mesh finite difference numerical solution is as follows:
Figure BDA0002536953790000081
when the formula (1-27) is satisfied, the numerical simulation result is accurate.
3) Frequency dispersion analysis
On the premise of keeping the difference order, the size of the grid can be reduced, and the wavelet dominant frequency suppression numerical value dispersion can be reduced.
4) Analysis of boundary conditions
The invention adopts a complete matching layer (PML) method, seismic waves can not be reflected when propagating to a boundary, the energy of the seismic waves can be absorbed and attenuated in a PML layer, and the absorption and attenuation function of the seismic waves is as follows:
Figure BDA0002536953790000082
in the formula (1-28), R is a theoretical reflection coefficient, vpIs the longitudinal wave velocity and is the thickness of the PML layer.
The difference equation added to PML is:
Figure BDA0002536953790000083
Figure BDA0002536953790000084
Figure BDA0002536953790000085
Figure BDA0002536953790000086
Figure BDA0002536953790000091
in summary, the advantages and positive effects of the invention are:
the numerical data simulation method of the undulating sea surface seismic wave field based on the sea wave spectrum has the advantages of simple equation use and higher efficiency, and the precision can be improved by a certain means. The three-dimensional seismic data acquisition construction cost is low, the sinking depths of the seismic source and the detector are moderate, and the phenomenon of generating false in-phase axes on seismic records can be effectively avoided.
In the actual seismic exploration process, longitudinal wave information is usually used most, the invention uses simpler sound waves to replace elastic waves for forward modeling, and uses the sound waves to research seismic waves, so that the invention is a reasonable approximate treatment for the seismic waves, namely, a sound wave fluctuation equation is used to replace an elastic wave fluctuation equation. In addition, the numerical calculation method based on the difference principle by utilizing the finite difference divides a solving area into a plurality of fine grids, solves a difference equation set on grid points to obtain the numerical solution of the equation, uses the finite difference staggered grid forward modeling to effectively avoid the influence of a boundary zone, and simultaneously conveniently uses the wave spectrum to model the launching sea surface.
Drawings
Fig. 1 is a flow chart of a method for simulating numerical data of an undulating sea surface seismic wave field based on a wave spectrum according to an embodiment of the present invention.
Fig. 2 is a schematic diagram of an interleaved trellis according to an embodiment of the present invention.
FIG. 3 is a diagram of a 30Hz Rake wavelet source waveform provided by an embodiment of the invention.
Fig. 4 is a schematic diagram of a PML absorption boundary provided by an embodiment of the present invention.
Fig. 5 is a schematic diagram of calculating a skew cable according to an embodiment of the present invention.
FIG. 6 is a schematic diagram of the Neumann spectra at different wind speeds provided by the embodiments of the present invention.
FIG. 7 is a schematic diagram of P-M spectra at different wind speeds provided by an embodiment of the present invention.
Fig. 8 is a schematic diagram of ITTC spectra at different wave heights provided by an embodiment of the present invention.
Fig. 9 is a schematic diagram of a two-parameter wave spectrum provided by an embodiment of the invention.
Fig. 10 is a schematic diagram of a swell and rise sea surface based on a wave spectrum provided by the embodiment of the invention.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is further described in detail with reference to the following embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.
Aiming at the problems in the prior art, the invention provides a numerical data simulation method of an undulating sea surface seismic wave field based on a sea wave spectrum, and the invention is described in detail below by combining the attached drawings.
As shown in fig. 1, a method for simulating numerical data of a seismic wave field of an undulating sea surface based on a wave spectrum according to an embodiment of the present invention includes the following steps:
s101: and obtaining a second-order acoustic wave equation through derivation, and performing order reduction processing on the second-order acoustic wave equation in a mode of simultaneously integrating time t at two sides of the equation to obtain an expression of a first-order velocity-stress acoustic wave equation set.
S102: and for the first-order acoustic wave equation obtained in the step S101, performing finite difference calculation on the acoustic wave equations of each grid point and each half grid point in time and space according to the calculation rule of the staggered grids.
S103: forward modeling of marine data; the PML layer above the model calculation region is not provided with an absorption/attenuation coefficient, and the difference calculation of the upper boundary is not performed.
S104: modeling the surge sea surface based on the wave spectrum; and obtaining a space domain model of the surge background undulating sea surface under the P-M spectrum condition by utilizing the P-M spectrum and combining with the inverse Fourier transform.
S105: and applying the obtained numerical value of the undulating sea surface to finite difference numerical simulation, and replacing the numerical value of the upper boundary of the calculation area of the model in the S104 by the numerical value of the undulating sea surface to obtain the simulation data under the surging and undulating sea surface.
The technical solution of the present invention is further described with reference to the following specific examples.
1. Derivation process of acoustic wave equation
Due to the complexity of elastic waves, in current practical seismic surveys, longitudinal wave seismic sources such as explosives, air guns, electric sparks, etc. are often used. In general, an ideal acoustic wave can be modeled numerically for seismic exploration.
The elastic wave equation in the two-dimensional homogeneous medium is derived from a stress equation, a strain equation and a stress-strain relationship, and can be expressed by the following formula:
Figure BDA0002536953790000111
equation (1-1) is an elastic wave equation of a two-dimensional homogeneous medium model, expressed in vector form as:
Figure BDA0002536953790000112
ideally, the shear stress τ is absent in the fluid, and has a value of 0, i.e.:
σxy＝σyz＝σzx＝0 (1-3)
the normal stress at this time is:
σxx＝σyy＝σzz＝-p (1-4)
in the above formula, p represents the average pressure of the fluid medium.
The following relationship exists in an ideal fluid:
σxx＝kθ＝λθ＝-p (1-5)
thus, in an ideal fluid, the elastic wave equation can be simplified as follows:
Figure BDA0002536953790000113
in the formula (1-6), where ρ represents density, divergence is taken simultaneously on both sides of the formula (1-6) and the order of differentiation is exchanged, and it can be obtained:
Figure BDA0002536953790000114
after simplification, the following can be obtained:
Figure BDA0002536953790000115
according to the relational expression
Figure BDA0002536953790000116
The formula (1-8) can be rewritten as
Figure BDA0002536953790000117
Wherein the content of the first and second substances,
Figure BDA0002536953790000118
when this is brought into the above formula (1-9), the acoustic wave equation required by the present invention is obtained, which is expressed as follows:
Figure BDA0002536953790000119
in the formula (1-10), V represents the acoustic velocity, ρ represents the density of the formation medium, and p represents the acoustic pressure.
In the actual forward modeling process, a seismic source function is added to the acoustic wave equation of the above formula to realize the forward modeling of the true acoustic wave equation. The density of the formation is typically assumed to be constant, so that the acoustic equations (1-10) can be expressed in the following relatively common form:
Figure BDA0002536953790000121
the above is the main derivation process of the conventional acoustic wave equation, and the equations (1-11) are acoustic wave equations expressed by pressure and can also be expressed by displacement.
In computer languages, data are discretized into a plurality of points, and in order to realize numerical simulation of seismic waves through programming, the invention needs to discretize the expression of an acoustic wave equation. Obviously, the acoustic wave equation given by the equation (1-11) is a second-order equation, which is inconvenient for discretizing the acoustic wave equation, therefore, before solving the high-order difference equation of the acoustic wave equation, the equation (1-11) should be reduced, and the two sides of the equation are integrated for time t at the same time, so as to obtain the expression of the first-order velocity-stress acoustic wave equation system, as follows:
Figure BDA0002536953790000122
in the formula (1-12), u represents a stress component, vx、vzRespectively representThe velocity components in the x and z directions, ρ represents the density of the formation and v represents the velocity.
2. Finite difference of sound wave equation staggered grid
After the acoustic wave equations of equations (1-12) are obtained, they can be discretized and finite difference calculations performed. The local accuracy of the staggered mesh is greatly improved, and the convergence speed is also increased. The staggered grid of acoustic wave equations is shown in fig. 2.
As can be seen from fig. 2, during the staggered grid calculation, the stress u is placed on each grid point, i.e. at spatial position (i, j); velocity component vxPut on the point that the i direction is an integer and the j direction is a half grid number, namely the space point (i, j + 1/2); similarly, the velocity component vzThe calculation is performed at a point where the i direction is a half mesh point number and the j direction is an integer, that is, a spatial point (i +1/2, j). In the differential calculation of the acoustic wave equation in time, when the calculation is performed by setting discrete points of the stress u on half-time grid points, and corresponding to this, the velocity value vxAnd vzThe calculations are performed at time integer grid points. Next, finite difference calculations need to be performed on each acoustic wave equation of each grid point and each half-grid point in time and space according to the calculation rule of the staggered grid.
(1) Time second order precision staggered grid finite difference
Assuming a velocity component v in the x-directionxIs derivable at an m-th order time, then at some time t is the pair
Figure BDA0002536953790000131
And
Figure BDA0002536953790000132
taylor-series expansion can be carried out to obtain the following formulas (1-13) and (1-14):
Figure BDA0002536953790000133
Figure BDA0002536953790000134
where Δ t is the time grid step, O (Δ t)m+1) Is an infinitesimal quantity of order m +1 with respect to Δ t.
Subtracting the two equations can obtain:
Figure BDA0002536953790000135
if v is paired in timexWhen the high-order difference calculation is carried out, firstly, the calculation amount is increased invisibly, and the calculation efficiency is greatly reduced; secondly, a large number of experiments prove that the influence of high-order time difference precision on numerical simulation is small, and the difference between the high-order difference precision and the second-order difference precision can be ignored, so that too much time does not need to be spent on the order of the time difference, thereby reducing the calculation amount and improving the calculation efficiency. Therefore, the present invention employs temporal second order differential precision. The formulae (1-15) can be simplified as:
Figure BDA0002536953790000136
for the same reason, for the velocity component v in the Z directionzBy performing the same calculation, it is possible to obtain:
Figure BDA0002536953790000137
the stress u is calculated at the point of the time half-grid, i.e. at
Figure BDA0002536953790000138
The Taylor-series expansion is carried out on t and t + delta t, and the same approximation is carried out, so that:
Figure BDA0002536953790000139
the formulae (1-16), the formulae (1-17) and the formulae (1-18) are each vx、vzSecond order difference common of sum u in timeFormula (II) is shown.
(2) Spatial 2N order finite difference format
According to the finite difference calculation rule of the space staggered grids, when calculation is carried out on the space grid points, each order partial derivative of the variable is paired at a space point x
Figure BDA0002536953790000141
And
Figure BDA0002536953790000142
the calculation is performed. Therefore, assume that any one of the functions f (x) has a partial derivative of order 2N +1, respectively
Figure BDA0002536953790000143
And
Figure BDA0002536953790000144
performing Taylor-series expansion:
Figure BDA0002536953790000145
Figure BDA0002536953790000146
Figure BDA0002536953790000147
Figure BDA0002536953790000148
Figure BDA00025369537900001417
Figure BDA0002536953790000149
Figure BDA00025369537900001410
where Δ x represents a spatial grid step length, and the difference between two formulas in the formula (2-21) can be obtained:
Figure BDA00025369537900001411
when the invention wants to calculate the first partial derivative of the function
Figure BDA00025369537900001412
Then, the coefficient may be 1, and the coefficients of the remaining terms may be zero, and in this case, the following formula (1-22) may be written:
Figure BDA00025369537900001413
at this time, the coefficient terms of the equations (1-23) can be written in the form of the following equation set:
Figure BDA00025369537900001414
solving equations (1-24) can be calculated using the Cramer's Law in linear algebra: when N is equal to 1, the compound is,
Figure BDA00025369537900001415
when the N is equal to 2, the N is not more than 2,
Figure BDA00025369537900001416
and so on.
In the seismic forward modeling process, in order to ensure the calculation accuracy and control the computation amount, the difference accuracy of the time 2 order, the space 10 order or the time 2 order, the space 12 order is generally used. The first six even-order difference coefficients, as shown in table 1:
TABLE 1 spatial difference coefficient tables of different orders
Figure BDA0002536953790000151
Dispersing partial derivatives in the formula (1-12) by using the difference method to obtain finite difference formulas of 2 orders of time and 2 Nth orders of space, wherein the finite difference formulas are shown in the formula (1-25):
Figure BDA0002536953790000152
wherein, Δ X and Δ Z respectively represent space grid step lengths in the X and Z directions, and i and j respectively represent corresponding space grid points; Δ t represents a time grid step, k represents a corresponding time grid point; u, V, W show stress u and velocity component vx、vzIn discrete form.
3. Finite difference conditional analysis
3.1 seismic Source Condition analysis
In the forward modeling of seismic waves, commonly used seismic source wavelets include Ricker wavelets, sine attenuating wavelets, Gauss wavelets, Shu wavelets, and the like. The zero-phase rake (Ricker) wavelet is narrower in main lobe and small in side lobe amplitude, is an ideal seismic source wavelet suitable for simulating seismic data, and has the following expression:
Figure BDA0002536953790000153
in the formula (1-26), fmIs the dominant frequency of the wavelet, t0Is the delay time of the seismic source. In order to ensure that the source wavelet can be described by the discretized grid points, the space sampling theorem needs to be satisfied, and if the space sampling is insufficient, serious dispersion is generated, that is, the higher the dominant frequency of the wavelet is, the smaller the grid is. Therefore, in the forward modeling process, the relationship between the grid size and the dominant frequency of the wavelet needs to be considered to satisfy the sampling theorem. The source wavelets used in this paper are all 30Hz zero-phase Rake wavelets, the waveform of which is shown in FIG. 3.
3.2 analysis of stability conditions
The stability condition is a measure of the accuracy of the difference format during forward modeling. In order to ensure the stability of the numerical simulation algorithm and the accuracy of the calculation result, and simultaneously ensure that the amplitude value is in a theoretical range so that the result does not generate large errors, it is necessary to introduce a stability condition.
The stability condition of the interleaved mesh finite difference numerical solution is as follows:
Figure BDA0002536953790000161
when the formula (1-27) is satisfied, the numerical simulation result is accurate.
3.3 frequency Dispersion analysis
The phenomenon that the waveform of seismic waves is changed and different frequency components have different phase velocities is called physical dispersion when the seismic waves are influenced by a propagation medium in the actual propagation process; the dispersion caused by replacing continuous points with discrete grids is called numerical dispersion, the finite difference is a process of putting continuous differential equations on the discretized grids for solving, and the numerical dispersion problem cannot be avoided.
Since the waveform is distorted by the numerical dispersion, the distortion is especially serious when the speed changes severely, even a false in-phase axis is generated, so that the accuracy of forward simulation is greatly reduced, and therefore, the dispersion needs to be suppressed in the forward simulation process. When the sampling theorem is satisfied, the frequency dispersion and the difference order are in negative correlation and in positive correlation with the grid size and the wavelet frequency. Therefore, on the premise of keeping the difference order, the grid size can be reduced, the wavelet dominant frequency is reduced, and the numerical dispersion is suppressed, but the calculation amount is increased by reducing the space grid size.
3.4 analysis of boundary conditions
Seismic waves are generally considered to propagate in an infinite half-space, but in the process of computer numerical simulation, due to various limitations, an infinite model cannot be simulated, and a boundary is inevitably added artificially, so that a boundary effect is caused, and a series of interferences are brought. The scholars propose a plurality of methods to solve the problem of numerical dispersion, such as an extended boundary method, which effectively avoids boundary reflection by enlarging a calculation area, but reduces the calculation efficiency; compared with the extended boundary method, the method for increasing the absorption attenuation effect for the boundary has the advantages of small calculation amount and high operation rate, and can be used for solving the problem of boundary reflection. The present invention employs a Perfect Matching Layer (PML) method, which is proposed by berenser.
The PML is an anisotropic medium essentially, when a wave reaches an interface between the PML and other spaces, the speed and the impedance of the wave are not changed, so that wave reflection can not be generated theoretically, wave transmission is only generated, meanwhile, the PML layer is provided with an absorption attenuation coefficient, the transmitted wave can be absorbed in the medium, and therefore any wave reaching the boundary can not be reflected by the boundary, and the effect of PML on processing the boundary problem is obvious. The PML perfect match layer is shown in fig. 4.
In fig. 4, the blue squares represent the model calculation regions, the remainder being the PML layers, assuming d (x) and d (z) represent the absorption attenuation parameters in the x and z directions, respectively. Both d (x) and d (z) are not 0 in region 1; d (x) is 0, d (z) is not 0 in region 2; in region 3, d (x) is not 0, and d (z) is 0. Therefore, the seismic wave does not reflect when propagating to the boundary, the energy of the seismic wave is absorbed and attenuated in the PML layer, and the absorption and attenuation functions are as follows:
Figure BDA0002536953790000171
in the formula (1-28), R is a theoretical reflection coefficient, vpIs the longitudinal wave velocity and is the thickness of the PML layer.
The derivation of the difference equation after the addition of the PML layer is similar to that before, and not too much derivation is made here, but only the difference equation with the PML added is given:
Figure BDA0002536953790000172
Figure BDA0002536953790000173
Figure BDA0002536953790000174
Figure BDA0002536953790000175
Figure BDA0002536953790000176
4. forward modeling of marine data
For marine exploration, the sea surface is a strong reflection interface, so in forward simulation of marine seismic exploration, the PML layer above the model calculation region is not provided with an absorption attenuation coefficient, and differential calculation of the upper boundary is not performed, so that the wave is reflected when propagating to the upper boundary and is absorbed and attenuated by the PML layer when propagating to other boundaries.
The conventional horizontal streamer acquisition is most widely applied, the sinking depth of a detector is fixed in the acquisition process, and the position of the detector is only required to be arranged on a row of fixed grid points in the forward simulation process.
For the variable depth cable, the depth of the detector increases along with the increase of the offset distance, the position of the detector may not be on a grid point, and the finite difference method can only calculate the wave field value of the grid point. Therefore, in the simulation of the variable depth cable, if the detector is not located at the grid point, the wave field value at the detector needs to be calculated by interpolating the wave fields at two grid points above and below the detector, as shown in fig. 5.
In order to calculate the wave field of the detector, it is necessary to interpolate the wave field value U (i, j) of the upper point of the detector and the wave field value U (i, j +1) of the lower point of the detector, so as to obtain the wave field value of the intermediate detector. If instead of interpolation, the wavefield values of the grid points closest to the detector are substituted, the seismic recording will exhibit a step-wise discontinuity, resulting in a discontinuity in the in-phase axis.
The interpolation mode adopted by the inclined cable simulation in the invention is as follows: and calculating weights according to the distances between the ith channel detector (ith column) and two adjacent grid points (i, j) and (i, j +1), and performing interpolation according to the weights to obtain the wave field value at the detector.
5. Common wave spectrum
The sea surface changes due to the influence of environmental factors, the sea water moves and changes at any moment, the influence of wind is the most huge, and the sea surface has the characteristic of random heaving because the instantaneous values of wind direction and wind speed are random. The fluctuation is difficult to describe by classical fluctuation theory and to define and describe by a certain deterministic function.
Many scholars have done relevant research on how to describe the important problem of undulating sea. In the marine science, it is generally regarded as a random process, in which the total variation range of the statistical properties, such as average state, amplitude (or wave height) and other elements, is studied, wherein there is a certain regularity, and for this reason, the concept of spectrum is introduced, forming the most commonly used wave spectrum in the marine science.
As early as the 50 s of the 20 th century, ocean waves were considered to be a superposition of many simple fluctuations that differ in amplitude, frequency, direction, and phase. The researchers studied the waves from this viewpoint and stipulate the amplitude or phase of these simple fluctuations as a random quantity, so the superposition result is also a random quantity.
The energy provided by the component waves with frequency f is represented by s (f), which is called the wave spectrum. The sea wave spectrum s (f) represents the distribution of the sea wave energy with respect to the frequency f, which is essentially a function describing the sea wave energy in a limited area and is also the power density spectrum of the sea waves. And multiplying a set of Gaussian-distributed random numbers by the wave spectrum function to obtain the undulating sea surface under the wave spectrum.
Common ocean wave spectra mainly include a Neumann spectrum (Noemann spectrum), a Pierson-Moskowitz spectrum (P-M spectrum), an ITTC spectrum and a two-parameter spectrum.
(1) Neumann spectrum
The Neumann spectrum is a sea wave spectrum proposed in 1952, and is a semi-empirical semi-theoretical formula deduced after observing the relation between wave height and period under different wind speed conditions and giving a series of assumed conditions. Although early relatively crude information was used and some controversial assumptions were given in the formulation derivation, it was practical to the extent that some engineering problems were still in use to date. The formula is as follows:
Figure BDA0002536953790000191
wherein U is the wind speed at 7.5 meters above the sea surface; c is a constant of 3.05m/s 2. The Neumann spectrum contains only one variable, namely the wind speed U. Fig. 6 is a Neumann wave spectrum curve at different wind speeds. From this it can be found that: theoretically, the wave spectrum contains all frequencies from zero to infinity, but the main components of the spectrum are concentrated in a narrow low-frequency range; secondly, with the increase of the wind speed, the total area under the spectrum curve is increased, the frequency band range is also increased, and the corresponding wave height and period are also increased; and thirdly, with the increase of the wind speed, the spectrum can move to the low-frequency direction.
(2) Pierson-Moskowitz (P-M) spectra
In 1964, Moskowitz performed 460 more spectral analyses and fitting on the Atlantic observations from 1955 to 1960, and obtained a semi-empirical theoretical model, i.e., P-M spectrum. The formula is as follows:
Figure BDA0002536953790000194
where g is the gravitational acceleration, the constant α is 0.0081, β is 0.74, U is the wind speed 19.5 meters above the sea surface, and the variable f is the frequency.
The P-M spectrum has a more sufficient observation data base and a more effective analysis method than the Neumann spectrum, is more convenient to use, replaces the Neumann spectrum in many occasions, and the P-M spectrum with the wind speed of 20M/s, 15M/s and 10M/s is shown in figure 7.
As can be seen from fig. 7, the larger the wind speed, the larger the frequency range involved, and the larger the energy of the wind waves; meanwhile, as the wind speed increases, the spectrum gradually moves to the low frequency direction. Comparing fig. 6 and fig. 7, it can be found that: the Neumann spectrum and the P-M spectrum are very similar. However, the accuracy of the data on which the Neumann spectrum is based is not high and the precision is low, and although the spectrum is consistent with the reality in a certain range, the method still has great controversy. Besides a few engineering problems still using the Noemann spectrum, the P-M spectrum is used for oceanographic calculation in most cases, and is also the most widely used wave spectrum at present. The invention also adopts the P-M spectrum which is widely used for carrying out the surge modeling.
(3) ITTC Spectrum
In 1972, an international towing tank conference (ITTC for short) modified a P-M spectrum to obtain an ITTC spectrum. According to the P-M spectrum, there are:
Figure BDA0002536953790000201
because of the following: 4 (m)0)1/2Namely:
Figure BDA0002536953790000202
therefore, the method comprises the following steps:
Figure BDA0002536953790000203
the ITTC spectrum can be obtained by substituting 0.0081 a, 0.78 g2 and B in the P-M spectrum into equation (3)
Figure BDA0002536953790000204
Among these, the average wave height is shown. Fig. 8 shows ITTC spectra for different wave heights, setting the wave heights to 4, 3 and 2 meters, respectively.
(4) Two parameter spectra
In 1978, at the 15 th ITTC conference, the ITTC spectrum was further optimized and improved, and a two-parameter ocean wave spectrum was derived, which is formulated as follows:
Figure BDA0002536953790000205
wherein, T is 5.127/B1/4, A is 1732/T4, and the calculation method of B is the same as ITTC spectrum. Fig. 9 shows a two-parameter wave spectrum with wave heights of 4 meters, 3 meters and 2 meters, respectively.
The above four wave spectrums are commonly used in physical oceanography, and are widely used in research of physical oceanography, and are briefly introduced. The P-M spectrum is selected for simulating the fluctuating sea surface of the surge background, and although the sea surface is random in a large range, the fluctuating sea surface simulated by the wave spectrum is more in line with the actual situation.
6. Surge sea surface modeling based on sea wave spectrum
The mathematical formula of the P-M spectrum is:
Figure BDA0002536953790000211
where α is 0.0081, β is 0.74, both of which are constants, g is the gravitational acceleration, ω is the circular frequency, and U is the wind speed above the sea surface, which is an expression of the wave spectrum in the frequency domain. In water, the following relationship exists between frequency and wavenumber:
ω2(k)＝gk (1-37)
the expression of the P-M spectrum in the wavenumber domain is:
Figure BDA0002536953790000212
after obtaining the equations 1-38, the equations are multiplied by a set of gaussian distributed random complex numbers and subjected to inverse fourier transform, so as to obtain a spatial domain model of the swell background undulating sea surface under the P-M spectrum condition, as shown in fig. 10.
In fig. 10, the abscissa represents the distance in the horizontal direction and the ordinate represents the wave height. The sea surface satisfies the distribution of the P-M wave spectrum, and better conforms to the actual sea surface condition than the common random number.
After the fluctuating sea surface numerical value is obtained by the method, the numerical value is applied to finite difference numerical simulation, and the numerical value of the upper boundary of the model calculation area is replaced by the series of numerical values, so that the simulation data under the surging and fluctuating sea surface can be obtained. When numerical simulation is performed, the size of the grid points cannot be a decimal, but the wave spectrum obtained by the method can be a decimal, so that the sea surface fluctuation value needs to be rounded.
The above description is only for the purpose of illustrating the preferred embodiments of the present invention and is not to be construed as limiting the invention, and any modifications, equivalents and improvements made within the spirit and principle of the present invention are intended to be included within the scope of the present invention.

Claims (8)

1. A numerical data simulation method of a fluctuating sea surface seismic wave field based on a sea wave spectrum is characterized by comprising the following steps:
step one, obtaining a second-order acoustic wave equation through derivation, and performing order reduction processing on the second-order acoustic wave equation in a mode of integrating time t at two sides of the equation at the same time to obtain an expression of a first-order speed-stress acoustic wave equation set;
step two, performing finite difference calculation on the first-order acoustic wave equation obtained in the step one according to the calculation rule of staggered grids on the acoustic wave equations of each grid point and each half grid point in time and space;
step three, forward modeling of marine data: the PML layer above the model calculation area is not provided with an absorption attenuation coefficient, and meanwhile, the difference calculation of an upper boundary is not carried out;
step four, modeling the surge sea surface based on the wave spectrum: obtaining a space domain model of the surge background undulating sea surface under the P-M spectrum condition by utilizing the P-M spectrum in combination with the inverse Fourier transform;
and step five, applying the obtained numerical value of the undulating sea surface to finite difference numerical simulation, and replacing the numerical value of the upper boundary of the calculation area of the model in the step four with the numerical value of the undulating sea surface to obtain simulation data under the surging and undulating sea surface.
2. The numerical data simulation method of the wave field of the undulating sea surface based on the wave spectrum of the sea as set forth in claim 1, wherein in the step one, the elastic wave equation in the two-dimensional uniform medium is derived according to the stress equation, the strain equation and the relationship between the stress and the strain, and is expressed by the following formula:
Figure FDA0002536953780000011
the formula is an elastic wave equation of a two-dimensional homogeneous medium model, and is expressed in a vector form as:
Figure FDA0002536953780000012
shear stress τ is absent in the fluid and has a value of 0:
σxy＝σyz＝σzx＝0；
the normal stress is:
σxx＝σyy＝σzz＝-p；
wherein p represents the average pressure of the fluid medium;
the following relationships exist in the fluid:
σxx＝kθ＝λθ＝-p；
in the fluid, the elastic wave equation is processed as follows:
Figure FDA0002536953780000021
where ρ represents density, taking the divergence for both sides of the formula and exchanging the differential order, yields:
Figure FDA0002536953780000022
and (3) after simplification:
Figure FDA0002536953780000023
according to the relational expression
Figure FDA0002536953780000024
Then will formula
Figure FDA0002536953780000025
Rewriting to;
Figure FDA0002536953780000026
wherein the content of the first and second substances,
Figure FDA0002536953780000027
bring it into the above formula
Figure FDA0002536953780000028
Obtaining an acoustic wave equation, wherein the expression is as follows:
Figure FDA0002536953780000029
in the formula, V represents the acoustic velocity, ρ represents the density of the formation medium, and p represents the acoustic pressure;
in the actual forward modeling process, a seismic source function needs to be added into the acoustic wave equation of the above formula, and the acoustic wave equation is assumed that the density value of the stratum is constant
Figure FDA00025369537800000210
Expressed in the following form:
Figure FDA00025369537800000211
in pair type
Figure FDA00025369537800000212
And (3) carrying out order reduction treatment, and integrating the time t at two sides of the equation to obtain an expression of a first-order velocity-stress sound wave equation set:
Figure FDA00025369537800000213
wherein u represents a stress component, vx、vzThe velocity components in the x and z directions are represented, respectively, with ρ representing the density of the formation and v representing the velocity.
3. The method for numerical data simulation of a seismic wave field at undulating sea based on sea wave spectrum according to claim 1, wherein in step two, the time second-order precision staggered grid finite difference comprises: velocity component v in the x-directionxIs derivable at an m-th order time, then is paired at some time t
Figure FDA0002536953780000031
And
Figure FDA0002536953780000032
and carrying out Taylor-series expansion to obtain:
Figure FDA0002536953780000033
Figure FDA0002536953780000034
where Δ t is the time grid step, O (Δ t)m+1) Is an infinitesimal quantity of order m +1 with respect to Δ t;
formula (II)
Figure FDA0002536953780000035
And formula
Figure FDA0002536953780000036
Subtracting the two equations to obtain:
Figure FDA0002536953780000037
the time second-order difference precision is adopted, and the formula is simplified as follows:
Figure FDA0002536953780000038
for the same reason, for the velocity component v in the Z directionzAnd calculating to obtain:
Figure FDA0002536953780000039
the stress u is calculated at the point of the time half-grid, and the stress u is calculated at
Figure FDA00025369537800000310
And performing Taylor series expansion on t and t + delta t, and performing the same approximation to obtain:
Figure FDA00025369537800000311
formula (II)
Figure FDA00025369537800000312
Formula (II)
Figure FDA00025369537800000313
And formula
Figure FDA00025369537800000314
Are respectively vx、vzAnd u is a second order difference formula over time.
4. The method for numerical data simulation of a seismic wave field at undulating sea based on sea wave spectrum according to claim 1, wherein in step two, the spatial 2N order finite difference format comprises: any function f (x) having a partial derivative of order 2N +1, respectively
Figure FDA00025369537800000315
And
Figure FDA00025369537800000316
performing Taylor-series expansion:
Figure FDA00025369537800000317
Figure FDA0002536953780000041
wherein, Δ x represents the space grid step length, and the difference between the two formulas is obtained:
Figure FDA0002536953780000042
calculating the first partial derivative of the formula
Figure FDA0002536953780000043
Let its coefficient be 1 and the coefficients of the remaining terms be zero, and then write as follows:
Figure FDA0002536953780000044
at this time, the coefficient terms of the equations are written in the form of the following equation set
Figure FDA0002536953780000045
The solution equation is calculated using the claime's law in linear algebra: when N is equal to 1, the compound is,
Figure FDA0002536953780000046
when the N is equal to 2, the N is not more than 2,
Figure FDA0002536953780000047
dispersing each partial derivative by using a difference method to obtain finite difference formulas of 2 orders of time and 2 Nth orders of space:
Figure FDA0002536953780000048
wherein, Deltax, Deltaz represent space grid step length in X and Z direction separately, i, j represent corresponding separatelyΔ t represents a time grid step, k represents the corresponding time grid point, U, V, W represents the stress u, velocity component v, respectivelyx、vzIn discrete form.
5. The method for simulating numerical data of a seismic wave field in an undulating sea surface based on a wave spectrum as claimed in claim 1, wherein in the third step, when the variable depth cable is simulated, if the geophone is not located at a grid point, the wave field value at the geophone is calculated by interpolation of wave fields of two grid points above and below the geophone;
the interpolation mode adopted by the skew cable simulation is as follows: and calculating the weight according to the distance between the ith detector and two adjacent grid points (i, j) and (i, j +1), and performing interpolation according to the weight to obtain the wave field value at the detector.
6. The numerical data simulation method of the seismic wavefield from undulating sea surface based on wave spectrum of claim 1, wherein in step four, the mathematical formula of the P-M spectrum is as follows:
Figure FDA0002536953780000051
wherein α is 0.0081, β is 0.74, both of which are constants, g is gravity acceleration, ω is circular frequency, and U is wind speed over the sea surface, which is an expression of a wave spectrum in a frequency domain; in water, the following relationship exists between frequency and wavenumber:
ω2(k)＝gk；
the expression of the P-M spectrum in the wavenumber domain is:
Figure FDA0002536953780000052
in obtaining a formula
Figure FDA0002536953780000053
Then, the formula is multiplied by a group of Gaussian random complex numbers and subjected to inverse Fourier transform to obtain the productA space domain model of a surge background fluctuating sea surface under the P-M spectrum condition; and after the fluctuating sea surface numerical value is obtained, applying the numerical value to finite difference numerical simulation, and replacing the numerical value of the upper boundary of the model calculation area with the series of numerical values to obtain simulation data under the surging and fluctuating sea surface.
7. The method according to claim 6, wherein the numerical data of the sea surface heave value is rounded during the numerical simulation.
8. The method for numerical data simulation of a seismic wavefield from an undulating sea surface based on a wave spectrum of claim 1, wherein the finite difference conditions include source conditions, stability conditions, dispersion, boundary conditions;
1) the seismic source wavelet of the analog seismic data selects a 30Hz zero-phase Rake wavelet, and the expression is as follows:
Figure FDA0002536953780000061
in the formula (f)mIs the dominant frequency of the wavelet, t0Is the delay time of the seismic source;
2) the stability condition of the interleaved mesh finite difference numerical solution is as follows:
Figure FDA0002536953780000062
3) maintaining differential order frequency dispersion analysis;
4) by adopting a PML method of a complete matching layer, seismic waves cannot be reflected when propagating to a boundary, the energy of the seismic waves can be absorbed and attenuated in the PML layer, and the absorption and attenuation function is as follows:
Figure FDA0002536953780000063
wherein R is the theoretical reflection coefficient, vpIs the velocity of longitudinal wavesDegree, is the PML layer thickness;
the difference equation added to PML is:
Figure FDA0002536953780000064
Figure FDA0002536953780000065
Figure FDA0002536953780000066
Figure FDA0002536953780000067
Figure FDA0002536953780000068
CN202010535639.5A 2020-06-12 2020-06-12 Numerical data simulation method for undulating sea surface seismic wave field based on sea wave spectrum Active CN111797552B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010535639.5A CN111797552B (en) 2020-06-12 2020-06-12 Numerical data simulation method for undulating sea surface seismic wave field based on sea wave spectrum

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010535639.5A CN111797552B (en) 2020-06-12 2020-06-12 Numerical data simulation method for undulating sea surface seismic wave field based on sea wave spectrum

Publications (2)

Publication Number Publication Date
CN111797552A true CN111797552A (en) 2020-10-20
CN111797552B CN111797552B (en) 2022-05-27

Family

ID=72803291

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010535639.5A Active CN111797552B (en) 2020-06-12 2020-06-12 Numerical data simulation method for undulating sea surface seismic wave field based on sea wave spectrum

Country Status (1)

Country Link
CN (1) CN111797552B (en)

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112749442A (en) * 2020-12-25 2021-05-04 青岛黄海学院 Automobile seismic source near-surface wave simulation algorithm
CN112883564A (en) * 2021-02-01 2021-06-01 中国海洋大学 Water body temperature prediction method and prediction system based on random forest
CN112882093A (en) * 2021-01-18 2021-06-01 中国测绘科学研究院 Method and system for calculating internal homoseismic deformation of elastic earth
CN113484914A (en) * 2021-07-13 2021-10-08 中海石油(中国)有限公司 Method, system, medium and equipment for manufacturing marine storm wave consistency influence gauge plate
CN113552628A (en) * 2021-07-20 2021-10-26 黄河勘测规划设计研究院有限公司 Earthquake surge height calculation method
CN114545505A (en) * 2022-01-12 2022-05-27 吉林大学 Dynamic ocean wave field information acquisition method based on particle grid thought
CN114745046A (en) * 2022-03-16 2022-07-12 中国科学院西安光学精密机械研究所 Method for analyzing pointing deviation of laser beam emitted from random fluctuation sea surface
CN116108720A (en) * 2023-02-17 2023-05-12 国家海洋环境预报中心 Wave forecasting method and system based on wave numerical mode of SCVT grid
CN116187064A (en) * 2023-02-14 2023-05-30 中国科学院国家空间科学中心 Numerical simulation method for second derivative of continuous signal time sequence

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101788683A (en) * 2009-12-29 2010-07-28 华东师范大学 Tsunami motion forecasting method based on multi-hierarchy interaction
US20130208565A1 (en) * 2012-02-09 2013-08-15 Pgs Geophysical As Methods and systems for correction of streamer-depth bias in marine seismic surveys

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101788683A (en) * 2009-12-29 2010-07-28 华东师范大学 Tsunami motion forecasting method based on multi-hierarchy interaction
US20130208565A1 (en) * 2012-02-09 2013-08-15 Pgs Geophysical As Methods and systems for correction of streamer-depth bias in marine seismic surveys

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
KUANGDAI LENG 等: "AxiSEM3D: broad-band seismic wavefields in 3-D global earth", 《GEOPHYSICAL JOURNAL INTERNATIONAL》 *
王春燕: "高阶交错网格有限差分地震波场计算", 《成都理工大学硕士学位论文》 *
赵天亮: "地震波的有限差分法正演模拟", 《地下水》 *

Cited By (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112749442A (en) * 2020-12-25 2021-05-04 青岛黄海学院 Automobile seismic source near-surface wave simulation algorithm
CN112882093A (en) * 2021-01-18 2021-06-01 中国测绘科学研究院 Method and system for calculating internal homoseismic deformation of elastic earth
CN112882093B (en) * 2021-01-18 2024-03-05 中国测绘科学研究院 Method and system for calculating internal co-vibration deformation of elastic earth
CN112883564A (en) * 2021-02-01 2021-06-01 中国海洋大学 Water body temperature prediction method and prediction system based on random forest
CN113484914A (en) * 2021-07-13 2021-10-08 中海石油(中国)有限公司 Method, system, medium and equipment for manufacturing marine storm wave consistency influence gauge plate
CN113484914B (en) * 2021-07-13 2023-09-12 中海石油(中国)有限公司 Method, system, medium and equipment for manufacturing marine storm consistency influence measuring plate
CN113552628B (en) * 2021-07-20 2023-08-15 黄河勘测规划设计研究院有限公司 Method for calculating height of earthquake wave
CN113552628A (en) * 2021-07-20 2021-10-26 黄河勘测规划设计研究院有限公司 Earthquake surge height calculation method
CN114545505A (en) * 2022-01-12 2022-05-27 吉林大学 Dynamic ocean wave field information acquisition method based on particle grid thought
CN114745046A (en) * 2022-03-16 2022-07-12 中国科学院西安光学精密机械研究所 Method for analyzing pointing deviation of laser beam emitted from random fluctuation sea surface
CN114745046B (en) * 2022-03-16 2023-09-01 中国科学院西安光学精密机械研究所 Method for analyzing pointing deviation of laser beam emitted from randomly-fluctuated sea surface
CN116187064A (en) * 2023-02-14 2023-05-30 中国科学院国家空间科学中心 Numerical simulation method for second derivative of continuous signal time sequence
CN116187064B (en) * 2023-02-14 2024-03-12 中国科学院国家空间科学中心 Numerical simulation method for second derivative of continuous signal time sequence
CN116108720A (en) * 2023-02-17 2023-05-12 国家海洋环境预报中心 Wave forecasting method and system based on wave numerical mode of SCVT grid
CN116108720B (en) * 2023-02-17 2023-08-25 国家海洋环境预报中心 Wave forecasting method and system based on wave numerical mode of SCVT grid

Also Published As

Publication number Publication date
CN111797552B (en) 2022-05-27

Similar Documents

Publication Publication Date Title
CN111797552B (en) Numerical data simulation method for undulating sea surface seismic wave field based on sea wave spectrum
CN111239802B (en) Deep learning speed modeling method based on seismic reflection waveform and velocity spectrum
CN101017204B (en) Three-dimensional two-way acoustic wave equation pre-stack imaging systems and method
CN104122585B (en) Seismic forward simulation method based on elastic wave field resolution of vectors and low-rank decomposition
Simmen et al. Wavefront folding, chaos, and diffraction for sound propagation through ocean internal waves
Liu et al. A mixed-grid finite element method with PML absorbing boundary conditions for seismic wave modelling
CN106873038B (en) A method of extracting Depth Domain seismic wavelet from Depth Domain seismic data
CN108594302B (en) A kind of extracting method and processing terminal of seismic wavelet
Lu et al. Optimal 15‐Point Finite Difference Forward Modeling in Frequency‐Space Domain
Lisitsa et al. Finite-difference algorithm with local time-space grid refinement for simulation of waves
CN108108331B (en) A kind of finite difference formulations method based on quasi- spatial domain equations for elastic waves
CN108051855B (en) A kind of finite difference formulations method based on quasi- spatial domain ACOUSTIC WAVE EQUATION
CN104459773A (en) Unconditionally stable seismic wave field continuation method based on staggered grid Lowrank decomposition
CN108508482A (en) A kind of subterranean fracture seismic scattering response characteristic analogy method
CN106772585A (en) Analysis method and device is intended in a kind of optimization based on elastic wave decoupling equation
Salatin et al. Effects of wave coherence on longshore variability of nearshore wave processes
CN109239776B (en) Seismic wave propagation forward modeling method and device
CN107102359A (en) Geological data protects width method for reconstructing and system
CN110146923A (en) A kind of efficient high accuracy depth domain methods of seismic wavelet extraction
CN105676280A (en) Two-phase medium geological data obtaining method and device based on rotationally staggered grids
Cao et al. Accelerating 2D and 3D frequency-domain seismic wave modeling through interpolating frequency-domain wavefields by deep learning
CN113552633B (en) Elastic wave frequency dispersion pressing method for optimizing difference coefficient and longitudinal and transverse wave separation FCT
Wolfson et al. Full-wave simulation of the forward scattering of sound in a structured ocean: A comparison with observations
Wang et al. 3-D acoustic wave equation forward modeling with topography
Li et al. Inversion of internal wave-perturbed sound-speed field via acoustic data assimilation

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant