Residual multiple construction and attenuation method based on Rayleigh integral forward modeling

Abstract

The invention discloses a remnant multiple construction and attenuation method based on Rayleigh integral forward modeling, and the method comprises the following steps: S1, building a high-resolution pre-stack time migration velocity spectrum, picking up multiple velocity, and building a multiple migration velocity field; s2, calculating three-dimensional pre-stack time migration based on a multiple migration velocity field; s3, constructing a residual multiple record based on time domain Rayleigh integral forward modeling; and S4, residual multiple suppression based on curvelet domain least square filtering. Through multiple record construction and adaptive subtraction processing of the method, residual multiples in the data can be further eliminated, and the signal-to-noise ratio of seismic data is remarkably improved. The method is beneficial supplement to a conventional multiple attenuation method, and a new way is provided for multiple suppression.

Description

Technical Field

[0001] The present invention belongs to the technical field of seismic data processing and analysis, and particularly relates to a method for constructing and attenuating residual multiple waves based on Rayleigh integral forward modeling. Background Art

[0002] At present, the seismic exploration method applied to oil and gas exploration is mainly the reflection wave method, which generally considers the primary reflection to be effective, while the direct wave, shallow refraction wave, multiple wave, etc. are all interference waves. For the interference of direct waves and shallow refraction waves, etc., they can generally be directly removed on the record, while the identification and suppression of multiple waves have been one of the hot issues that people have been focusing on for a long time. For exploration areas with multiple strong wave impedance interfaces, since both the sea surface and the sea floor are strong wave impedance interfaces, and there will be strong shielding layers of basement and limestone in the middle and deep parts, the corresponding seismic records usually contain a large number of multiple wave interferences of various types and strong amplitudes, completely covering the weak amplitude primary reflection signals in the deep part.

[0003] The current multiple wave elimination methods mainly include predictive deconvolution technology, multiple wave suppression methods based on apparent velocity differences, and multiple wave suppression methods based on wave theory, etc. Predictive deconvolution has the advantage of high efficiency, but the period of multiple waves shows time-varying and space-varying characteristics. Therefore, an accurate prediction step length cannot be given, resulting in residual multiple waves or damage to effective signals, and it is difficult to effectively attenuate long-period multiple waves.

[0004] The multiple wave suppression method based on apparent velocity differences has been widely used in actual seismic data processing. However, due to the small time difference in the near-offset traces, it often leads to serious residuals of the in-phase axis of strong amplitude multiple waves, and the additional processing steps such as NMO and de-NMO introduced often cause waveform distortion in the processing results.

[0005] In the multiple wave suppression method based on wave theory, the free surface multiple wave attenuation method (SRME) has become one of the preferred methods for seismic data processing, but it is difficult to effectively suppress multiple waves in the far-offset traces.

[0006] Therefore, for seismic data with multiple strong wave impedance interfaces, even after the combined suppression of the above-mentioned multiple wave attenuation methods, there will still be obvious residual multiple waves in the data. If not eliminated, it will seriously affect the authenticity and reliability of seismic imaging, and further mislead the subsequent seismic geological interpretation. Summary of the Invention

[0007] The problem to be solved by the present invention is to provide a method for constructing and attenuating residual multiple waves based on Rayleigh integral forward modeling. The construction of multiple wave records and the adaptive subtraction processing of this method can further eliminate the residual multiple waves in the data and significantly improve the signal-to-noise ratio of seismic data.

[0008] To solve the above technical problems, the technical solution adopted by the present invention is: a method for constructing and attenuating residual multiple waves based on Rayleigh integral forward modeling, comprising the following steps:

[0009] S1: Establish a high-resolution prestack time migration velocity spectrum, pick up the multiple wave velocity, and establish a multiple wave migration velocity field;

[0010] S2: Calculate the 3D prestack time migration based on the multiple wave migration velocity field;

[0011] S3: Construct a residual multiple wave record based on time-domain Rayleigh integral forward modeling simulation;

[0012] S4: Suppress the residual multiple waves based on curvelet domain least squares filtering.

[0013] Further, in the S1, the establishment of the high-resolution prestack time migration velocity spectrum includes the following steps:

[0014] S11: Determine the starting migration velocity v0, velocity interval Δv, and maximum velocity v according to the velocity range of the multiple waves in the seismic data. max , then the series of velocities v k The formula is as follows:

[0015] v k = v0 + kΔv (Formula 1)

[0016] In the formula, k is the velocity serial number, 0 ≤ k ≤ K, K is the serial number corresponding to the maximum velocity, and satisfies the following relational expression:

[0017]

[0018] S12: The prestack time migration velocity spectrum is the migration result of the series of velocities. Based on the series of velocities defined in Formula (1), input the single-shot record d i (j,t) for Kirchhoff prestack time migration in the shot domain. The formula is as follows:

[0019]

[0020] In the formula, IG k is the imaging gather migrated using the velocity v k ; P is the coordinate of the imaging point corresponding to the spatial position of the velocity spectrum and the travel time within the spectrum; W M is the weight factor of Kirchhoff migration, τ represents the sum of the travel times of the ray from the shot point to the imaging point and from the imaging point to the geophone, (S i , G j ) is the coordinate of a shot-receiver pair; d iIt is a single-shot record, where i and j are the shot number and trace number respectively, t represents the travel time, and I and J are the number of shots and traces applied during migration;

[0021] S13: The weight factor W M (v k , S i , G j , P) and the travel time τ(v k , S i , G j , P) both depend on the velocity v k , and they are respectively expressed as follows,

[0022]

[0023] In the formula, r0 represents the spatial distance from the shot point S i to the imaging point P, and r is the spatial distance from the imaging point P to the geophone G j .

[0024] S14: Based on formula (3), the common - reflection - surface stacking weighting process is introduced into the velocity spectrum calculation. The formula for calculating the common - reflection - surface stacking weighting factor for the imaging gather is as follows,

[0025]

[0026] In the formula, F(v k , P) represents the common - reflection - surface stacking weighting factor with respect to the velocity v k ; i is the trace number in the imaging gather, I represents the total number of traces; l is the time - window length; ζ is a constant to ensure that the denominator is not zero;

[0027] S15: For each imaging gather corresponding to the scanning velocity value v k , the data of each trace are horizontally stacked and the absolute value is taken to obtain the migration velocity spectrum. Using formula (5), the weighting factor is used to perform common - reflection - surface stacking weighting on the above result, and this process is expressed as,

[0028]

[0029] In the formula, V(v k , P) represents the high - resolution migration velocity spectrum, P is the coordinate of the imaging point corresponding to the spatial position of the velocity spectrum and the travel time within the spectrum; i is the trace number in the imaging gather, and I represents the total number of traces.

[0030] Furthermore, in the above - mentioned S1, the establishment of the multiple - wave migration velocity field includes the following steps,

[0031] S16: For the offset velocity spectrum at a certain spatial position, plot a two-dimensional image with the velocity value v on the abscissa and the travel time t on the ordinate. Identify the energy of the residual multiple waves in the velocity spectrum and pick out the corresponding positions to obtain the velocity curve of the residual multiple waves.

[0032] S17: Interpolate in time and space according to the velocity curves at different positions to create a multiple-wave migration velocity field.

[0033] Further, the S2 includes the following steps.

[0034] S21: According to the multiple-wave migration velocity field, use the prestack seismic data to generate a multiple-wave imaging data volume. The Kirchhoff prestack migration in the shot gather domain is expressed as an integral summation process in the following form.

[0035]

[0036] In the formula, M represents the multiple-wave prestack time migration data volume, P is the data volume sample point coordinate; W M is the weight factor of the Kirchhoff integral migration, τ represents the sum of the travel times of the ray from the shot point to the imaging point and from the imaging point to the geophone, (S i , G j ) is the coordinate of a shot-geophone pair; d i is the shot gather record containing the residual multiple waves, i and j are the shot number and trace number respectively, t represents the travel time, and I and J are the number of shots and traces applied during migration.

[0037] S22: Simplify the calculation processes of the weight factor W M (S i , G j , P) and the travel time τ(S i , G j , P). The two are expressed as follows.

[0038]

[0039] In the formula: v(P) is the migration velocity of the multiple waves, cosθ represents the tilt factor, r0 is the distance from the shot point to the imaging point, and r is the distance from the imaging point to the shot point.

[0040] Further, the S3 includes the following steps.

[0041] S31: Input the multiple-wave data volume M(P) generated according to formula (7), introduce the direct path ray tracing process based on the migration velocity model, and the forward modeling process of the time-domain Rayleigh integral is expressed as follows.

[0042]

[0043] where m is the multiple wave record constructed by forward modeling using the time-domain Rayleigh integral, (S i , G j ) is the coordinate of a shot-receiver pair, i and j represent the shot number and trace number respectively; M(P) represents the multiple wave imaging data volume, P is the profile sample point coordinate; ds represents the integration surface element; R is the weight factor of the time-domain Rayleigh integral, τ represents the sum of the travel times of the ray from the shot point S i to point P and from point P to the receiver point G j ;

[0044] S32: In the forward modeling process of the time-domain Rayleigh integral shown in formula (9), the integration weight factor R(P, S i , G j ) adopts the exact diffraction superposition calculation, and the formula is as follows,

[0045]

[0046] where r0 and r are the propagation distances of the incident ray and the diffracted ray respectively; θ0 and θ are the angles between r0, r and the normal vector n of the surface element; v(P) represents the migration velocity of the multiple waves, and P is the data volume sample point coordinate.

[0047] Furthermore, the said S4 includes the following steps,

[0048] S41: Perform adaptive subtraction of multiple waves in the curvelet domain with sparse characteristics. Transform the original seismic record d(x, t) and the predicted multiple wave record m(x, t) into the curvelet domain to obtain their curvelet coefficients c d (j, l, k) and c m (j, l, k), and the formula is as follows,

[0049]

[0050] where represents the curvelet basis function; c d (j, l, k) and c m (j, l, k) represent the curvelet coefficients of the original seismic record and the multiple wave record respectively, and j, l, k represent the scale, direction and position of the curvelet respectively; x and t represent the position of the seismic trace and the travel time of the reflected signal respectively;

[0051] S42: For the sample point at position k in the two-dimensional curvelet coefficient matrix of any scale j and direction l, set a rectangular window centered on it, intercept the data block therein for least squares filtering, and the result of subtracting the multiple waves is expressed as follows,

[0052] c p (n) = c d (n) - f * c m(n) (Formula 12)

[0053] Wherein, c p (n) is the recorded sample point after attenuating the multiple waves; f is the filtering factor to be obtained; c d (n), c m (n) are respectively the sample points in the intercepted record block, and n is the sample point serial number in the intercepted record block;

[0054] S43: For the ghost wave adaptive subtraction processing based on the L2 norm, the filtering factor is determined by minimizing the sum of the squared error e in the following formula. The formula is as follows,

[0055] e = ||c d - c m f||2 (Formula 13)

[0056] Wherein, the vector c d and c m respectively represent the vectors composed of the original record and the multiple wave record data block. The filtering factor f is obtained by solving Formula (13). When the length of the filtering factor f is 1, the solution of the linear equation system can be avoided, and the filtering factor can be directly obtained. The formula is as follows,

[0057]

[0058] Wherein, T represents the transpose of the matrix (or vector);

[0059] S44: Substitute the filtering factor f obtained based on Formula (14) into Formula (12) to obtain the curvelet coefficient record block c p (n) after suppressing the residual multiple waves, and assign its central sample point to the central position of the original rectangular window in the new curvelet coefficient c p (j, l, k). Perform the suppression processing of S42 - S43 on each sample point in c p (j, l, k) to obtain the seismic record curvelet coefficient c p (j, l, k) after attenuating the residual multiple waves. Then, obtain the residual multiple wave suppression result p(x, t) in the space - time domain through the inverse curvelet transform. The formula is as follows,

[0060]

[0061] Wherein, p(x, t) is the seismic record after attenuating the residual multiple waves; represents the curvelet basis function; c p (j, l, k) is the curvelet coefficient after attenuating the residual multiple waves. j, l, and k respectively represent the scale, direction, and position of the curvelet; x and t respectively represent the position of the seismic trace and the travel time of the reflection signal.

[0062] Further, the method further includes the following steps:

[0063] Input the single-cable data in the original shot gather record and the single-cable data in the residual multiple shot gather record into Equation (11) respectively, and use Equation (11) to transform them into the curvelet domain.

[0064] For the samples at each scale, direction, and position in the curvelet domain, perform multiple wave matching attenuation based on the least squares filtering process described in Equations (12)-(14).

[0065] Use Equation (15) to transform the curvelet domain record with residual multiples eliminated back to the spatio-temporal domain.

[0066] Further, the present invention provides a device for running the above data processing method.

[0067] Further, the present invention provides a device including a memory, a processor, and an algorithm stored in the memory and executable on the processor. When the processor executes the computer program, the above data processing method is implemented.

[0068] Further, the present invention provides a computer-readable storage medium storing a computer algorithm, and when the computer algorithm is executed by a processor, the above data processing is implemented.

[0069] The present invention has the following advantages and positive effects:

[0070] For exploration areas with multiple strong wave impedance interfaces, seismic records usually contain a large number of multiple wave interferences of various types and strong amplitudes, completely masking the weak amplitude primary reflection signals in the deep part. After the combined attenuation of multiple methods such as predictive deconvolution technology, multiple wave suppression methods based on apparent velocity differences, and multiple wave suppression methods based on wave theory, obvious multiple wave residuals still exist in the data. Through the construction and adaptive subtraction processing of multiple wave records by the method of the present invention, the residual multiple waves in the data can be further eliminated, significantly improving the signal-to-noise ratio of seismic data. The present invention is a beneficial supplement to conventional multiple wave attenuation methods and provides a new way for multiple wave suppression. BRIEF DESCRIPTION OF THE DRAWINGS

[0071] Figure 1 is a schematic diagram of the overall process flow of an embodiment of the present invention.

[0072] Figure 2 is an example diagram of a shot gather record containing residual multiple waves in an embodiment of the present invention.

[0073] Figure 3 is a schematic diagram of the ray tracing process of the "direct ray" model in prestack time migration in an embodiment of the present invention.

[0074] Figure 4 It is the pre-stack time migration velocity spectrum containing residual multiple wave energy and the picked multiple wave velocity curve in the embodiment of the present invention.

[0075] Figure 5 It is an example diagram of a profile along the main survey line direction in the residual multiple wave migration data volume in the embodiment of the present invention.

[0076] Figure 6 It is a schematic diagram of the ray tracing process of the forward simulation of the Rayleigh integral in the time domain in the embodiment of the present invention.

[0077] Figure 7 It is the residual multiple wave record constructed by using the forward simulation of the Rayleigh integral in the time domain in the embodiment of the present invention.

[0078] Figure 8 It is a comparison diagram of single-cable data before and after attenuating the residual multiple waves in the embodiment of the present invention. Detailed implementation manners

[0079] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are some, but not all, of the embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0080] The embodiments of the present invention will be further described below with reference to the accompanying drawings:

[0081] As Figure 1 shown, a method for constructing and attenuating residual multiple waves based on forward simulation of the Rayleigh integral includes the following steps.

[0082] S1: Establish a high-resolution pre-stack time migration velocity spectrum, pick the multiple wave velocity, and establish a multiple wave migration velocity field. Specifically, S1 includes the following steps.

[0083] S11: Determine the starting migration velocity v0, velocity interval Δv, and maximum velocity v according to the velocity range of multiple waves in the seismic data. max , then the series of velocities v k The formula is as follows.

[0084] v k = v0 + kΔv (Formula 1)

[0085] In the formula, k is the velocity serial number, 0 ≤ k ≤ K, K is the serial number corresponding to the maximum velocity, and the following relational expression is satisfied.

[0086]

[0087] S12: Essentially, the prestack time migration velocity spectrum is the migration result of a series of velocities. Based on the series of velocities defined in formula (1), the single-shot record d i (j,t) is subjected to Kirchhoff prestack time migration in the shot domain, and the formula is as follows,

[0088]

[0089] where IG k is the image gather migrated using velocity v k ; P is the coordinate of the imaging point corresponding to the spatial position of the velocity spectrum and the travel time within the spectrum; W M is the weight factor for Kirchhoff migration, τ represents the sum of the travel times of the ray from the shot point to the imaging point and from the imaging point to the geophone, (S i ,G j ) is the coordinate of a shot-geophone pair; d i is the single-shot record, i and j are the shot number and trace number respectively, t represents the travel time, and I and J are the number of shots and traces applied during migration.

[0090] S13: Both the weight factor W M (v k ,S i ,G j ,P) and the travel time τ(v k ,S i ,G j ,P) depend on velocity v k . Since the ray tracing process of Kirchhoff prestack time migration adopts the "direct ray" model, as Figure 3 shown, they are respectively expressed as follows,

[0091]

[0092] where r0 represents the spatial distance of the ray from the shot point S i to the imaging point P, and r is the spatial distance from the imaging point P to the geophone G j .

[0093] S14: Formula (3) gives the migration result of the single-shot record based on velocity v k , which has low velocity resolution. It is necessary to introduce the in-phase weighting process into the velocity spectrum calculation. The formula for the in-phase weighting factor of the image gather is as follows,

[0094]

[0095] where F(v k ,P) represents with respect to velocity v kThe in-phase weighting factor; i is the trace number in the imaging gather, and I represents the total number of traces; l is the window length; ζ is a constant to ensure that the denominator is not zero, generally taking 0.01 - 0.001 of the average amplitude.

[0096] S15: For each scanning velocity value v k The corresponding imaging gather, horizontally stack the data of each trace and take the absolute value to obtain the migration velocity spectrum, and use formula (5) to calculate the weighting factor to perform in-phase weighting on the above results. This process is expressed as

[0097]

[0098] In the formula, V(v k ,P) represents the high-resolution migration velocity spectrum, P is the coordinate of the imaging point corresponding to the spatial position of the velocity spectrum and the travel time within the spectrum; i is the trace number in the imaging gather, and I represents the total number of traces.

[0099] S15: For the migration velocity spectrum at a certain spatial position, it can be plotted as a Figure 4 two-dimensional image with the velocity value v on the abscissa and the travel time t on the ordinate as shown. Identify the energy of the residual multiple waves in the velocity spectrum, pick up the corresponding positions, and obtain the curve of the velocity (v) of the residual multiple waves changing with the travel time (t), that is, the velocity curve.

[0100] S16: Interpolate in time and space according to the velocity curves at different positions to create a multiple-wave migration velocity field.

[0101] S2: Calculate the three-dimensional prestack time migration based on the multiple-wave migration velocity field. Specifically, S2 includes the following steps

[0102] S21: Kirchhoff prestack time migration is a mature migration imaging method in current production. In the case of creating a multiple-wave migration velocity field, use the prestack seismic data to generate a multiple-wave imaging data volume. The Kirchhoff prestack migration in the shot gather domain is expressed as an integral summation process in the following form

[0103]

[0104] In the formula, M represents the multiple-wave prestack time migration data volume, P is the data volume sample point coordinate; W M is the weighting factor of the Kirchhoff integral migration, τ represents the sum of the travel times of the ray from the shot point to the imaging point and from the imaging point to the geophone, (S i ,G j ) is the coordinate of a shot - geophone pair; d i is the shot gather record containing residual multiple waves, i and j are the shot number and trace number respectively, t represents the travel time, and I and J are the number of shots and traces applied during migration respectively.

[0105] S22: The ray tracing process of Kirchhoff prestack time migration adopts the "direct ray" model. As shown in Figure 3 , therefore, the weight factor W in formula (7) M (S i , G j , P) and the calculation process of the travel time τ(S i , G j , P) are significantly simplified. The two are expressed as

[0106]

[0107] where: v(P) is the migration velocity of the multiple wave, cosθ represents the tilt factor, r0 is the distance from the shot point to the imaging point, and r is the distance from the imaging point to the shot point.

[0108] S3: Construct the residual multiple wave record based on the forward simulation of the Rayleigh integral in the time domain. Specifically, S3 includes the following steps

[0109] S31: Compared with the conventional de-migration technique, the forward simulation of the Rayleigh integral in the time domain, which is a quasi-de-migration method, has higher computational efficiency and is more suitable for constructing the residual multiple wave shot gather records of 3D seismic exploration. Input the multiple wave data volume M(P) generated according to formula (7), introduce the direct path ray tracing process based on the migration velocity model, and the forward simulation process of the Rayleigh integral in the time domain is expressed as follows

[0110]

[0111] where m is the multiple wave record constructed by the forward simulation of the Rayleigh integral in the time domain, (S i , G j ) is the coordinate of a shot-receiver pair, i and j represent the shot number and the trace number respectively; M(P) represents the multiple wave imaging data volume, P is the profile sample point coordinate; ds represents the integral surface element; R is the weight factor of the Rayleigh integral in the time domain, and τ represents the sum of the travel times of the ray from the shot point S i to the point P and from the point P to the receiver point G j .

[0112] S32: In the forward simulation process of the Rayleigh integral in the time domain shown in formula (9), the integral weight factor R(P, S i , G j ) adopts the accurate diffraction stacking calculation, as shown in Figure 6 , that is

[0113]

[0114] In the formula, r0 and r are the propagation distances of the incident ray and the diffracted ray, respectively; θ0 and θ are the angles between r0, r and the normal vector n of the surface element; v(P) represents the migration velocity of the multiple wave, and P is the sample point coordinate of the data volume.

[0115] S4: Residual multiple wave suppression based on curvelet domain least squares filtering. Specifically, the S4 includes the following steps:

[0116] S41: To avoid damaging the effective signal, the multiple wave adaptive subtraction is selected in the curvelet domain with sparse characteristics. Therefore, the original seismic record d(x,t) and the predicted multiple wave record m(x,t) need to be transformed into the curvelet domain to obtain the curvelet coefficients c d (j, l, k) and c m (j, l, k), and the formula is as follows:

[0117]

[0118] In the formula: represents the curvelet basis function; c d (j, l, k) and c m (j, l, k) represent the curvelet coefficients of the original seismic record and the multiple wave record, respectively. j, l, and k represent the scale, direction, and position of the curvelet, respectively; x and t represent the position of the seismic trace and the travel time of the reflection signal.

[0119] S42: To ensure the stability of the filtering process, for the sample point at position k in the two-dimensional curvelet coefficient matrix of any scale j and direction l, a rectangular window is set centered on it, and the data block therein is intercepted for least squares filtering. The result of subtracting the multiple wave is expressed as follows:

[0120] c p (n) = c d (n) - f * c m (n) (Formula 12)

[0121] In the formula, c p (n) is the recorded sample point after attenuating the multiple wave; f is the filtering factor to be obtained; c d (n), c m (n) are the sample points in the intercepted record block respectively, n is the sample point serial number in the intercepted record block. Generally, to ensure the stability of equation solving, it is usually required that n ≥ 4.

[0122] S43: For the ghost wave adaptive subtraction processing based on the L2 norm, the filtering factor is determined by minimizing the sum of squared error energy e in the following formula. The formula is as follows:

[0123] e = ||c d - c mf||2 (Formula 13)

[0124] Among them, the vector c d and c m respectively represent the vectors composed of the original record and the multiple wave record data block. The filtering factor f is obtained by solving Formula (13). When the length of the filtering factor f is 1, the solution of the linear equation system can be avoided, and the filtering factor can be directly obtained, that is

[0125]

[0126] In the formula, T represents the transpose of the matrix (or vector).

[0127] S44: Substitute the filtering factor f obtained according to Formula (14) into Formula (12) to obtain the curvelet coefficient record block c p (n) after suppressing the residual multiples, and assign its central sample point to the new curvelet coefficient c p (j, l, k) at the central position of the original rectangular window. Perform the suppression processing of S42 - S43 on each sample point in c p (j, l, k) to obtain the seismic record curvelet coefficient c p (j, l, k) after attenuating the residual multiples. Then, obtain the suppression result p(x, t) of the residual multiples in the space - time domain through the inverse curvelet transform, that is

[0128]

[0129] In the formula, p(x, t) is the seismic record after attenuating the residual multiples; represents the curvelet basis function; c p (j, l, k) are the curvelet coefficients after attenuating the residual multiples. j, l, and k respectively represent the scale, direction, and position of the curvelet; x and t respectively represent the position of the seismic trace and the travel time of the reflected signal.

[0130] The following combines specific embodiments to specifically elaborate on the present invention:

[0131] Sea area A is a shallow - water area with a relatively flat seabed, a water depth of about 50m, and there is a strong shielding layer between 500 - 1000ms underground. The seismic waves are continuously excited in a moving - ship mode by an air - gun source towed at the stern of the ship, and 6 streamers are used to receive seismic signals with a streamer spacing of 100m; the sampling interval and record length of the seismic data are 0.004 seconds and 8 seconds respectively. Due to the existence of multiple strong interfaces, there are a large number of free - surface multiples and inter - layer multiples with strong amplitudes in the seismic data, and diffracted multiples are formed at the positions where the underground strong interfaces fluctuate violently. These multiples seriously affect the imaging quality of the underground structure and will mislead the subsequent geological interpretation and analysis.

[0132] The specific steps are as follows

[0133] S1: Establishment of high-resolution prestack time migration velocity spectrum and picking of multiple velocities.

[0134] Specifically, determine the starting migration velocity v0, velocity interval Δv, and maximum velocity v according to the velocity range of multiples in the seismic data of the work area max , and let the corresponding values be 1000 m / s, 25 m / s, and 7000 m / s respectively, and calculate a series of velocities v using the following formula k ,

[0135] v k = v0 + kΔv (Formula 1)

[0136] In the formula, k is the velocity serial number, 0 ≤ k ≤ K, K is the serial number corresponding to the maximum velocity, and there is

[0137]

[0138] The prestack time migration velocity spectrum is the migration result of a series of velocities. Based on the series of velocities defined in formula (1), input the single-shot record d i (j,t), as Figure 2 shown, perform Kirchhoff prestack time migration in the shot domain.

[0139]

[0140] In the formula, IG k is the imaging gather migrated using velocity v k ; P is the coordinate of the imaging point corresponding to the spatial position of the velocity spectrum and the travel time within the spectrum; W M is the weight factor of Kirchhoff migration, τ represents the sum of the travel times of the ray from the shot point to the imaging point and from the imaging point to the geophone, (S i , G j ) is the coordinate of a shot-geophone pair; d i is the single-shot record, i and j are the shot number and trace number respectively, t represents the travel time, and I and J are the number of shots and traces applied during migration respectively.

[0141] The weight factor W M (v k , S i , G j , P) and the travel time τ(v k , S i , G j , P) both depend on velocity v k , because the ray tracing process of Kirchhoff prestack time migration adopts the "direct ray" model, as Figure 3 shown, the two can be respectively expressed as,

[0142]

[0143] In the formula, r0 represents the spatial distance from the shot point S of the ray i to the imaging point P, and r is the spatial distance from the imaging point P to the geophone G j .

[0144] Formula (3) gives the migration result of a single-shot record based on the velocity v k , which has a low velocity resolution. It is necessary to introduce the in-phase weighting process into the velocity spectrum calculation. The formula for the in-phase weighting factor for the imaging gather is as follows

[0145]

[0146] In the formula, F(v k , P) represents the in-phase weighting factor with respect to the velocity v k ; i is the trace number in the imaging gather, and I represents the total number of traces; l is the time window length; ζ is a constant to ensure that the denominator is not zero, generally taking 0.01 - 0.001 of the average amplitude.

[0147] For the imaging gather corresponding to each scanned velocity value v k , by horizontally stacking the data of each trace with equal τ and taking the absolute value, the migration velocity spectrum can be obtained, and the above result is weighted in phase using the weighting factor obtained by formula (5). This process can be expressed as

[0148]

[0149] In the formula, V(v k , P) represents the high-resolution migration velocity spectrum, P is the coordinate of the imaging point corresponding to the spatial position of the velocity spectrum and the travel time within the spectrum; i is the trace number in the imaging gather, and I represents the total number of traces.

[0150] For the migration velocity spectrum at a certain spatial position, it can be plotted as a two-dimensional image with the velocity value v on the abscissa and the travel time t on the ordinate as shown in Figure 4 . Identify the energy of the residual multiple waves in the velocity spectrum, and pick up the corresponding positions to obtain the curve of the velocity (v) of the residual multiple waves varying with the travel time (t) (i.e., the velocity curve, see the white curve in Figure 4 ), and then create the multiple-wave migration velocity field by interpolating in time and space according to the velocity curves at different positions.

[0151] S2: 3D prestack time migration based on the multiple-wave migration velocity field.

[0152] Specifically, Kirchhoff prestack time migration is a mature migration imaging method in current production. In the case of creating a multiple wave migration velocity field, it uses prestack seismic data to generate a multiple wave imaging data volume. The Kirchhoff prestack migration in the shot gather domain can be expressed as an integral summation process in the following form:

[0153]

[0154] In the formula, M represents the multiple wave prestack time migration data volume, P is the sample point coordinate of the data volume; W M is the weight factor of Kirchhoff integral migration, τ represents the sum of the travel times of the ray from the shot point to the imaging point and from the imaging point to the geophone, (S i , G j ) is the coordinate of a shot - geophone pair; d i is the shot gather record containing residual multiples, i and j are the shot number and trace number respectively, t represents the travel time, and I and J are the number of shots and traces applied during migration respectively.

[0155] The ray tracing process of Kirchhoff prestack time migration adopts the "direct ray" model. As Figure 3 shown, therefore, the calculation processes of the weight factor W M (S i , G j , P) and the travel time τ(S i , G j , P) are significantly simplified and can be expressed respectively as

[0156]

[0157] In the formula, v(P) is the migration velocity of multiples, cosθ represents the dip factor, r0 is the distance from the shot point to the imaging point, and r is the distance from the imaging point to the shot point.

[0158] According to Figure 4 the velocity curve picked up in, create a multiple wave migration velocity field, input the shot gather record shown in Figure 2 , and establish a residual multiple imaging data volume based on the Kirchhoff prestack time migration process described by formula (7) and formula (8). The profile along the main survey line is as Figure 4 shown, which contains residual multiple reflection signals below the strong reflection interface in the work area.

[0159] S3: Construct a residual multiple record based on forward modeling of the time - domain Rayleigh integral.

[0160] Specifically, compared with the conventional anti-offset technology, the forward modeling of the time-domain Rayleigh integral, a quasi-anti-offset method, has higher computational efficiency and is more suitable for constructing the residual multiple shot gather records of 3D seismic exploration. By inputting the multiple wave data volume M(P) generated according to formula (7) and introducing the direct-path ray tracing process based on the migration velocity model, the forward modeling process of the time-domain Rayleigh integral can be expressed as follows:

[0161]

[0162] where m is the multiple wave record constructed by the forward modeling of the time-domain Rayleigh integral, (S i , G j ) are the coordinates of a shot-receiver pair, i and j represent the shot number and trace number respectively; M(P) represents the multiple wave imaging data volume, P is the profile sample point coordinate; ds represents the integration surface element; R is the weight factor of the time-domain Rayleigh integral, and τ represents the sum of the travel times of the ray from the shot point S i to point P and from point P to the receiver point G j .

[0163] In the forward modeling process of the time-domain Rayleigh integral shown in formula (9), the integration weight factor R(P, S i , G j ) should be calculated using accurate diffraction summation, as shown in Figure 6 , that is

[0164]

[0165] where r0 and r are the propagation distances of the incident ray and the diffracted ray respectively; θ0 and θ are the angles between r0, r and the normal vector n of the surface element; v(P) represents the migration velocity of the multiple wave, and P is the data volume sample point coordinate.

[0166] Based on the migration velocity field of the residual multiple waves, by inputting the residual multiple wave data volume shown in Figure 5 , the shot gather record shown in Figure 7 is constructed using the forward modeling of the time-domain Rayleigh integral given by formula (9) and formula (10), which only contains the residual multiple waves below the strong reflection interface.

[0167] S4: Suppression of residual multiple waves based on curvelet domain least squares filtering.

[0168] To avoid damaging the effective signals, the multiple wave adaptive subtraction is selected in the curvelet domain with sparse characteristics. Therefore, the original seismic record d(x, t) and the predicted multiple wave record m(x, t) need to be transformed into the curvelet domain to obtain their curvelet coefficients c d (j, l, k) and c m (j, l, k), that is

[0169]

[0170] In the formula, represents the curvelet basis function; c d (j, l, k) and c m (j, l, k) represent the curvelet coefficients of the original seismic record and the multiple wave record respectively. j, l, and k represent the scale, direction, and position of the curvelet respectively; x and t represent the position of the seismic trace and the travel time of the reflection signal respectively.

[0171] To ensure the stability of the filtering process, for the sample point at position k in the two-dimensional curvelet coefficient matrix of any scale j and direction l, a rectangular window is set centered on it, and the data block therein is intercepted for least squares filtering. Then the result of subtracting the multiple wave can be expressed as

[0172] c p (n) = c d (n) - f * c m (n) (Formula 12)

[0173] In the formula, c p (n) is the recorded sample point after attenuating the multiple wave; f is the filtering factor to be obtained; c d (n), c m (n) are the sample points in the intercepted record block respectively, n is the sample point serial number in the intercepted record block. To ensure the stability of equation solving, it is usually required that n ≥ 4.

[0174] For the ghost wave adaptive subtraction processing based on the L2 norm, the filtering factor is determined by minimizing the sum of the squared error energy e in the following formula

[0175] e = ||c d - c m f||2 (Formula 13)

[0176] where the vector c d and c m represent the vectors composed of the data blocks of the original record and the multiple wave record respectively. The filtering factor f can be obtained by solving Equation (13). When the length of the filtering factor f is 1, the solution of the linear equation system can be avoided, and the filtering factor can be directly obtained, that is

[0177]

[0178] In the formula, T represents the transpose of the matrix (or vector).

[0179] Substituting the filtering factor f obtained based on Formula (14) into Formula (12), the curvelet coefficient record block c p (n) after suppressing the residual multiple wave can be obtained. Assign its central sample point to the new curvelet coefficient c p(j, l, k) is the central position of the original rectangular window in the medium. For c p Perform the suppression processing of steps (12)-(14) on each sample point in (j, l, k), and the seismic record curve wave coefficient c after attenuating the residual multiple waves can be obtained p (j, l, k), and then obtain the suppression result p(x, t) of the residual multiple waves in the space-time domain through the inverse curve wave transform, that is

[0180]

[0181] In the formula, p(x, t) is the seismic record of the attenuated residual multiple waves; Represents the curve wave basis function; c p (j, l, k) is the curve wave coefficient of the attenuated residual multiple waves, and j, l, k represent the scale, direction, and position of the curve wave respectively; x and t represent the position of the seismic trace and the travel time of the reflection signal respectively.

[0182] The original shot gather record of the work area is as Figure 2 shown, and the constructed residual multiple wave shot gather record is as Figure 7 shown, both containing 6-cable data. To improve the efficiency of the curve wave transform, input the single-cable data in the two records, and transform it into the curve wave domain using formula (11); then, for the sample points at each scale, direction, and position in the curve wave domain, perform multiple wave matching attenuation based on the least square filtering process described by formula (12)-formula (14); finally, use formula (15) to transform the curve wave domain record with residual multiple waves eliminated back to the space-time domain. Figure 8 The comparison of the single-cable data of the residual multiple wave suppression effect is given. The left figure is the shot gather record before attenuating the residual multiple waves, and the right figure is the shot gather record after attenuating the residual multiple waves. Through comparison, it can be seen that after the attenuation processing of this method, the residual multiple waves in the original record are significantly suppressed, and the signal-to-noise ratio of the data is significantly improved.

[0183] The advantages and positive effects of the present invention are:

[0184] For exploration areas with multiple strong wave impedance interfaces, the seismic records usually contain a large number of multiple wave interferences of various types and strong amplitudes, completely covering the weak amplitude primary reflection signals in the deep. After the combined attenuation of multiple methods such as predictive deconvolution technology, multiple wave suppression methods based on apparent velocity differences, and multiple wave suppression methods based on wave theory, obvious multiple wave residuals still exist in the data. Through the construction and adaptive subtraction processing of the multiple wave records of the present invention method, the residual multiple waves in the data can be further eliminated, and the signal-to-noise ratio of the seismic data can be significantly improved. The present invention is a beneficial supplement to the conventional multiple wave attenuation methods and provides a new way for multiple wave suppression.

[0185] The above has described in detail an embodiment of the present invention, but the above content is only a preferred embodiment of the present invention and cannot be considered as defining the scope of implementation of the present invention. All equivalent changes and improvements made according to the scope of application of the present invention shall still fall within the scope covered by the patent of the present invention.

Claims

1. A method for constructing and attenuating residual multiple waves based on forward modeling of Rayleigh integral, characterized in that: It includes the following steps: S1: Establish a high-resolution prestack time migration velocity spectrum, pick up the multiple velocity, and establish a multiple migration velocity field; S2: Calculate the 3D prestack time migration based on the multiple migration velocity field; S3: Construct a residual multiple record based on the forward modeling of the time-domain Rayleigh integral; S4: Suppress the residual multiples based on curvelet-domain least squares filtering.

2. A method for constructing and attenuating residual multiple waves based on Rayleigh integral forward modeling according to claim 1, characterized in that: In the above S1, the establishment of the high-resolution prestack time migration velocity spectrum includes the following steps: S11: Determine the initial migration velocity v0, velocity interval Δv, and maximum velocity v according to the velocity range of multiples in seismic data max , then the series of velocities v k The formula is as follows v k = v0 + kΔv (Equation 1) In the formula, k is the velocity serial number, 0 ≤ k ≤ K, where K is the serial number corresponding to the maximum velocity, and the following relational expression is satisfied: S12: The prestack time migration velocity spectrum is the migration result of the series of velocities. Based on the series of velocities defined in formula (1), the single-shot record d i (j, t) is subjected to Kirchhoff prestack time migration in the shot domain, and the formula is as follows, Where, IG k is the imaging gather migrated using velocity v k ; P is the coordinate of the imaging point corresponding to the spatial position of the velocity spectrum and the travel time within the spectrum; W M is the weight factor of Kirchhoff migration, τ represents the sum of the travel times of the ray from the shot point to the imaging point and from the imaging point to the geophone, (S i , G j ) is the coordinate of a shot - geophone pair; d i is the single - shot record, i and j are the shot number and trace number respectively, t represents the travel time, and I and J are the number of shots and traces used during migration respectively; S13: The weight factor W M (v k , S i , G j , P) and the travel time τ(v k , S i , G j , P) both depend on the speed v k , and the two are expressed as follows, where \(r_0\) represents the spatial distance from the shot point \(S\) of the ray i to the imaging point \(P\), and \(r\) is the spatial distance from the imaging point \(P\) to the geophone point \(G\) j ; S14: Based on formula (3), introduce the in-phase weighting process into the velocity spectrum calculation. The calculation formula for the in-phase weighting factor of the imaging gather is as follows: where F(v k , P) represents the in-phase weighting factor with respect to the velocity v k ; i is the trace number in the imaging gather, I represents the total number of traces; l is the time window length; ζ is a constant to ensure that the denominator is not zero; S15: Each scanning speed value v k For the corresponding imaging gather, the data of each trace are horizontally stacked and the absolute value is taken to obtain the migration velocity spectrum. Using formula (5), the weighted factor is obtained to perform in-phase weighting on the above results. This process is expressed as where V(v k , P) represents the high-resolution migration velocity spectrum, P is the coordinate of the imaging point corresponding to the spatial position of the velocity spectrum and the travel time within the spectrum; i is the trace number in the imaging gather, and I represents the total number of traces.

3. A method for constructing and attenuating residual multiple waves based on Rayleigh integral forward modeling according to claim 2, characterized in that: In the above S1, the establishment of the multiple migration velocity field includes the following steps: S16: For the migration velocity spectrum at a certain spatial position, draw a two-dimensional image with the velocity value v on the abscissa and the travel time t on the ordinate. Identify the energy of the residual multiples in the velocity spectrum, and pick up the corresponding positions to obtain the velocity curve of the residual multiples; S17: Interpolate in time and space according to the velocity curves at different positions to create a multiple migration velocity field.

4. A method for constructing and attenuating residual multiple waves based on Rayleigh integral forward modeling according to any one of claims 1 to 3, characterized in that: The above S2 includes the following steps: S21: According to the multiple migration velocity field, use the prestack seismic data to generate a multiple imaging data volume. The Kirchhoff prestack migration in the shot gather domain is expressed as an integral summation process in the following form: In the formula, M represents the prestack time migration data volume of multiple waves, P is the sample point coordinate of the data volume; W M is the weight factor of Kirchhoff integral migration, τ represents the sum of the travel times of the ray from the shot point to the imaging point and from the imaging point to the geophone, (S i , G j ) is the coordinate of a shot - geophone pair; d i is the shot gather record containing residual multiple waves, i and j are the shot number and trace number respectively, t represents the travel time, and I and J are the number of shots and traces applied during migration; S22: Simplify the weight factor W in formula (7) M (S i ,G j ,P) and the calculation process of travel time τ(S i ,G j ,P), and both are expressed as In the formula: v(P) is the migration velocity of the multiples, cosθ represents the tilt factor, r0 is the distance from the shot point to the imaging point, and r is the distance from the imaging point to the shot point.

5. A method for constructing and attenuating residual multiple waves based on Rayleigh integral forward modeling according to claim 4, characterized in that: The above S3 includes the following steps: S31: Input the multiple data volume M(P) generated according to formula (7), introduce a direct-path ray tracing process based on the migration velocity model, and the forward modeling process of the time-domain Rayleigh integral is expressed as follows: where m is the multiple wave record constructed by forward modeling using the time-domain Rayleigh integral, (S i , G j ) is the coordinate of a shot-receiver pair, i and j represent the shot number and trace number respectively; M(P) represents the multiple wave imaging data volume, P is the profile sample point coordinate; ds represents the integration surface element; R is the weight factor of the time-domain Rayleigh integral, τ represents the sum of the travel times of the ray from the shot point S i to the point P and from the point P to the receiver point G j . S32: During the forward simulation of the time-domain Rayleigh integral shown in formula (9), the integral weight factor R(P, S i , G j ) uses precise diffraction superposition calculation, and the formula is as follows, In the formula, r0 and r are the propagation distances of the incident ray and the diffracted ray respectively; θ0 and θ are the angles between r0, r and the normal vector n of the bin; v(P) represents the migration velocity of the multiples, and P is the sample point coordinate of the data volume.

6. A method for constructing and attenuating residual multiple waves based on Rayleigh integral forward modeling according to claim 5, characterized in that: The above S4 includes the following steps: S41: Perform multiple-wave adaptive subtraction in the curvelet domain with sparse features, transform the original seismic record d(x,t) and the predicted multiple-wave record m(x,t) to the curvelet domain, and obtain the curvelet coefficients c d (j, l, k) and c m (j, l, k), and the formula is as follows In the formula, represents a curvelet basis function; c d (j, l, k) and c m (j, l, k) represent the curvelet coefficients of the original seismic record and the multiple wave record respectively, where j, l, and k represent the scale, direction, and position of the curvelet respectively; x and t represent the position of the seismic trace and the travel time of the reflection signal respectively; S42: For the sample point at position k in the two-dimensional curvelet coefficient matrix at any scale j and direction l, set a rectangular window centered on it, intercept the data block therein for least squares filtering, and the result of subtracting the multiples is expressed as follows: c p y(n) = c d y(n) - f * c m y(n) (Equation 12) where c p (n) is the recorded sample point after attenuating the multiple wave; f is the filter factor to be obtained; c d (n), c m (n) are the sample points in the intercepted record block respectively, and n is the sample point serial number in the intercepted record block; S43: For the ghost wave adaptive subtraction processing based on the L2 norm, determine the filtering factor by minimizing the sum of squared error energy e in the following formula. The formula is as follows: e = ||c d -c m f||2 (Formula 13) Among them, vector c d and c m represent the vectors composed of the original record and the multiple record data blocks respectively. The filtering factor f is obtained by solving formula (13). When the length of the filtering factor f is 1, the solution of the linear equation system can be avoided, and the filtering factor can be directly obtained. The formula is as follows In the formula, T represents the transpose of the matrix (or vector); S44: Substitute the filtering factor f obtained based on formula (14) into formula (12) to obtain the curvelet coefficient record block c after suppressing the residual multiple waves. p (n), and assign its central sample to the new curvelet coefficient c p (j, l, k) at the central position of the original rectangular window in, and perform the suppression processing of S42 - S43 on each sample in c p (j, l, k) to obtain the curvelet coefficient c of the seismic record after attenuating the residual multiple waves p (j, l, k), and then obtain the suppression result p(x, t) of the residual multiple waves in the space-time domain through the inverse curvelet transform. The formula is as follows. Wherein, p(x,t) is the seismic record of the attenuated residual multiple wave; represents the curvelet basis function; c p (j, l, k) are the curvelet coefficients of the attenuated residual multiple wave, and j, l, and k respectively represent the scale, direction, and position of the curvelet; x and t respectively represent the position of the seismic trace and the travel time of the reflected signal.

7. A method for constructing and attenuating residual multiple waves based on Rayleigh integral forward modeling according to claim 6, characterized in that: It further includes the following steps: Input the single-cable data in the original shot gather record and the single-cable data in the residual multiple shot gather record into formula (11) respectively, and use formula (11) to transform them into the curvelet domain; For the sample points at each scale, direction and position in the curvelet domain, perform multiple matching attenuation based on the least squares filtering process described by formula (12) - formula (14); Use formula (15) to transform the curvelet domain record with residual multiples eliminated back to the spatio-temporal domain.

8. A device, characterized in that: Run the data processing method according to any one of claims 1 to 7.

9. An apparatus, comprising a memory, a processor, and an algorithm stored in the memory and executable on the processor, characterized in that: When the processor executes the computer program, it implements the data processing method according to any one of claims 1 to 7.

10. A computer-readable storage medium storing a computer algorithm, characterized in that, When the computer algorithm is executed by the processor, it implements the data processing according to any one of claims 1 to 7.

Images