Water body wave phenomenon identification method based on wave propagation ray theory simulation in seawater

Abstract

The application provides a water body body wave phenomenon identification method based on wave propagation ray theory simulation in seawater, and the method comprises the following steps: step 1, obtaining the depth values of all OBN nodes on a single line; step 2, simulating the water body propagation wave field through a trial shooting method; step 3, realizing the tracking identification of the wave field based on Markov optimal decision; and step 4, outputting the water body wave field travel time acquisition result. The water body body wave phenomenon identification method based on wave propagation ray theory simulation in seawater is simple in technical operation, overcomes the problem that the water body wave field is difficult to identify, helps to extract the wave field caused by the water layer alone, and is used for fine modeling of the seabed surface form; lays a foundation for subsequent establishment of a fine water body velocity model, and provides a reference for OBN seismic data processing and interpretation work.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of seismic data processing technology of seismic exploration, and particularly relates to a water body wave phenomenon identification method based on wave propagation ray theory simulation in seawater. BACKGROUND

[0002] Due to the existence of the two water body strong reflection interfaces of seawater-air interface and seawater-seabed interface, the water body wave field appears as strong energy in the OBN record, the direct wave, the multi-order water body related multiple waves and the effective reflection wave field from below the water body interface are intertwined together, and the multiple reflection and refraction strong energy interference from the water body multiple oscillation also appears at the far offset, which brings great difficulty to the subsequent processing. The complex wave phenomenon caused by the complex shallow water layer is one of the important interference waves, and if the water body wave phenomenon is identified, it is beneficial to the fine tomographic inversion of the water body model, and it is also convenient to extract the wave field caused by the water layer alone and highlight the deep reflection wave. However, the water body wave phenomenon identification problem has less research at present.

[0003] In the Chinese patent application with the application number: CN202010535639.5, a fluctuating sea surface seismic wave field numerical data simulation method based on sea wave spectrum is involved, a first-order velocity-stress acoustic wave equation set expression is obtained through derivation and order reduction processing; the acoustic wave equation items of each grid point and half grid point in time and space are respectively calculated according to the calculation rules of the staggered grid; offshore data forward modeling; the spatial domain model of the swell background fluctuating sea surface under the P-M spectrum condition is obtained; the obtained fluctuating sea surface numerical value is applied to the finite difference numerical simulation, and the sea surface numerical value is used to replace the upper boundary numerical value of the above-mentioned model calculation area, so that the simulation data under the fluctuating sea surface of the swell can be obtained. The numerical simulation method equation provided by the invention is simple to use, efficient, and the precision can be improved by certain means; the three-dimensional seismic data acquisition construction cost is low, the source and the receiver are sunk to a moderate depth, and the occurrence of false events can be effectively avoided.

[0004] In the Chinese patent application No. CN201811292058.2, a method for quickly establishing a three-dimensional near-seabed velocity model in a shallow beach area is disclosed, which comprises: extracting first arrival information and extracting common midpoint gathers; calculating the slope of the first arrival travel time in the common midpoint gathers; determining whether the shot point is located in the seawater layer, and if so, correcting the shot point to the seabed; calculating the velocity and depth of the turning point corresponding to each offset distance; performing linear interpolation on the velocity at each depth to obtain a one-dimensional velocity model of the common midpoint; after processing all common midpoint gathers, using linear interpolation method, using the one-dimensional velocity models of all common midpoints to calculate the velocity values of the remaining grid points, and outputting the velocity model. The method for quickly establishing a three-dimensional near-seabed velocity model in a shallow beach area overcomes the influence of the seawater layer on ray propagation, eliminates the influence of multiple excitation and reception modes on the inversion result, and avoids the problem of requiring multiple ray tracing forward modeling in the conventional method.

[0005] In the Chinese patent application No. CN201910816102.3, a method and device for collecting ocean data and a storage medium are disclosed to solve the technical problem of how to optimize the OBN observation system to make the collected ocean data meet the requirements of oil and gas exploration and development in related technologies. The method for collecting ocean data comprises: constructing a three-dimensional block geological model according to the geological measured data of the work area to be collected; simulating the layout of multiple OBN observation systems in the three-dimensional block geological model; simulating Gaussian ray irradiation from the shot point as the source point to the target layer in the three-dimensional block geological model; obtaining the bin migration energy analysis result of the target layer according to the illumination energy of the Gaussian ray irradiation received by the geophone in the simulated OBN observation system and the illumination energy reflected by the target layer; selecting a target OBN observation system from the multiple OBN observation systems according to the bin migration energy analysis result; and deploying collection equipment to the work area to be collected according to the target OBN observation system to collect ocean data.

[0006] In the Chinese patent application No. CN202210592015.6, a two-dimensional sea surface ghost wave water body imaging measurement method, system, terminal and flow measuring equipment are disclosed. The frequency domain forward transform data of the seismic data is obtained; for any frequency component and ray parameter, the delay time of the virtual reflection and the first reflection is calculated to obtain the first downlink wave field frequency domain Radon transform result; the average delay time and the average amplitude difference coefficient of the virtual reflection and the first reflection corresponding to the ray parameter are obtained and calculated again, and the iteration is alternately performed until the convergence condition is met to obtain the final frequency domain Radon transform result of the virtual reflection, and then the frequency domain Radon inverse transform and the time domain Fourier inverse transform are performed to obtain the final virtual reflection prediction data. The invention can perform full waveform inversion on the useful information such as the virtual reflection signal received by the ADCP to obtain the property information of the surface seawater, and accurately image the water body.

[0007] The above prior art is quite different from the present application, and cannot solve the technical problems we want to solve, therefore we have invented a new water body wave phenomenon identification method based on wave propagation ray theory simulation in seawater. SUMMARY

[0008] The purpose of the present application is to provide a water body wave phenomenon identification method based on wave propagation ray theory simulation in seawater, which simulates the propagation path of the wave field in the water body by trial shooting method, and realizes the tracking identification of the water body wave field based on Markov optimal decision under model constraints.

[0009] The purpose of the present application can be achieved by the following technical measures: a water body wave phenomenon identification method based on wave propagation ray theory simulation in seawater, which comprises:

[0010] Step 1, obtaining the depth values of all OBN nodes on a single line;

[0011] Step 2, simulating the water body propagation wave field by trial shooting method;

[0012] Step 3, realizing the tracking identification of the wave field based on Markov optimal decision;

[0013] Step 4, outputting the water body wave field travel time acquisition result.

[0014] The purpose of the present application can also be achieved by the following technical measures:

[0015] In step 1, OBN common receiver gather reading is performed, the trace header information of P component data is obtained, and the depth values of all OBN nodes on a single line are obtained.

[0016] Step 1 comprises:

[0017] A1.1: OBN shot line reading is performed, the trace header information is obtained, including line number, shot coordinates, receiver coordinates, elevation, offset distance, sampling rate, and sampling time;

[0018] A1.2: The depth values of all OBN nodes on a single line are recorded.

[0019] In step 2, the water body model seabed interface is expressed in the form of an analytical function according to the OBN node depth values; and the theoretical first-order and second-order water body related multiple wave arrival times of each shot and each OBN node are calculated by trial shooting method based on the initial seabed model.

[0020] Step 2 comprises:

[0021] A2.1: The interface of the water body model is expressed in the form of an analytical function according to the OBN node depth values;

[0022] A2.2: Given the ray starting point position and assuming a ray exit direction, the ray will be reflected according to the geometric properties after reaching the interface, thus tracing a ray, calculating its intersection point when reflected at the interface and its final landing point at the interface;

[0023] A2.3: Adjust the initial exit direction of the ray according to a certain exit angle interval, and find the left and right two initial exit angles of the ray closest to the target receiving point coordinates;

[0024] A2.4: Further narrow the angle range of the target ray by bisection method, so that the calculated theoretical exit ray reaches the interface intersection point coordinates approaching the coordinates of the target receiving point, and records the intersection point when the ray passes through the interface;

[0025] A2.5: Assuming a constant speed background, calculate the travel time of the wave reaching the receiving point according to the length of the ray path;

[0026] A2.6: Get the theoretical time-distance relationship of the first and second order water wave field on the common receiver gather by data rearrangement.

[0027] In step 2.1, for a two-dimensional geological body, the coordinates and depth of each OBN node are known, and the depth value of any position on the interface can be calculated by linear interpolation, polynomial interpolation and cubic spline function interpolation. Assuming that there are n+1 OBN nodes, which are divided into n intervals, the depth of a point x∈[x k ,x k+1 ] in a certain interval is described by a cubic spline function:

[0028] y=a i +b i x+c i x 2 +d i x 3

[0029] Where y is the depth of the interpolation point, a i ,b i ,c i ,d i are the constant term, linear term coefficient, quadratic term coefficient and cubic term coefficient of the interpolation function respectively

[0030] All points of the cubic spline function need to satisfy the interpolation condition: f i (x i )=y i (i=0,1,...,n), and secondly the curve needs to be smooth, that is, f i (x i ),f i '(x i ),f i ”(xi ) continuous; therefore, all interior endpoints need to satisfy the cubic equation for both left and right segments; n-1 interior points first derivative continuous equation: f i '(x i ) = f i ' +1 (x i ), f i ”(x i ) = f i ” +1 (x i ); assuming natural boundary conditions, i.e. f0”(x0) = f n ”(x n ) = 0, then a total of 4n equations are composed, by derivation to get linear equations:

[0031]

[0032] where h i = x i+1 - x i , m i = f i ”(x i ), h i denotes the step size, and m i is the second derivative at the node. Construct linear equations

[0033]

[0034] Solve by Jacobi iteration method, and then substitute m back into the coefficients of the cubic spline function; there are the following formulas:

[0035]

[0036] For three-dimensional geological bodies, the interface is described by triangular patches or B-spline surfaces.

[0037] In step 2.2, given a ray starting point position and the initial direction of the ray, the downgoing wave process and the upgoing wave process of the ray propagation path are calculated; given the downgoing ray slope k1 and intercept b1, the intersection coordinates of the downgoing ray and the water body model seabed interface are calculated by simultaneous equations;

[0038]

[0039] where x2, y2 are the coordinates of the interpolation point, a i , b i , c i , d i are the constant term, linear term coefficient, quadratic term coefficient, and cubic term coefficient of the interpolation function, respectively.

[0040] Given the ray incidence point coordinates and the interface intersection point coordinates, the direction vector of the incident ray is determined according to the slope k1 of the incident ray The direction vector of the outgoing ray after reflection also needs to be determined According to the law of reflection, the incidence angle is equal to the emergence angle. The direction vector of the incident ray and the direction vector of the reflected ray are symmetrical with respect to the normal vector. The direction vector of the outgoing ray is determined as follows:

[0041]

[0042] Wherein is the local slope of the interface at the intersection point, the ray incidence point coordinates are (x1, y1), the interface intersection point coordinates are (x2, y2), and the direction of the incident ray is The direction of the outgoing ray is

[0043] According to the local slope of the interface, the slope of the upward outgoing ray can be derived, and according to the interface intersection point and the slope of the upward ray, the intersection point coordinates with the free surface and the ray slope when the free surface is incident can be calculated. Similarly, the downward ray slope and the intersection point coordinates to the sea bottom interface can be calculated when the free surface is reflected.

[0044] In step 2, the processing of step 2 is implemented for each trace data of each common shot gather, and then the theoretical time-distance relationship of the first-order water-related multiple of each common shot gather can be drawn; the simulation of the second-order water-related multiple is similar to the first-order, and one upward wave and one downward wave simulation are added on the basis of the first-order; and the theoretical time-distance relationship of the water wave field on the common receiver gather can be obtained through data rearrangement.

[0045] In step 3, based on the calculated simulation of the OBN data direct wave, the first-order and the second-order water-related multiple, the direct wave, the first-order and the second-order water-related multiple in the actual data common receiver gather are identified based on Markov optimal decision under model constraint in each common receiver gather when the direct wave, the first-order and the second-order water-related multiple arrive.

[0046] Step 3 includes:

[0047] A3.1: Generate the feature attributes of each common receiver gather, including energy attributes, envelope attributes, and waveform similarity attributes;

[0048] A3.2: According to the time-distance relationship of the direct wave, the distribution area of the direct wave arrival time can be formed; the attributes related to the characteristics of the direct wave are generated; the cumulative reward value of the tracking position is evaluated according to the state transition probability and the instantaneous reward function, and the path with the maximum cumulative reward value is finally obtained; the direct wave tracking detection is performed in the constraint area from the common receiver gather with the smallest offset.

[0049] A3.3: Based on the time-distance relationship of first-order and second-order water-related multiples, the distribution areas of first-order and second-order water-related multiples can be formed; attributes related to the characteristics of direct waves are generated; the cumulative reward value of the tracking position is evaluated according to the state transition probability and the instantaneous reward function, and finally the path with the maximum cumulative reward value is obtained; starting from the gather with the smallest shot-receiver distance, direct wave tracking and detection are performed within the constrained area.

[0050] In step 3.1, since the energy of the direct wave decays slowly in the water, it is the first strong energy wave to reach the OBN node. The energy characteristics of the seismic traces before and after the direct wave should change significantly. A time window is set and divided into two parts. The arrival time of the strong energy wave field is determined by the energy ratio before and after the time window.

[0051] For any real signal f(t), there is a corresponding analytic signal. The envelope property can be obtained through Hilbert; the envelope of the data reflects the macroscopic changes of the waveform in the time domain.

[0052]

[0053]

[0054] Where f(t) is the time-domain signal, H represents the Hilbert transform of f(t), K represents the Cauchy principal value, and E(t) is the Hilbert transform envelope of f(t).

[0055] Suppose there is a time window containing N channels centered on the analysis point, and use the similarity coefficient r to measure the overall similarity of multiple channels within the time window;

[0056]

[0057] Where N represents the number of channels contained within the local window, x n The nth coordinate is represented by q, the local scanning tilt angle of the in-phase axis within the time window is represented by u, the local time window length is represented by M, and the time sampling interval is Δt. The similarity coefficient is used to measure the continuity of a waveform segment.

[0058] Markov decision process is a learning process that an agent takes action according to the environment to change its state and obtain rewards; it is generally composed of a five-tuple <S, A, T, P, r>; for the tracking of direct wave, wherein the state set S represents the current time part t of performing decision, and is used to describe the information of the selected direct wave travel time position; the action A represents that under the given model and time-distance relationship constraint, the direct wave position of the next trace is estimated according to the action probability π(a|s); the action generally has only two, one is to move one sampling point in the direction of travel time decrease, and the other is to move one sampling point in the direction of travel time increase, and the action probability π can be regarded as random; the state transition matrix P represents that according to the attribute information of the current state s, the possibility of transferring to the state s' can be achieved by the action A; according to the control of the travel time difference between the two adjacent traces, the overall distribution presents that the greater the difference from the theoretical travel time difference, the lower the transition probability; the reward function R represents the instantaneous reward value r obtained after the transition of A, which is constructed by the characteristic attributes:

[0059]

[0060] wherein M represents the total number of attributes used, w represents the weight of the attribute, and f represents the attribute value of the current position;

[0061] The time sequence T represents the entire evolution process; the entire state transition decision process is to find a decision sequence that can produce the maximum cumulative reward value;

[0062]

[0063] wherein γ represents a discount factor that decays over time, p is a state transition matrix, a is a current action, S represents an attribute used to describe the information of the selected water body wave field travel time position, and v n (S) is the cumulative reward value obtained at the current time step, v n+1 (S) is the cumulative reward value obtained at the next time step.

[0064] In step 3.2, according to the best decision principle, the direct wave path is tracked within the theoretical range of the direct wave travel time; starting from the trace set with the minimum offset, the attribute value at each time point is calculated within the constraint region, and the action probability is randomly moved up or down to another time point, if the expected cumulative reward value is greater, the time step is selected as the next state, and the process is repeated to track and detect the direct wave; the introduction of the intertrace travel time difference constraint between the adjacent two traces ensures the smoothness of the travel time as a whole.

[0065] In step 3.3, according to the best decision principle, the first-order and second-order water body related multiple wave paths are tracked within the first-order and second-order water body related multiple wave travel time theory range; the direct wave tracking detection is carried out in the constraint area from the minimum offset gather; the adjacent trace travel time difference constraint is introduced between the adjacent two traces to guarantee the smoothness of the travel time.

[0066] The water body body wave phenomenon identification method based on the wave propagation ray theory simulation in seawater in the application completes the water body body wave phenomenon identification based on the wave propagation ray theory simulation in seawater, and contains direct waves, first-order and second-order water body related multiples. The application simulates the wave field propagation path in the water body by using the trial shooting method, and realizes the tracking identification of the wave field based on the Markov best decision under the model constraint. The application has the advantages of simple technical operation, overcoming the problem that the water body wave field is difficult to identify, and helping to extract the wave field caused by the water layer alone, which is used for fine modeling of the seafloor surface form; that is, laying a foundation for subsequent establishment of a fine water body velocity model, and providing a reference for OBN seismic data processing and interpretation. Compared with the prior art, the technical scheme provided by the application has the advantages that (1) it provides complex OBN data processing and interpretation reference for researchers; (2) it helps to extract the wave field caused by the water layer alone, which is used for fine modeling of the seafloor surface form and inversion of the water body model; (3) it helps to extract the wave field caused by the water layer alone, which can highlight the deep reflection wave, and is convenient for subsequent multiple wave removal, extraction of related wave phenomena and migration imaging processing. BRIEF DESCRIPTION OF DRAWINGS

[0067] Figure 1 The flowchart of a specific embodiment of the water body body wave phenomenon identification method based on the wave propagation ray theory simulation in seawater in the application;

[0068] Figure 2 The water body model schematic diagram and the trial shooting method simulation first-order water body related multiple wave schematic diagram in a specific embodiment 1 of the application;

[0069] Figure 3 The calculation common receiver gather similarity coefficient attribute schematic diagram in a specific embodiment 1 of the application;

[0070] Figure 4 The first-order water body related multiple wave identification and tracking schematic diagram in actual OBN data in a specific embodiment 1 of the application;

[0071] Figure 5 The first-order water body related multiple wave identification and tracking travel time schematic diagram in actual OBN data in a specific embodiment 1 of the application.

[0072] Figure 6 The water body model schematic diagram and the trial shooting method simulation first-order water body related multiple wave schematic diagram in a specific embodiment 2 of the application;

[0073] Figure 7 Figure 2 is a schematic diagram of calculating the similarity coefficient attribute according to the fluctuating sea bottom water body model shot gather in a specific embodiment 2 of the present application;

[0074] Figure 8 Figure 3 is a schematic diagram of the first-order multiple wave travel time pickup interval obtained according to the water body model constraint in a specific embodiment 2 of the present application;

[0075] Figure 9 Figure 4 is a schematic diagram of the first-order water body related multiple wave identification and travel time tracking in the simulated shot gather in a specific embodiment 2 of the present application;

[0076] Figure 10 Figure 5 is a schematic diagram of the first-order water body related multiple wave simulation based on the observation system by the trial shooting method in a specific embodiment 3 of the present application;

[0077] Figure 11 Figure 6 is a schematic diagram of the first-order multiple wave travel time pickup interval obtained according to the water body model constraint in a specific embodiment 3 of the present application;

[0078] Figure 12 Figure 7 is a schematic diagram of the direct wave and the first-order water body related multiple wave travel time tracking in the actual OBN data in a specific embodiment 3 of the present application. DETAILED DESCRIPTION

[0079] It should be noted that the following detailed description is merely exemplary in nature and is intended to provide further description of the application. Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs.

[0080] It is also to be understood that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting. As used herein, the singular forms "a", "an" and "the" are intended to include the plural forms as well, unless the context clearly indicates otherwise. It will be further understood that the terms "comprises" and / or "comprising," when used in this specification, specify the presence of stated features, steps, operations, elements, and / or components, but do not preclude the presence or addition of one or more other features, steps, operations, elements, components, and / or groups thereof.

[0081] As shown in Figure 1 Figure 8 is a flowchart of the water body body wave phenomenon identification method based on the wave propagation ray theory simulation in seawater of the present application. The water body body wave phenomenon identification method based on the wave propagation ray theory simulation in seawater includes: Figure 1

[0082] Step 1: OBN common geophone gather reading, obtaining the trace header information of the P component data thereof, and obtaining the depth value of all OBN nodes on a single line.

[0083] ​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​Step 2: According to the OBN node depth value. The water body model interface is expressed in the form of an analytical function. Based on the initial seabed surface model, the theoretical first-order and second-order water body related multiple wave arrival times of each shot and each OBN node are calculated by trial shooting method.

[0084] Step 3: Based on the calculated simulated OBN data direct wave, first-order and second-order water body related multiple wave arrival times, the actual data direct wave and first-order and second-order water body related multiple waves in each common receiver gather are identified based on Markov optimal decision under model constraint. The water body wave field travel time acquisition result is output.

[0085] Further, step 1 includes the following detailed steps:

[0086] A1.1: OBN shot line reading, obtaining trace header information, including line number, shot coordinates, receiver coordinates, elevation, offset, sampling rate, sampling time, etc.

[0087] A1.2: Record the depth values of all OBN nodes on a single line.

[0088] Further, for each single shot, step 2 simulates the propagation of the ray, and step 2 includes the following detailed steps:

[0089] A2.1: According to the OBN node depth value. The interface of the water body model is expressed in the form of an analytical function.

[0090] A2.2: Given the starting point position of the ray and assuming a ray emission direction, the ray will reflect according to the geometric properties after reaching the interface, thus tracing a ray, calculating its intersection point when reflecting on the interface and its final landing point on the interface.

[0091] A2.3: Adjust the initial emission direction of the ray at certain emission angle intervals to find the left and right two initial emission angles of the ray closest to the target receiver point coordinates.

[0092] A2.4: Further narrow the angle range of the target ray by bisection method, so that the calculated theoretical emission ray intersection point coordinates on the interface approach the coordinates of the target receiver point. Record the intersection point when the ray passes through the interface.

[0093] A2.5: Assuming a constant speed background, calculate the travel time of the wave reaching the receiver point according to the length of the ray path.

[0094] A2.6: Obtain the theoretical time-distance relationship of the first-order and second-order water body wave field on the common receiver gather by data rearrangement.

[0095] Further, step 3 is implemented for each common receiver gather data on the line, and step 3 for waveform tracking includes the following detailed steps:

[0096] A3.1: Generating the characteristic attributes of each common receiver gather, including energy attributes, envelope attributes, and waveform similarity attributes.

[0097] A3.2: According to the direct wave time-distance relationship, a direct wave arrival time distribution area can be formed, and attributes related to the characteristics of the direct wave are generated. According to the state transition probability and the instantaneous reward function, the cumulative reward value of the tracking position is evaluated, and finally the path with the maximum cumulative reward value is obtained. Starting from the gather with the smallest offset, direct wave tracking detection is performed within the constraint area.

[0098] A3.3: According to the first-order and second-order water-related multiple wave time-distance relationship, a first-order and second-order water-related multiple wave distribution area can be formed, and attributes related to the characteristics of the direct wave are generated. According to the state transition probability and the instantaneous reward function, the cumulative reward value of the tracking position is evaluated, and finally the path with the maximum cumulative reward value is obtained. Starting from the gather with the smallest offset, direct wave tracking detection is performed within the constraint area.

[0099] The following are several specific embodiments of the application

[0100] Embodiment 1

[0101] In a specific embodiment 1 of the application, the water body wave phenomenon recognition method based on the wave propagation ray theory simulation in seawater mainly involves the recognition of direct waves, first-order and second-order water-related multiple waves in OBN observation data. The shot gather data used in embodiment 1 is shown in Figure 2 , which is forward modeling data generated based on a horizontal seabed model. The model size is 10000m*600m, 121 shots are excited, and 91 OBN nodes are simulated. The shot spacing is 50m, the receiver spacing is 100m, the maximum offset is 9000m, and the water depth is 250m.

[0102] Before execution, the OBN gather data file directory, the output file directory, the gather basic information, and the prior parameters designed by the processing personnel need to be provided. The water body wave phenomenon recognition method based on the wave propagation ray theory simulation in seawater includes:

[0103] The first step is 1: input of OBN common receiver gather, and acquisition of trace header information.

[0104] Step 1 specifically includes the following steps:

[0105] A1.1: OBN common receiver gather reading, acquisition of P component data trace header information, including line number, shot coordinates, receiver coordinates, elevation, offset, sampling rate, sampling time, etc.

[0106] A1.2: Record and arrange the depth values of all OBN nodes on a single line.

[0107] The ray tracing in step 2 is performed by the trial method, and the first-order water body related multiple wave calculation is taken as an example, and step 2 specifically comprises the following steps:

[0108] A2.1: Description of the interface

[0109] The upper interface of the water body model is regarded as a horizontal free surface, the lower interface is expressed by an analytical function, and the internal model is regarded as a uniform continuous medium. For a two-dimensional geological body, the coordinates and depths of each OBN node are known, and the depth value of any position of the interface is obtained by cubic spline function interpolation and the like.

[0110] The model interpolation interface in example 1 is shown in Figure 2 .

[0111] According to the local slope of the interface, the slope of the upgoing emergent ray can be derived, and according to the intersection point of the interface and the slope of the upgoing ray, the intersection coordinate of the free surface and the ray slope when the free surface is incident can be calculated, and similarly the downgoing ray slope and the intersection coordinate of the sea bottom interface when the free surface is reflected can be calculated.

[0112] A2.3: Adjust the initial emergent direction of the ray according to a certain emergent angle interval, and obtain the left and right two initial emergent angles of the ray located at the target receiving point coordinate.

[0113] A2.4: Further narrow the angle range of the target ray by the dichotomy method, so that the calculated theoretical emergent ray reaches the interface intersection coordinate close to the coordinate of the target receiving point. Record the ray path from the shot point to the receiving point and the intersection point when passing through the interface.

[0114] A2.5: Assuming a constant speed background, the travel time of the wave reaching the receiving point is calculated according to the length of the ray path.

[0115] The first-order multiple wave ray path simulated according to the trial method in example 1 is shown in Figure 2 .

[0116] A3.1: Calculate the common receiver point gather characteristic attribute, including energy attribute, envelope attribute, and waveform similarity attribute.

[0117] For any real signal f(t), there is a corresponding analytical signal The envelope attribute can be obtained by Hilbert, and the envelope of the data reflects the macroscopic change of the waveform in the time domain.

[0118]

[0119]

[0120] Where H represents the Hilbert transform of f(t), and K represents the Cauchy principal value.

[0121] Suppose a time window containing N traces is centered at an analysis point, the similarity coefficient r is used to measure the similarity of the whole traces in the time window;

[0122]

[0123] Where N represents the number of traces contained in the local window, x n represents the n-th trace coordinate, q represents the local scanning dip angle of the event in the time window, u represents the amplitude, M represents the length of the local time window, and Δt is the time sampling interval; the similarity coefficient is used to measure the continuity of a section of waveforms; Figure 3 is a schematic diagram of calculating the similarity coefficient attribute of the common receiver point gather in a specific embodiment of the application.

[0124] A Markov decision state transition decision process is constructed to find a decision sequence that can generate the maximum cumulative reward value.

[0125]

[0126] Where γ represents a discount factor that decays over time, p is a state transition matrix, a is the current action, S represents an attribute of information used to describe the selected water body wave field travel time position, v n (S) is the cumulative reward value obtained at the current time step, v n+1 (S) is the cumulative reward value obtained at the next time step.

[0127] A3.2: According to the best decision principle, the direct wave path is tracked within the range of direct wave travel time theory. Starting from the gather with the smallest offset, the attribute value at each time point is calculated within the constraint region, and the action probability is randomly moved to another time point upward or downward, and if the expected cumulative reward value is greater, the time step is selected as the next state, and the process is repeated to perform direct wave tracking detection; the introduction of the intertrace travel time difference constraint between adjacent two traces ensures the smoothness of the travel time obtained.

[0128] Figure 4 The first-order multiple wave travel time picking interval obtained according to the water body model constraint in Example 1 is shown.

[0129] A3.3: According to the best decision principle, the first-order and second-order water body related multiple wave paths are tracked within the range of the first-order and second-order water body related multiple wave travel time theory. Starting from the gather with the smallest offset, direct wave tracking detection is performed within the constraint region; the intertrace travel time difference constraint is introduced between adjacent two traces to ensure the smoothness of the travel time obtained.

[0130] Figure 5 The first-order water body related multiple wave identification and travel time tracking schematic diagram in the simulation shot gather of Example 1 is shown.

[0131] The model data used in Example 2 is as follows: Figure 6 The image shows forward modeling data generated based on a model of an undulating seabed. The model size is 10000m * 801m, with 121 shots fired and 91 simulated OBN nodes. The shot spacing is 50m, the receiver spacing is 100m, the maximum shot-receiver distance is 9000m, the time sampling rate is 0.4ms, and the average water depth is 300m. The travel times of first-order multiple waves in the water wave field were tracked using the aforementioned wave phenomenon identification process.

[0132] Figure 6 This is a schematic diagram of the water body model in Example 2 and a schematic diagram of the first-order water body related multiple waves simulated by the test firing method;

[0133] Figure 7 This is a schematic diagram illustrating the calculation of similarity coefficient attributes based on the shot gather of the undulating seabed water body model in Example 2 of this embodiment;

[0134] Figure 8 This is a schematic diagram of the first-order multiple wave travel time picking interval obtained based on the water body model constraints in Embodiment 2.

[0135] Figure 9 This is a schematic diagram of the identification and travel time tracking of first-order water-related multiples in the simulated shot concentration in Example 2.

[0136] The shot collection data used in Example 3 are as follows: Figure 2 The data shown is actual OBN data. The target water body model is a shallow water model with a depth between 60 and 80 meters. The OBN spacing is approximately 50 meters, the maximum offset is approximately 8000 meters, and the time sampling rate is 0.1 ms. The travel time of the first-order multiple waves in the water wave field is tracked using the above wave phenomenon identification process.

[0137] Figure 10 This is a schematic diagram of simulating first-order water-related multiples using a test firing method based on an observation system in Embodiment 3 of this example;

[0138] Figure 11 This is a schematic diagram of the first-order multiple wave travel time picking interval obtained based on the water body model constraints in specific embodiment 3;

[0139] Figure 12 This is a schematic diagram of the actual OBN data direct wave and first-order water body related multiple wave travel time tracking in Example 3.

[0140] Finally, it should be noted that the above only describes the preferred embodiments of the present application and is not intended to limit the present application. Although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art will appreciate that modifications can be made to the technical solutions described in the foregoing embodiments, or some of the technical features thereof can be replaced equivalently, without departing from the spirit and principle of the present application. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present application should be included in the protection scope of the present application.

[0141] All that is not described in the specification is known to those skilled in the art.

Claims

1. A method for identifying water body wave phenomena based on wave propagation ray theory simulation in seawater, characterized in that, The water body wave phenomenon identification method based on the ray theory simulation of wave propagation in seawater comprises: Step 1, obtaining the depth values of all OBN nodes on a single line; Step 2, simulating the water body propagation wave field by trial shooting; Step 3, realizing the tracking identification of the wave field based on the Markov optimal decision: Step 4, outputting the water body wave field travel time acquisition result; Step 2 comprises: A2.1: according to the OBN node depth values, expressing the interface of the water body model in the form of an analytical function; A2.2: given the ray starting point position and assuming a ray emission direction, the ray will reflect according to the geometric properties after reaching the interface, thus tracking a ray, calculating the intersection point when the ray reflects on the interface and the final landing point on the interface; A2.3: adjusting the initial emission direction of the ray at certain emission angle intervals to obtain the left and right two initial emission angles of the ray closest to the target receiving point coordinates; A2.4: further narrowing the angle range of the target ray by the dichotomy method, so that the calculated theoretical emission ray reaches the interface intersection point coordinates close to the coordinates of the target receiving point, and the intersection point when the ray passes through the interface is recorded; A2.5: assuming a constant speed background, calculating the wave travel time to the receiving point according to the length of the ray path; A2.6: obtaining the theoretical time-distance relationship of the first-order and second-order water body waves on the common receiver point gather by data rearrangement; Step 3 comprises: A3.1: generating the characteristic attributes of each common receiver point gather, including energy attributes, envelope attributes, and waveform similarity attributes; A3.2: according to the direct wave time-distance relationship, a direct wave arrival time distribution area can be formed; the attributes related to the direct wave characteristics are generated; the cumulative reward value of the tracking position is evaluated according to the state transition probability and the instantaneous reward function, and finally the maximum path of the cumulative reward value is obtained; starting from the gather with the smallest offset, the direct wave tracking detection is performed within the constraint area; A3.3: according to the time-distance relationship of the first-order and second-order water body related multiples, the first-order and second-order water body related multiple distribution areas can be formed; the attributes related to the direct wave characteristics are generated; the cumulative reward value of the tracking position is evaluated according to the state transition probability and the instantaneous reward function, and finally the maximum path of the cumulative reward value is obtained; starting from the gather with the smallest offset, the direct wave tracking detection is performed within the constraint area.

2. The method according to claim 1, wherein, In step 1, OBN common receiver point gather is read to obtain the trace header information of P component data and the depth values of all OBN nodes on a single line.

3. The method according to claim 2, wherein, Step 1 comprises: A1.1: OBN shot line reading is performed to obtain the trace header information, including line number, shot coordinates, receiver point coordinates, elevation, offset, sampling rate, and sampling time; A1.2: recording the depth values of all OBN nodes on a single line.

4. The method according to claim 1, wherein, In step 2, the seabed interface of the water body model is expressed in the form of an analytical function according to the OBN node depth values; based on the initial seabed surface model, the theoretical first-order and second-order water body related multiples arriving time of each shot and each OBN node is calculated by trial shooting.

5. The method according to claim 1, wherein, In step 2.1, for a two-dimensional geological body, the coordinates and depth of each OBN node are known, and the depth value at any position of the interface is obtained by linear interpolation, polynomial interpolation, and cubic spline function interpolation. Assuming that there are n+1 OBN nodes, which are divided into n intervals, the depth of a point x∈[x k ,x k+1 ] in a certain interval is described by a cubic spline function. y = a i + b i x + c i x 2 + d i x 3 where y is the depth of the interpolation point, a i , b i , c i , d i are the constant term, the linear term coefficient, the quadratic term coefficient, the cubic term coefficient of the interpolation function, respectively. All points in a cubic spline function must satisfy the interpolation condition: f i (x i )=y i (i = 0, 1, ..., n), and secondly, the curve needs to be smooth, that is, f at the endpoints. i (x i ),f i '(x i ),f i (x) i ) Continuous; therefore, all internal endpoints must satisfy the cubic equations of the left and right segments; the equation for the continuity of the first and second derivatives of the n-1 internal points: f' i (x i )=f' i+1 (x i ),f” i (x i )=f” i+1 (x i Assume natural boundary conditions, i.e., f”0(x0)=f” n (x n If ) = 0, then a total of 4n equations are formed, and the linear equation is derived through derivation: where h i = x i+1 - x i , m i = f i (x i ), h i denotes the step size, m i is the second derivative at the node; a system of linear equations is constructed: Solve with Jacobi iteration method, then return m i to the coefficients of cubic spline function; There is the following formula: For a three-dimensional geological body, the interface is described in the form of triangular patches or B-spline surfaces.

6. The method according to claim 1, wherein, In step 2, the step 2 processing is implemented for each trace of each common shot gather, and then the theoretical time-distance relationship of the first-order water-related multiple of each common shot gather can be drawn; the simulation of the second-order water-related multiple is similar to the first-order, and one upgoing wave and one downgoing wave simulation are added on the basis of the first-order; and the theoretical time-distance relationship of the water wave field on the common receiver gather can be obtained through data rearrangement.

7. The method according to claim 1, wherein, In step 3, based on the arrival time of the simulated OBN data direct wave, the first-order and the second-order water-related multiple, the direct wave and the first-order and the second-order water-related multiple in the actual data common receiver gather are identified in each common receiver gather based on the Markov best decision under the model constraint.

8. The method according to claim 7, wherein, In step 3.1, since the direct wave energy decays slowly in the water body, it is the strong energy wave phenomenon that first arrives at the OBN node, and the energy characteristics of the seismic traces before and after the direct wave should change greatly; a time window is set, the time window is divided into two parts before and after, and the arrival time of the strong energy wave field is determined by the energy ratio before and after the time window; For any time-domain real signal f(t), there is a corresponding analytic signal The envelope property can be obtained by Hilbert, and the envelope of data reflects the macroscopic changes of the waveform in the time domain; Where H represents the Hilbert transform of f(t), K represents the Cauchy principal value, and E(t) is the Hilbert transform envelope of f(t); Suppose a time window containing N traces centered on the analysis point, the similarity coefficient r is used to measure the similarity of the whole in the time window; where N represents the number of traces contained in the local window, x n represents the coordinate of the nth trace, q represents the local scan angle of the event in the time window, U n represents the amplitude, M represents the length of the local time window, and Δt is the time sampling interval; the similarity coefficient is used to measure the continuity of a section of waveforms; The Markov decision process is a learning process in which an agent changes its state to obtain a reward according to the environment; It is generally composed of five tuples <S, A, T, P, r>; for the tracking of direct waves, the state set S represents the current time part t of the decision execution, which is used to describe the information of the selected direct wave travel time position; the action A represents the estimation of the direct wave position of the next trace under the given model and time-distance relationship constraint according to the action probability π(a|s), and the action generally has only two, one is to move one sampling point in the direction of decreasing travel time, and the other is to move one sampling point in the direction of increasing travel time, and the action probability π can be regarded as random; the state transition matrix P represents the possibility of transferring to the state s' according to the attribute information of the current state s when the action A is implemented, and the overall distribution presents that the greater the difference between the theoretical travel time difference and the travel time difference, the lower the transition probability; the reward function R represents the instantaneous reward value r obtained after the transition of A, which is constructed by the characteristic attributes: where M t denotes the total number of attributes used, w m denotes the weight of the attribute, f m (s) denotes the attribute value of the current position, r(s) is the current point instantaneous reward value; The time sequence T represents the entire evolution process; The entire state transition decision process is to find a decision sequence that can produce the maximum cumulative reward value. where γ represents a discount factor that decays over time, a is the current action, S represents an attribute of information describing the selected water body wave field travel time position, v n (S) is the cumulative reward value obtained at the current time step, v n+1 (S) is the cumulative reward value obtained at the next time step.

9. The method according to claim 7, wherein, In step 3.2, according to the best decision principle, the direct wave path is tracked within the theoretical range of the direct wave travel time; starting from the shot gather with the smallest offset, the attribute value at each time point is calculated within the constraint region, and the action probability is randomly moved to another time point, if the expected cumulative reward value is greater, the time step is selected as the next state, and the process is repeated to perform direct wave tracking detection; the introduction of the travel time difference constraint between adjacent traces ensures the smoothness of the travel time.

10. The method according to claim 7, wherein, In step 3.3, according to the best decision principle, the first-order and second-order water-related multiple path is tracked within the first-order and second-order water-related multiple travel time theory range; the direct wave tracking detection is carried out in the constraint area from the minimum offset gather; the adjacent trace travel time difference constraint is introduced between the adjacent two traces to ensure the smoothness of the travel time as a whole.

Images