Water body wave phenomenon identification method based on seawater medium wave propagation ray theoretical simulation

Abstract

The invention provides a water body wave phenomenon identification method based on seawater medium wave propagation ray theory simulation. The water body wave phenomenon identification method based on seawater medium wave propagation ray theory simulation comprises the following steps: step 1, obtaining depth values of all OBN nodes on a single line; step 2, simulating a water body propagation wave field through a test shooting method; step 3, realizing tracking identification of a wave field based on a Markov optimal decision; and step 4, outputting a water body wave field travel time acquisition result. The water body wave phenomenon identification method based on seawater medium wave propagation ray theoretical simulation is simple and convenient in technical operation, overcomes the problem that a water body wave field is difficult to identify, is beneficial to independently extracting a wave field caused by a water layer, and is used for fine modeling of a seabed surface form; a foundation is laid for subsequently establishing a fine water body speed model, and reference is provided for OBN seismic data processing and interpretation work.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic data processing for seismic exploration, and in particular to a method for identifying water body wave phenomena based on theoretical simulation of wave propagation rays in seawater. Background Art

[0002] Due to the existence of two strong reflection interfaces of water bodies, namely the seawater-air interface and the seawater-seabed interface, the water body wave field appears as strong energy in the OBN record. Direct waves, multi-order water body related multiple waves are intertwined with the effective reflection wave field from below the water body interface. At far offsets, multiple reflection and refraction strong energy interference from multiple oscillations in the water body also appears, which brings great difficulties to subsequent processing. The complex wave phenomenon caused by complex shallow water layers is one of the important interference waves. If the water body wave phenomenon is identified, it is conducive to the refined tomographic inversion of the water body model, and it is also convenient to extract the wave field caused by the water layer separately in the future, highlighting the deep reflection wave. However, there has been little research on the identification of water body wave phenomena.

[0003] In the Chinese patent application with application number: CN202010535639.5, a method for simulating numerical data of seismic wave fields of fluctuating sea surface based on wave spectrum is involved. The first-order velocity-stress acoustic wave equation group expression is obtained by derivation and order reduction processing; the acoustic wave equations of each grid point and half grid point in time and space are calculated by finite difference according to the calculation rules of staggered grids; forward simulation of offshore data is performed; the spatial domain model of the fluctuating sea surface under the condition of PM spectrum is obtained; the fluctuating sea surface values ​​obtained are applied to finite difference numerical simulation, and the upper boundary values ​​of the above-mentioned model calculation area are replaced by sea surface values, so as to obtain the simulated data under the fluctuating sea surface. The numerical simulation method equation provided by the invention is simple to use, efficient, and the accuracy can be improved by certain means; the construction cost of three-dimensional seismic data acquisition is low, and the sinking depth of the source and the detector is moderate, which can effectively avoid the occurrence of false phase axis phenomenon.

[0004] In the Chinese patent application with application number: CN201811292058.2, a method for quickly establishing a three-dimensional near-bottom seabed velocity model in shallow sea areas is involved, including: extracting first arrival information and extracting common center point gathers; calculating the slope of the first arrival travel time in the common center point gathers; judging whether the shot point is located in the seawater layer, if so, correcting the shot point to the seabed; calculating the velocity and depth of the inflection point corresponding to each offset; linearly interpolating the velocity at each depth to obtain the one-dimensional velocity model of the common center point; after processing all the common center point domain gathers, using the linear interpolation method, using the one-dimensional velocity model of all the common center points, calculate the velocity values ​​of the remaining grid points, and output the velocity model. The method for quickly establishing a three-dimensional near-bottom seabed velocity model in shallow sea areas overcomes the influence of the seawater layer on ray propagation, eliminates the influence of multiple excitation and reception methods on the inversion results, and avoids the problem of conventional methods requiring multiple ray tracing forward simulations.

[0005] In the Chinese patent application with application number: CN201910816102.3, a method, device and storage medium for collecting ocean data are involved, which are used to solve the technical problem of how to optimize the OBN observation system in the related technology so that the ocean collection data meets the requirements of oil and gas exploration and development. The method for collecting ocean data includes: constructing a three-dimensional block geological model according to the measured geological data of the work area to be collected; simulating the layout of multiple OBN observation systems in the three-dimensional block geological model; using the shot point as the source point to simulate Gaussian ray irradiation to the target layer in the three-dimensional block geological model; obtaining the facet offset energy analysis result of the target layer according to the illumination energy of the Gaussian ray irradiation received by the detector in the simulated OBN observation system and the illumination energy reflected by the target layer; selecting the target OBN observation system from multiple OBN observation systems according to the facet offset energy analysis result; and placing the 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 with application number: CN202210592015.6, a two-dimensional sea surface ghost wave water body imaging measurement method, system, terminal and flow measurement equipment are involved. Obtain the frequency domain forward transformation data of the seismic data; for any frequency component and ray parameter, calculate the delay time between the virtual reflection and the primary reflection, and solve to obtain the frequency domain Radon transform result of the first downlink wave field; obtain the average delay time and average amplitude difference coefficient of the virtual reflection and the primary reflection corresponding to the ray parameter and calculate again, and continuously alternate and iterate until the convergence condition is met to obtain the frequency domain Radon transform result corresponding to the final virtual reflection, and then perform frequency domain Radon inverse transform and time domain Fourier inverse transform to obtain the final virtual reflection prediction data. The invention can perform full waveform inversion on useful information such as the virtual reflection signal received by ADCP to obtain the property information of the surface seawater and realize accurate imaging of the water body.

[0007] The above existing technologies are all significantly different from the present invention and fail to solve the technical problem we want to solve. For this reason, we have invented a new method for identifying water body wave phenomena based on theoretical simulation of wave propagation rays in seawater. Summary of the invention

[0008] The purpose of the present invention is to provide a method for identifying body wave phenomena in water based on the ray theory simulation of wave propagation in seawater, which uses a test shooting method to simulate the propagation path of wave fields in water bodies and realizes tracking and identification of body wave fields in water based on Markov optimal decision under model constraints.

[0009] The object of the present invention can be achieved by the following technical measures: a method for identifying water body wave phenomena based on theoretical simulation of wave propagation rays in seawater, the method for identifying water body wave phenomena based on theoretical simulation of wave propagation rays in seawater comprises:

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

[0011] Step 2, simulating the wave field propagating in water by the test shooting method;

[0012] Step 3: Track and identify the wave field based on Markov optimal decision;

[0013] Step 4: Output the water body wave field travel time to obtain the result.

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

[0015] In step 1, the OBN common detection point gather is read to obtain the header information of its P component data and the depth values ​​of all OBN nodes on the single line.

[0016] Step 1 includes:

[0017] A1.1: Read the OBN shot line and obtain the trace header information, including the line number, shot coordinates, detection point coordinates, elevation, offset, sampling rate, and sampling time;

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

[0019] 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 value; based on the initial seabed model, the theoretical first-order and second-order water-related multiple wave arrival times of each shot and each OBN node are calculated by the trial firing method.

[0020] Step 2 includes:

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

[0022] A2.2: Given the starting point of the ray and an assumed direction of ray emission, the ray will be reflected according to the geometric properties after reaching the interface. In this way, a ray is traced, and its intersection point when reflected on the interface and its final landing point on the interface are calculated;

[0023] A2.3: Adjust the initial emission direction of the ray at a certain emission angle interval, and find the initial emission angles of the two rays closest to the target receiving point;

[0024] A2.4: Use the dichotomy method to further reduce the angle range of the target ray so that the coordinates of the intersection point of the theoretical outgoing ray reaching the interface are close to the coordinates of the target receiving point, and record the intersection point when the ray passes through the interface;

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

[0026] A2.6: The theoretical time-distance relationship of the first-order and second-order water wave fields on the common detection point gathers is obtained 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 obtained by linear interpolation, polynomial interpolation, and cubic spline function interpolation; assuming that there are n+1 OBN nodes, divided into n intervals, and a cubic spline function is used to describe a point x∈[x k ,x k+1 Depth of:

[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 They 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, the endpoints f i (x i ),f i '(x i ),f i ”(xi ) is continuous; therefore, all internal endpoints need to satisfy the cubic equations of the two left and right segments; the equations for the first and second order derivatives of n-1 internal points are continuous: f i '(x i ) = f i ' +1 (x i ),f i ”(x i ) = f i ” +1 (x i ); Assume natural boundary conditions, that is, f0”(x0)=f n ”(x n )=0, then a total of 4n equations are formed, and the linear equation is obtained by deduction:

[0031]

[0032] where h i =x i+1 -x i ,m i =f i ”(x i ), h i Indicates the step length, m i is the second-order derivative at the node. Construct a linear system of equations

[0033]

[0034] Use Jacobi iteration method to solve, and then substitute m back into the coefficients of the cubic spline function; the formula is as follows:

[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 ray emission direction, the down-going wave process and the up-going wave process of the ray propagation path are calculated; given the down-going ray slope k1 and intercept b1, the coordinates of the intersection of the down-going ray and the seabed interface of the water body model are calculated by the simultaneous formula;

[0038]

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

[0040] Given the coordinates of the incident point of the ray and the coordinates of the interface intersection point, according to the slope k1 of the incident ray, the direction vector of the incident ray is It is also necessary to find the direction vector of the outgoing ray after reflection According to the law of reflection, the incident angle is equal to the outgoing angle. The direction vector of the incident ray and the direction vector of the reflected ray are symmetric with respect to the normal vector. To find the direction vector of the outgoing ray:

[0041]

[0042] in is the local slope of the interface at the intersection, the coordinates of the ray incident point are (x1, y1), the coordinates of the interface intersection 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. According to the intersection point of the interface and the slope of the upward ray, the coordinates of the intersection with the free surface and the slope of the ray when it enters the free surface can be calculated. Similarly, the slope of the downward ray when reflected from the free surface and the coordinates of the intersection when it reaches the seabed interface can be calculated.

[0044] In step 2, the processing of step 2 is implemented for each data of each common shot gather, and the theoretical time-distance relationship of the first-order water-related multiple waves of each common shot gather can be plotted; the simulation of the second-order water-related multiple waves is similar to the first-order, and one upgoing wave and one downgoing wave are added to the first-order simulation; and the theoretical time-distance relationship of the water wave field on the common detection point gather can be obtained by data rearrangement.

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

[0046] Step 3 includes:

[0047] A3.1: Generate characteristic attributes of each common detection point gather, including energy attributes, envelope attributes, and waveform similarity attributes;

[0048] A3.2: The distribution area of ​​the direct wave arrival time can be formed according to the time-distance relationship of the direct wave; 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 path with the maximum cumulative reward value is obtained; starting from the gather with the smallest shot offset, the direct wave tracking detection is carried out in the constrained area;

[0049] A3.3: The distribution area of ​​the first-order and second-order water-related multiple waves can be formed according to the time-distance relationship of the first-order and second-order water-related multiple waves; 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 path with the maximum cumulative reward value is obtained; starting from the data set with the smallest shot offset, the direct wave tracking detection is carried out in the constrained area.

[0050] In step 3.1, since the energy of the direct wave decays slowly in the water body, it is the first strong energy wave to arrive at the OBN node. The energy characteristics of the seismic traces before and after the direct wave should change greatly. Set a time window, divide the time window into two parts, and determine the arrival time of the strong energy wave field by the energy ratio between the front and rear parts of the time window.

[0051] For any real signal f(t) there is a corresponding analytic signal The envelope property can be obtained through Hilbert, and 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] Assume a time window with N channels centered at the analysis point, and use the similarity coefficient r to measure the overall similarity of multiple channels in the time window;

[0056]

[0057] Where N represents the number of channels contained in the local window, x n represents the coordinate of the nth track, q represents the local scanning inclination of the event axis 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 waveform;

[0058] The Markov decision process is a learning process in which an agent takes actions according to the environment to change its state and obtain rewards; 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 execution decision, which is used to describe the information of the selected direct wave travel time position; action A represents the estimation of the next direct wave position according to the action probability π(a|s) under the given model and time-distance relationship constraints. There are generally only two actions, one is to move a sampling point in the direction of decreasing travel time, and the other is to move a sampling point in the direction of increasing travel time. The action probability π can be regarded as random; the state transfer matrix P represents the possibility of transferring to state s' according to the attribute information of the current state s when implementing action A. According to the theoretical travel time difference control of two adjacent paths, the overall distribution shows that the greater the difference from the theoretical travel time difference, the lower the transfer probability; the reward function R represents the instantaneous reward value r obtained after the A transfer is realized, which is constructed by the characteristic attributes:

[0059]

[0060] Where 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 series 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] where γ represents the discount factor that decays over time, p is the state transfer matrix, a is the current action, S represents the attribute used to describe the information of the selected water 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 in 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 track set with the smallest shot offset, the attribute value of each time point is calculated in the constraint area, and it is randomly moved up or down to another time point according to the action probability. If the expected cumulative reward value is larger, the time step is selected as the next state, and this process is repeatedly executed to perform direct wave tracking detection; the introduction of the adjacent track travel time difference constraint between two adjacent tracks ensures the overall smoothness of the travel time.

[0065] In step 3.3, according to the best decision principle, the paths of the first-order and second-order water-related multiple waves are tracked within the theoretical range of the first-order and second-order water-related multiple wave travel times; starting from the track set with the smallest shot offset, direct wave tracking detection is performed within the constraint area; and the adjacent track travel time difference constraint is introduced between two adjacent tracks to ensure the overall smoothness of the acquired travel time.

[0066] The method for identifying water body wave phenomena based on the theoretical simulation of wave propagation rays in seawater in the present invention completes the identification of water body wave phenomena based on the theoretical simulation of wave propagation rays in seawater, including direct waves, first-order and second-order water body related multiple waves. The present invention uses the test shot method to simulate the propagation path of the wave field in the water body, and realizes the tracking and identification of the wave field based on the Markov optimal decision under the model constraint. The technology of the present invention is simple to operate, overcomes the problem that the water body wave field is difficult to identify, and is helpful to extract the wave field caused by the water layer separately for the refined modeling of the seabed surface morphology; that is, it lays the foundation for the subsequent establishment of a refined water body velocity model, and also provides a reference for the processing and interpretation of OBN seismic data. Compared with the prior art, the technical solution provided by the present invention has the following advantages: (1) It provides researchers with reference for the processing and interpretation of complex OBN data. (2) It is helpful to extract the wave field caused by the water layer separately for the refined modeling of the seabed surface morphology and the inversion of the water body model. (3) It is helpful to extract the wave field caused by the water layer separately, which can highlight the deep reflection wave, and is also convenient for the subsequent removal of multiple waves, extraction of related wave phenomena for offset imaging and other processing. BRIEF DESCRIPTION OF THE DRAWINGS

[0067] Figure 1 It is a flow chart of a specific embodiment of the method for identifying water body wave phenomena based on theoretical simulation of wave propagation rays in seawater of the present invention;

[0068] Figure 2 A schematic diagram of a water body model in a specific embodiment 1 of the present invention and a schematic diagram of a first-order water body related multiple wave simulation by a test firing method;

[0069] Figure 3 It is a schematic diagram of calculating the similarity coefficient attribute of the common detection point gather in a specific embodiment 1 of the present invention;

[0070] Figure 4 It is a schematic diagram of identifying and tracking first-order water-related multiple waves in actual OBN data in a specific embodiment 1 of the present invention;

[0071] Figure 5 It is a schematic diagram of identifying and tracking the travel time of first-order water-related multiple waves in actual OBN data in a specific embodiment 1 of the present invention.

[0072] Figure 6 A schematic diagram of a water body model in a specific embodiment 2 of the present invention and a schematic diagram of a first-order water body related multiple wave simulation by a test firing method;

[0073] Figure 7 It is a schematic diagram of calculating similarity coefficient attributes according to a shot gather of an undulating seabed water body model in a specific embodiment 2 of the present invention;

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

[0075] Fig. 9 It is a schematic diagram of the identification and travel time tracking of first-order water-related multiple waves in a simulated shot gather in a specific embodiment 2 of the present invention;

[0076] Fig.10 It is a schematic diagram of simulating first-order water-related multiple waves by a test-firing method based on an observation system in a specific embodiment 3 of the present invention;

[0077] Fig.11 It is a schematic diagram of the first-order multiple wave travel time picking interval obtained according to the water body model constraint in a specific embodiment 3 of the present invention;

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

[0079] It should be noted that the following detailed descriptions are exemplary and are intended to provide further explanation of the present invention. Unless otherwise specified, all technical and scientific terms used herein have the same meanings as those commonly understood by those skilled in the art to which the present invention belongs.

[0080] It should be noted that the terms used herein are only for describing specific embodiments and are not intended to limit exemplary embodiments according to the present invention. As used herein, unless the context clearly indicates otherwise, the singular form is also intended to include the plural form. In addition, it should be understood that when the terms "comprising" and / or "including" are used in this specification, it indicates the presence of features, steps, operations and / or combinations thereof.

[0081] like Figure 1 As shown, Figure 1 The flowchart of the method for identifying water body wave phenomena based on the theoretical simulation of wave propagation rays in seawater of the present invention. The method for identifying water body wave phenomena based on the theoretical simulation of wave propagation rays in seawater comprises:

[0082] Step 1: Read the OBN common detection point gather, obtain the header information of its P component data, and obtain the depth values ​​of all OBN nodes on the single line.

[0083] Step 2: According to the OBN node depth value, the seabed interface of the water body model is expressed in the form of an analytical function. Based on the initial seabed 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 the trial firing method.

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

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

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

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

[0088] Furthermore, the process of simulating propagation rays in step 2 is performed on each single shot. Step 2 includes the following detailed steps:

[0089] A2.1: According to the OBN node depth value, express the interface of the water model in the form of an analytical function.

[0090] A2.2: Given the starting point of the ray and an assumed direction of ray emission, the ray will be reflected according to the geometric properties after reaching the interface. In this way, a ray is traced, and its intersection point when reflected on the interface and its final landing point on the interface are calculated.

[0091] A2.3: Adjust the initial emission direction of the ray at a certain emission angle interval, and calculate the initial emission angles of the two left and right rays closest to the target receiving point coordinates.

[0092] A2.4: Use the dichotomy method to further reduce the angle range of the target ray so that the coordinates of the intersection point of the theoretical outgoing ray reaching the interface are close to the coordinates of the target receiving point. Record the intersection point when the ray passes through the interface.

[0093] A2.5: Assuming a constant velocity background, calculate the travel time of the wave to the receiving point based on the length of the ray path.

[0094] A2.6: The theoretical time-distance relationship of the first-order and second-order water wave fields on the common detection point gathers is obtained by data rearrangement.

[0095] The common detection point gather data on each line are further aligned to implement the processing of step 3. Step 3 for waveform tracking includes the following detailed steps:

[0096] A3.1: Generate characteristic attributes for each common detection point gather, including energy attributes, envelope attributes, and waveform similarity attributes.

[0097] A3.2: The distribution area of ​​the direct wave arrival time can be formed according to the time-distance relationship of the direct wave; 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 path with the maximum cumulative reward value is obtained; starting from the gather with the smallest shot offset, the direct wave tracking detection is performed in the constrained area.

[0098] A3.3: The distribution area of ​​the first-order and second-order water-related multiple waves can be formed according to the time-distance relationship of the first-order and second-order water-related multiple waves; 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 path with the maximum cumulative reward value is obtained; starting from the gather with the smallest shot offset, the direct wave tracking detection is carried out in the constrained area.

[0099] The following are several specific embodiments of the present invention.

[0100] Example 1

[0101] In a specific embodiment 1 of the present invention, a method for identifying water body wave phenomena based on the theoretical simulation of wave propagation rays in seawater is provided, which mainly involves the identification of direct waves, first-order and second-order water body related multiple waves in OBN observation data. The shot gather data used in embodiment 1 is as follows: Figure 2 As shown in the figure, it is the forward modeling data generated based on the horizontal seabed model. The model size is 10000m*600m, 121 shots are fired, and there are 91 simulated OBN nodes. The shot spacing is 50m, the receiver spacing is 100m, the maximum shot offset is 9000m, and the water depth is 250m.

[0102] Before execution, it is necessary to provide the OBN gather data file directory, output file directory, basic gather information, and prior parameters designed by the processor. The water body wave phenomenon identification method based on the simulation of wave propagation ray theory in seawater includes:

[0103] The first step is 1: input of OBN common detection point gathers and acquisition of trace header information.

[0104] Step 1 specifically includes the following steps:

[0105] A1.1: Read the OBN common detection point gather to obtain the header information of the P component data, including line number, shot coordinates, detection point 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] In step 2, ray tracing is performed using the trial shooting method. Taking the calculation of first-order water-related multiple waves as an example, step 2 specifically includes the following steps:

[0108] A2.1: How to describe 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 interior of the model is regarded as a uniform continuous medium. 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 is obtained by using cubic spline function interpolation.

[0110] The model interpolation interface in Example 1 is as follows Figure 2 as shown in .

[0111] According to the local slope of the interface, the slope of the upward outgoing ray can be derived. According to the intersection point of the interface and the slope of the upward ray, the coordinates of the intersection with the free surface and the slope of the ray when it enters the free surface can be calculated. Similarly, the slope of the downward ray when reflected from the free surface and the coordinates of the intersection when it reaches the seabed interface can be calculated.

[0112] A2.3: Adjust the initial emission direction of the ray at a certain emission angle interval, and calculate the initial emission angles of the two left and right rays closest to the target receiving point coordinates.

[0113] A2.4: Use the dichotomy method to further reduce the angle range of the target ray so that the coordinates of the intersection point of the theoretical outgoing ray reaching the interface are close to the coordinates 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 velocity background, calculate the travel time of the wave to the receiving point based on the length of the ray path.

[0115] The first-order multiple wave ray path simulated by the test shooting method in Example 1 is as follows Figure 2 shown.

[0116] A3.1: Calculate the characteristic attributes of the common detection point gathers, including energy attributes, envelope attributes, and waveform similarity attributes.

[0117] For any real signal f(t) there is a corresponding analytic signal The envelope property can be obtained through Hilbert, and the envelope of the data reflects the macroscopic changes 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] Assume a time window with N channels centered at the analysis point, and use the similarity coefficient r to measure the overall similarity of multiple channels in the time window;

[0122]

[0123] Where N represents the number of channels contained in the local window, x n represents the coordinate of the nth track, q represents the local scanning inclination of the event axis 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 waveform; Figure 3 It is a schematic diagram of calculating the similarity coefficient attributes of common detection point gathers in a specific embodiment of the present invention.

[0124] Construct a Markov decision state transition decision process and find a decision sequence that can produce the maximum cumulative reward value.

[0125]

[0126] where γ represents the discount factor that decays over time, p is the state transfer matrix, a is the current action, S represents the attribute used to describe the information of the selected water 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 in the next time step.

[0127] A3.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 track set with the smallest shot distance, the attribute value of each time point is calculated in the constraint area, and it is randomly moved up or down to another time point according to the action probability. If the expected cumulative reward value is larger, the time step is selected as the next state, and this process is repeated to perform direct wave tracking detection; the travel time difference constraint between adjacent tracks is introduced to ensure the overall smoothness of the travel time.

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

[0129] A3.3: According to the best decision principle, the paths of the first-order and second-order water-related multiple waves are tracked within the theoretical range of the first-order and second-order water-related multiple wave travel times. Starting from the trace set with the smallest shot offset, direct wave tracking detection is performed within the constraint area; the travel time difference constraint of adjacent traces is introduced between two adjacent traces to ensure the overall smoothness of the travel time.

[0130] Figure 5 Shown is a schematic diagram of the identification and travel time tracking of first-order water-related multiple waves in the simulated shot set of Example 1.

[0131] The model data used in Example 2 is as follows: Figure 6 As shown in the figure, it is the forward modeling data generated based on the model of the undulating seabed. The model size is 10000m*801m, 121 shots are fired, and there are 91 simulated OBN nodes. The shot spacing is 50m, the detection point spacing is 100m, the maximum shot offset is 9000m, the time sampling rate is 0.4ms, and the average water depth is 300m. The travel time of the first-order multiple waves in the water wave field is tracked using the above wave phenomenon identification process.

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

[0133] Figure 7 It is a schematic diagram of calculating similarity coefficient attributes according to the shot gather of the undulating seabed water body model in this embodiment 2;

[0134] Figure 8 Schematic diagram of the first-order multiple wave travel time picking interval obtained according to the water body model constraint in this embodiment 2;

[0135] Fig. 9 Schematic diagram of the identification and travel time tracking of first-order water-related multiple waves in the simulated shot gather in Example 2.

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

[0137] Fig.10 This is a schematic diagram of simulating first-order water-related multiple waves by the test-firing method based on the observation system in Example 3;

[0138] Fig.11 Schematic diagram of the first-order multiple wave travel time picking interval obtained according to the water body model constraint in specific embodiment 3;

[0139] Fig.12 Schematic diagram of the actual OBN data direct wave and first-order water-related multiple wave travel time tracking in Example 3.

[0140] Finally, it should be noted that the above is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art can still modify the technical solutions recorded in the aforementioned embodiments or replace some of the technical features therein with equivalents. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

[0141] Except for the technical features described in the specification, all other technical features are known technologies to those skilled in the art.

Claims

1. A method for identifying water body wave phenomena based on the simulation of wave propagation ray theory in seawater, characterized in that: The method for identifying water body wave phenomena based on theoretical simulation of wave propagation rays in seawater includes: Step 1, obtain the depth values ​​of all OBN nodes on a single line; Step 2, simulating the wave field propagating in water by the test shooting method; Step 3: Track and identify the wave field based on Markov optimal decision: Step 4: Output the water body wave field travel time to obtain the result.

2. The method for identifying water body wave phenomena based on theoretical simulation of wave propagation rays in seawater according to claim 1 is characterized in that: In step 1, the OBN common detection point gather is read to obtain the header information of its P component data and the depth values ​​of all OBN nodes on the single line.

3. The method for identifying water body wave phenomena based on theoretical simulation of wave propagation rays in seawater according to claim 2 is characterized in that: Step 1 includes: A1.1: Read the OBN shot line and obtain the trace header information, including the line number, shot coordinates, detection point coordinates, elevation, offset, sampling rate, and sampling time; A1.2: Record the depth values ​​of all OBN nodes on a single line.

4. The method for identifying water body wave phenomena based on theoretical simulation of wave propagation rays in seawater according to claim 1 is characterized in that: 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 value; based on the initial seabed model, the theoretical first-order and second-order water-related multiple wave arrival times of each shot and each OBN node are calculated by the trial firing method.

5. The method for identifying water body wave phenomena based on theoretical simulation of wave propagation rays in seawater according to claim 4 is characterized in that: Step 2 includes: A2.1: According to the OBN node depth value, the interface of the water model is expressed in the form of an analytical function; A2.2: Given the starting point of the ray and an assumed direction of ray emission, the ray will be reflected according to the geometric properties after reaching the interface. In this way, a ray is traced, and its intersection point when reflected on the interface and its final landing point on the interface are calculated; A2.3: Adjust the initial emission direction of the ray at a certain emission angle interval, and find the initial emission angles of the two rays closest to the target receiving point; A2.4: Use the dichotomy method to further reduce the angle range of the target ray so that the coordinates of the intersection point of the theoretical outgoing ray reaching the interface are close to the coordinates of the target receiving point, and record the intersection point when the ray passes through the interface; A2.5: Assuming a constant velocity background, calculate the travel time of the wave to the receiving point based on the length of the ray path; A2.6: The theoretical time-distance relationship of the first-order and second-order water wave fields on the common detection point gathers is obtained by data rearrangement.

6. The method for identifying water body wave phenomena based on theoretical simulation of wave propagation rays in seawater according to claim 5, characterized in that: 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 obtained by linear interpolation, polynomial interpolation, and cubic spline function interpolation; assuming that there are n+1 OBN nodes, divided into n intervals, and a cubic spline function is used to describe a point x∈[x k ,x k+1 Depth of: 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 They are the constant term, linear term coefficient, quadratic term coefficient, and cubic term coefficient of the interpolation function respectively; 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, the endpoints f i (x i ),f i '(x i ),f i ”(x i ) is continuous; therefore, all internal endpoints need to satisfy the cubic equations of the two left and right segments; the equations for the first and second order derivatives of n-1 internal points are continuous: f i '(x i ) = f i ' +1 (x i ),f i ”(x i ) = f i ” +1 (x i ); Assume natural boundary conditions, that is, f0”(x0)=f n ”(x n )=0, then a total of 4n equations are formed, and the linear equation is obtained by deduction: where h i =x i+1 -x i ,m i =f i ”(x i ), h i Indicates the step length, m i is the second-order derivative at the node; construct a system of linear equations: Use Jacobi iteration method to solve, and then substitute m back into the coefficients of the cubic spline function; the formula is as follows: For three-dimensional geological bodies, the interface is described by triangular patches or B-spline surfaces.

7. The method for identifying water body wave phenomena based on theoretical simulation of wave propagation rays in seawater according to claim 5, characterized in that: In step 2.2, given a ray starting point position and the initial ray emission direction, calculate the down-going wave process and the up-going wave process of the ray propagation path; given the down-going ray slope k1 and intercept b1, calculate the intersection coordinates of the down-going ray and the seabed interface of the water model by the simultaneous formula: Where y is the depth of the interpolation point, a i ,b i ,c i ,d i They are the constant term, linear term coefficient, quadratic term coefficient, and cubic term coefficient of the interpolation function respectively; Given the coordinates of the incident point of the ray and the coordinates of the interface intersection point, according to the slope k1 of the incident ray, the direction vector of the incident ray is It is also necessary to find the direction vector of the outgoing ray after reflection According to the law of reflection, the incident angle is equal to the outgoing angle; the direction vector of the incident ray and the direction vector of the reflected ray are symmetric with respect to the normal vector; find the direction vector of the outgoing ray: in is the local slope of the interface at the intersection, the coordinates of the ray incident point are (x1, y1), the coordinates of the interface intersection are (x2, y2), and the direction of the incident ray is The direction of the outgoing ray is According to the local slope of the interface, the slope of the upward outgoing ray can be derived. According to the intersection point of the interface and the slope of the upward ray, the coordinates of the intersection with the free surface and the slope of the ray when it enters the free surface can be calculated. Similarly, the slope of the downward ray when reflected from the free surface and the coordinates of the intersection when it reaches the seabed interface can be calculated.

8. The method for identifying water body wave phenomena based on theoretical simulation of wave propagation rays in seawater according to claim 5, characterized in that: In step 2, the processing of step 2 is implemented for each data of each common shot gather, and the theoretical time-distance relationship of the first-order water-related multiple waves of each common shot gather can be plotted; the simulation of the second-order water-related multiple waves is similar to the first-order, and one upgoing wave and one downgoing wave are added to the first-order simulation; and the theoretical time-distance relationship of the water wave field on the common detection point gather can be obtained by data rearrangement.

9. The method for identifying water body wave phenomena based on theoretical simulation of wave propagation rays in seawater according to claim 1, characterized in that: In step 3, when the direct wave, first-order and second-order water-related multiple waves of the calculated simulated OBN data arrive, the direct wave and first-order and second-order water-related multiple waves in the common detection point gathers of the actual data are identified based on the Markov optimal decision under the model constraint in each common detection point gather.

10. The method for identifying water body wave phenomena based on theoretical simulation of wave propagation rays in seawater according to claim 9, characterized in that: Step 3 includes: A3.1: Generate characteristic attributes of each common detection point gather, including energy attributes, envelope attributes, and waveform similarity attributes; A3.2: The distribution area of ​​the direct wave arrival time can be formed according to the time-distance relationship of the direct wave; 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 path with the maximum cumulative reward value is obtained; starting from the gather with the smallest shot offset, the direct wave tracking detection is carried out in the constrained area; A3.3: The distribution area of ​​the first-order and second-order water-related multiple waves can be formed according to the time-distance relationship of the first-order and second-order water-related multiple waves; 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 path with the maximum cumulative reward value is obtained; starting from the data set with the smallest shot offset, the direct wave tracking detection is carried out in the constrained area.

11. The method for identifying water body wave phenomena based on theoretical simulation of wave propagation rays in seawater according to claim 10, characterized in that: In step 3.1, since the energy of the direct wave decays slowly in the water body, it is the first strong energy wave to arrive at the OBN node. The energy characteristics of the seismic traces before and after the direct wave should change greatly. Set a time window, divide the time window into two parts, and determine the arrival time of the strong energy wave field by the energy ratio between the front and rear parts of the time window. For any real signal f(t) there is a corresponding analytic signal The envelope property can be obtained through Hilbert, and the envelope of the data reflects the macroscopic changes of the waveform in the time domain; 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); Assume a time window with N channels centered at the analysis point, and use the similarity coefficient r to measure the overall similarity of multiple channels in the time window; Where N represents the number of channels contained in the local window, x n represents the coordinate of the nth track, q represents the local scanning inclination of the event axis 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 waveform; The Markov decision process is a learning process in which an agent takes actions based on the environment to change its state and obtain rewards; It is generally composed of five elements<S,A,T,P,r> For the tracking of direct waves, the state set S represents the current time part t of the execution decision, which is used to describe the information of the selected direct wave travel time position; action A represents the estimation of the next direct wave position according to the action probability π(a|s) under the given model and time-distance relationship constraints. There are generally only two actions, one is to move a sampling point in the direction of decreasing travel time, and the other is to move a sampling point in the direction of increasing travel time. The action probability π can be regarded as random; the state transfer matrix P represents the possibility of transferring to state s' according to the attribute information of the current state s when implementing action A. According to the theoretical travel time difference control of two adjacent paths, the overall distribution shows that the greater the difference from the theoretical travel time difference, the lower the transfer probability; the reward function R represents the instantaneous reward value r obtained after the A transfer is realized, which is constructed by the characteristic attributes: Where M represents the total number of attributes used, w represents the weight of the attribute, f represents the attribute value of the current position, and r is the instantaneous reward value of the current point; The time series 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 the discount factor that decays over time, p is the state transfer matrix, a is the current action, S represents the attribute used to describe the information of the selected water 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 in the next time step.

12. The method for identifying water body wave phenomena based on theoretical simulation of wave propagation rays in seawater according to claim 10, characterized in that: 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 track set with the smallest shot offset, the attribute value of each time point is calculated in the constraint area, and it is randomly moved up or down to another time point according to the action probability. If the expected cumulative reward value is larger, the time step is selected as the next state, and this process is repeatedly executed to perform direct wave tracking detection; the introduction of the adjacent track travel time difference constraint between two adjacent tracks ensures the overall smoothness of the travel time.

13. The method for identifying water body wave phenomena based on theoretical simulation of wave propagation rays in seawater according to claim 10, characterized in that: In step 3.3, according to the best decision principle, the paths of the first-order and second-order water-related multiple waves are tracked within the theoretical range of the first-order and second-order water-related multiple wave travel times; starting from the track set with the smallest shot offset, direct wave tracking detection is performed in the constrained area; and the adjacent track travel time difference constraint is introduced between two adjacent tracks to ensure the overall smoothness of the acquired travel time.

Images