Backward wave field synthesis processing method and device, computer equipment and storage medium

Abstract

The invention relates to a return wave field synthesis processing method and device, computer equipment and a storage medium. The method comprises the following steps: respectively acquiring seismic source return wave data and seismic source wavelet data; determining a sparse linear filter coefficient according to the seismic source return wave data and the seismic source wavelet data; and performing weighted stacking on the forward transmission wave field time delay by using the numerical value and delay of the coefficient of the non-zero sparse linear filter, and synthesizing a backward transmission wave field. According to the return wave field synthesis processing method, the seismic source return data is used as the convolution of the seismic source wavelet data and the sparse linear filter, the coefficient of the sparse linear filter is determined through the seismic source return data and the seismic source wavelet data, and then according to the superposition principle, weighted superposition is performed on the forward wave field time delay to synthesize the return wave field. The return wave field can be obtained without relying on the return data to solve the wave equation, the computing power is greatly reduced, the processing efficiency can be improved, the processing time can be shortened, the memory consumption can be reduced, and thus the computing power cost can be reduced.

Description

Technical Field

[0001] The present application relates to the technical field of earthquake source data return data processing, and in particular to a return wave field synthesis processing method, device, computer equipment and storage medium. Background Art

[0002] In the field of velocity modeling and seismic imaging technology in oil and gas exploration and development, the analysis and processing of source return waves are usually involved, which are used as the basic data for exploration and analysis. Source return waves refer to the signals generated by the seismic source in seismic exploration, which propagate underground and then return to the surface and are received by the seismic detectors. This signal contains information about the underground geological structure and is one of the basic data for seismic exploration. The characteristics of source return waves are closely related to factors such as the type of source, the layout of seismic detectors, and the underground geological structure. In seismic exploration, the reflection coefficient of the underground stratum can be obtained by collecting and processing the source return waves, and then the properties and structure of the stratum can be analyzed, providing important data support for the exploration of mineral resources such as oil and natural gas.

[0003] Geophysical methods based on wave equations, such as reverse time migration and full waveform inversion, often need to calculate the forward wave field of the earthquake source and the return wave field of the seismic data. The calculation of conventional return wave fields is to use the return seismic data as the earthquake source and use differential equation solving methods such as finite difference method and finite element method in the discrete grid of the underground physical parameter model to solve the wave equation. In the actual application of geophysical methods based on pre-stack shot records, this conventional return wave field calculation method has two disadvantages: 1. It is necessary to solve the wave field of all model grid points within the calculation range at the same time, so a large amount of memory is required; 2. The complexity of return wave field calculation increases sharply with the increase of model dimension and complexity of wave equation, thus hindering the practical application of more advanced geophysical methods. In other words, the existing synthetic processing methods for return waves have the problem of large amounts of computational processing data, which requires solving the wave fields of all model grid points within the calculation range and occupies a large amount of memory. The processing efficiency is low and the processing time is long. In addition, the processing process is complex and increases dramatically with the increase of model dimensions and the complexity of the wave equation. Therefore, it hinders the actual application of more advanced geophysical methods and limits the application scenarios. Summary of the invention

[0004] Based on this, it is necessary to provide a return wave field synthesis processing method, device, computer equipment and storage medium that can improve processing efficiency, reduce processing time, reduce memory consumption, reduce computing power costs and improve universality.

[0005] In a first aspect, the present application provides a return wave field synthesis processing method, comprising the following steps:

[0006] Acquire the earthquake source return wave data and earthquake source wavelet data respectively;

[0007] Determining sparse linear filter coefficients according to the earthquake source return wave data and the earthquake source wavelet data;

[0008] The forward wavefield is time-delayed and weightedly superimposed using the non-zero values ​​and delays of the sparse linear filter coefficients to synthesize the return wavefield.

[0009] In one embodiment, the step of determining the sparse linear filter coefficients according to the earthquake source return wave data and the earthquake source wavelet data comprises:

[0010] According to the earthquake source return wave data and the earthquake source wavelet data, sparse linear filter coefficients are determined based on least squares matching.

[0011] In one embodiment, in the step of determining the sparse linear filter coefficients based on least squares matching, the sparse linear filter coefficients are determined using the following formula:

[0012]

[0013] Among them, φ(w) is the least square matching error, w(i) is the filter coefficient, f(τ(i)) is the source wavelet data after a time delay of τ(i), and d is the source return wave data.

[0014] In one of the embodiments, in the least squares matching formula, the effective set method is used to solve the sparse linear filter coefficients.

[0015] In one embodiment, in the step of using the non-zero sparse linear filter coefficient value and delay to perform time-delay weighted superposition on the forward wavefield to synthesize the return wavefield, the return wavefield is determined using the following formula:

[0016]

[0017] in, is the earthquake source position x s Excitation, at detector position x r The wave field value at the model space position x at time t transmitted by the received source record; U f (x,x′ s =x r , t-τ(i)) is the earthquake source at x r The wave field value of the forward wave field excited at the spatial position x after a time delay of t-τ(i).

[0018] In a second aspect, the present application further provides a return wave field synthesis processing device, comprising:

[0019] An acquisition module, used to respectively acquire earthquake source return wave data and earthquake source wavelet data;

[0020] A filter coefficient determination module, used to determine the sparse linear filter coefficient according to the earthquake source return wave data and the earthquake source wavelet data acquired by the acquisition module;

[0021] The synthesis module is used to use the values ​​and delays of the non-zero sparse linear filter coefficients given by the filter coefficient determination module to perform weighted superposition on the forward wave field delay to synthesize the return wave field.

[0022] In one of the embodiments, the filter coefficient determination module is used to determine the sparse linear filter coefficients based on the least squares matching according to the earthquake source return wave data and the earthquake source wavelet data.

[0023] In one embodiment, the filter coefficient determination module determines the sparse linear filter coefficient using the following formula:

[0024]

[0025] Among them, φ(w) is the least square matching error, w(i) is the filter coefficient, f(τ(i)) is the source wavelet data after a time delay of τ(i), and d is the source return wave data.

[0026] In a third aspect, the present application further provides a server, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein when the processor executes the program, the steps of the return wave field synthesis processing method described in any of the above embodiments are implemented.

[0027] In a fourth aspect, the present application further provides a computer-readable storage medium having a computer program stored thereon, and when the computer program is executed by a processor, the steps of the return wave field synthesis processing method described in any of the above embodiments are implemented.

[0028] The above-mentioned return wave field synthesis processing method is based on the convolution of the source return data as the source wavelet data and the sparse linear filter. The sparse linear filter coefficients are determined by the source return data and the source wavelet data. Then, according to the superposition principle, the forward wave field is time-delayed and weighted superimposed according to the delay and value of the non-zero coefficient of the sparse linear filter to synthesize the return wave field. In the traditional return wave field synthesis method, it is necessary to rely on the return data to solve the wave equation to obtain the return wave field. Compared with the traditional return wavelength synthesis processing method, the present application does not rely on solving the wave equation during the return, and can directly synthesize the return wave field value at any spatial position and any time. The computational complexity of this method is only related to the number of non-zero coefficients of the sparse filter to be solved, and has nothing to do with the model dimension, boundary conditions, grid order, and complexity of the wave equation. Therefore, when synthesizing the return wave field of a three-dimensional high-order complex wave equation, this method has a significant advantage of low computational cost compared to traditional methods. The backhaul wave field synthesis processing method of the present application does not need to rely on the backhaul data to solve the wave equation to obtain the backhaul wave field. The method of determining the backhaul wavelength through the sparse linear filter coefficient and the forward wavelength can greatly reduce the computing power, improve the processing efficiency, reduce the processing time, and reduce the memory consumption, thereby reducing the computing power cost. Moreover, the computational complexity is only related to the number of non-zero coefficients of the sparse filter to be solved, but is not related to the model dimension, boundary conditions, grid order, and the complexity of the wave equation. Therefore, it is more universal and can adapt to more application scenarios. BRIEF DESCRIPTION OF THE DRAWINGS

[0029] Figure 1 It is a flow chart of a processing method of a return wave field synthesis processing method in an embodiment;

[0030] Figure 2a An amplitude diagram of returned seismic data in a returned wave field synthesis processing method in an embodiment;

[0031] Figure 2b A diagram showing the amplitude of sparse filter coefficients in a backhaul wave field synthesis processing method in an embodiment;

[0032] Figure 2c A diagram showing the amplitude comparison between the comparative seismic data and the source wavelet subjected to time-delay superposition using a filter in a return wave field synthesis processing method in an embodiment;

[0033] Figure 3a , Figure 3c , Figure 3e It is a synthetic wavefield diagram applied to shot gather data testing in a return wavefield synthesis processing method in an embodiment;

[0034] Figure 3b , Figure 3d , Figure 3fIt is a synthetic wavefield diagram used in the traditional return wavefield synthesis method for shot gather data testing;

[0035] Figure 4 It is a structural schematic diagram of a processing device for a return wave field synthesis processing method in an embodiment;

[0036] Figure 5 Schematic diagram of the structure of a computer device in one embodiment. DETAILED DESCRIPTION

[0037] In order to facilitate the understanding of the present application, in order to make the above-mentioned purposes, features and advantages of the present application more obvious and easy to understand, the specific implementation methods of the present application are described in detail below in conjunction with the accompanying drawings. In the following description, many specific details are set forth to facilitate a full understanding of the present application, and the preferred implementation methods of the present application are given in the accompanying drawings. However, the present application can be implemented in many different forms and is not limited to the implementation methods described herein. On the contrary, the purpose of providing these implementation methods is to make the disclosure of the present application more thoroughly and comprehensively understood. The present application can be implemented in many other ways different from those described herein, and those skilled in the art can make similar improvements without violating the connotation of the present application, so the present application is not limited to the specific embodiments disclosed below. Unless otherwise defined, all technical and scientific terms used herein have the same meaning as those generally understood by those skilled in the art of the present application. The terms used herein are only for the purpose of describing specific implementation methods and are not intended to limit the present application. The term "and / or" used herein includes any and all combinations of one or more related listed items.

[0038] Example 1

[0039] In the first aspect, the present application provides a method for processing return wave field synthesis, see Figure 1 The return wave field synthesis processing method includes the following steps:

[0040] S210: Acquire earthquake source return wave data and earthquake source wavelet data respectively;

[0041] In this application, the source return wave data and the source wavelet data are both instrument received data. Optionally, the source wavelet refers to the seismic wave generated when the earthquake source is excited and propagates in the viscoelastic medium. Its waveform will change with time, and finally form a signal with a certain starting time, limited energy and a certain continuation length, which is the source wavelet. The return wave data refers to the data of the seismic wave returning to the ground after being reflected underground; when the seismic wave encounters a concave interface with a large interface depth or a situation with a very large velocity gradient, a rotation phenomenon will occur to form a rotation wave, and the time-distance curve of the rotation wave is overlapping and cross-shaped. This is because the reflected waves on both sides of the concave interface and the flat interface reach the same point on the ground, and the order of the points on the time-distance curve is opposite to the order of the corresponding reflection points; the characteristic of the rotation wave is that the apparent velocity is small and the amplitude increases at the turning point. It is always accompanied by normal reflection waves and interferes with each other. Optionally, the source wavelet data and the return wave data can be received by a seismograph, which is a highly sensitive seismic wave receiving device that can record the vibration of seismic waves and convert them into analyzable data. The source wavelet data can record the start time and energy of the seismic waves, as well as the changes in the waveform, through the seismograph. The return wave data can also be recorded by the seismograph to record the signal of the seismic wave after it is reflected underground, and obtain information such as its arrival time and amplitude. Of course, it is not limited to this. It should be noted that how to obtain the source return wave data and the source wavelet data based on the relevant instruments, please refer to the prior art, and this application will not go into details here.

[0042] S220: Determine sparse linear filter coefficients according to the earthquake source return wave data and the earthquake source wavelet data;

[0043] In this embodiment, the seismic return data is equal to the convolution of the source wavelet data and the sparse linear filter, and thus the sparse linear filter coefficients are determined based on the source return wave data and the source wavelet data.

[0044] Specifically, the earthquake return data d and the source wavelet data f satisfy the following relationship:

[0045] d=wf

[0046] Where w is the sparse linear filter coefficient.

[0047] In order to better obtain the sparse linear filter coefficients, in one embodiment, the step of determining the sparse linear filter coefficients according to the earthquake source return wave data and the earthquake source wavelet data includes:

[0048] According to the earthquake source return wave data and the earthquake source wavelet data, sparse linear filter coefficients are determined based on least squares matching.

[0049] It should be noted that in the traditional synthesis of return wavelengths, it is necessary to solve the wave equation through seismic return data to obtain seismic return waves, and to solve the wave equation through the source wavelet data to obtain the source wavelet. In particular, in the process of solving the wave equation through seismic return data to obtain seismic return waves, geophysical methods based on wave equations, such as reverse time migration and full waveform inversion, often need to calculate the forward wave field of the source and the return wave field of seismic data. The calculation of conventional return wave fields is to use the return seismic data as the source, and use differential equation solving methods such as finite difference method and finite element method in the discrete grid of the underground physical parameter model to solve the wave equation. In the actual application of geophysical methods based on pre-stack shot records, this conventional return wave field calculation method has two disadvantages: 1. It is necessary to solve the wave field of all model grid points within the calculation range at the same time, so a large amount of memory is required; 2. The complexity of return wave field calculation increases sharply with the increase of model dimension and complexity of wave equation, thus hindering the practical application of more advanced geophysical methods. In this embodiment, in order to reduce the amount of calculation for solving complex wave equations with return wave data, by avoiding the strategy of solving wave equations, by combining seismic return data equal to the convolution of source wavelet data and sparse linear filter, the sparse linear filter coefficients are obtained, and then the return wave field is obtained by convolution of the forward wave field of the source wavelet field and the sparse linear filter coefficients. This method is based on the principle of energy superposition, and uses the time delay superposition of the forward wave field to synthesize the return wave field for the return seismic data near the source. This method does not need to solve the wave equation when synthesizing the return wave field, and can synthesize the return wave field at any time at any position in the model, and its computational complexity does not increase with the increase of the model dimension and the complexity of the wave equation. The return wave field synthesized by the present invention does not rely on solving the wave equation during the return, and can directly synthesize the return wave field value at any spatial position and any time. The return wave field calculated by this method has high least squares accuracy. The computational complexity of this method is only related to the number of non-zero coefficients of the sparse filter being solved, and has nothing to do with the model dimension, boundary conditions, grid order, and complexity of the wave equation. Therefore, when synthesizing the return wave field of a three-dimensional high-order complex wave equation, this method has a significant advantage of low computational cost over traditional methods.

[0050] In one embodiment, in the step of determining the sparse linear filter coefficients based on least squares matching, the sparse linear filter coefficients are determined using the following formula:

[0051]

[0052] Among them, φ(w) is the least square matching error, w(i) is the sparse linear filter coefficient, f(τ(i)) is the source wavelet data after time delay τ(i), and d is the source return wave data.

[0053] Optionally, least squares matching has an iteration number. When using least squares matching to solve a system of linear equations, there is a convergence precision, which is a key factor in determining the number of iterations. When solving matrix decomposition using alternating least squares matching, the maximum number of iterations is determined by the dimension of the score matrix and the degree of the coefficients of the score matrix. Specifically, if the dimension of the original matrix is ​​low or the degree of its coefficients is low (more sparse), then fewer iterations are required; conversely, if the dimension of the original matrix is ​​high or the degree of its coefficients is high (more dense), then more iterations are required.

[0054] In the present application, the least square matching error and the number of iterations can be defined to obtain more appropriate sparse linear filter coefficients.

[0055] In one of the embodiments, in the least squares matching formula, the effective set method is used to solve the sparse linear filter coefficients.

[0056] In a specific embodiment, the process of solving the sparse linear filter coefficients using the effective set method is as follows:

[0057] (1) Input the source return wave data d, source wavelet data f, number of iterations n, and least squares matching error φ(w);

[0058] (2) Define the search vector e, which satisfies the following relationship:

[0059] e=F T dF T FWf

[0060] Where F is the Toblitz matrix of the source wavelet data, T is the transposed matrix, and w is the filter coefficient;

[0061] (3) When the number of iterations does not reach the maximum number of iterations, and the maximum value of the e vector is greater than the least squares matching error φ(w), define:

[0062]

[0063] Move the element number of m in vector e out of R and add it to P,

[0064] s P =[(F P ) T F P ] -1 (F P ) T d

[0065] Where, R is the number of sampling points, R = {1, ..., ns}, ns is the number of seismic data sample points;

[0066] In this embodiment, when the maximum value of the e vector is greater than the least squares matching error φ(w), the element number of m in the vector e is moved out of R and added to P, so as to continuously adjust s. P In this embodiment, F P is the Toblitz matrix F of the source wavelet based on the specific sampling point data in the corresponding sampling point R.

[0067] (4) When s P When it is less than 0, it is defined as:

[0068] α=-min n∈P [w n / (w n -s n )]

[0069] w=w+α(sw)

[0070] When w is less than or equal to 0, move the element numbers less than 0 in w out of P and add them to R. Define:

[0071] s P =[(F P ) T F P ] -1 (F P ) T d

[0072] s R =0

[0073] In this embodiment, s P and R , respectively corresponding to the corresponding sampling point data.

[0074] definition:

[0075] w=s

[0076] e=F T dF T FWf

[0077] When the maximum value of the e vector is less than or equal to the least square matching error φ(w), the iteration is terminated, otherwise the loop continues and returns to execute step (3).

[0078] In a more specific embodiment, the process algorithm for solving the coefficients of the sparse linear filter using the effective set method is as follows:

[0079] Input: F, d, n, err

[0080]

[0081] R={1,...,ns}, ns is the number of seismic data samples

[0082] w=0

[0083] s=0

[0084] e=F T dF T FWf

[0085] iter=1

[0086] while md max i∈R (e i )>errand iter<ndo:

[0087]

[0088] Move the element number of m in vector e out of R and add it to P

[0089] s P =[(F P ) T F P ] -1 (F P ) T d

[0090] whilemin(s P )≤0do:

[0091] α=-min n∈P [w n / (w n -s n )]

[0092] w=w+α(sw)

[0093] Move the element numbers less than 0 in w out of P and add them to R

[0094] s P =[(F P ) T F P ] -1 (F P ) T d

[0095] s R =0

[0096] w=s

[0097] e=F T dF T FWf

[0098] iter=iter+1

[0099] return w

[0100] In this algorithm, F is the Toblitz matrix of the source wavelet data, n is the maximum number of algorithm iterations, and w also defines the upper limit of the number of zero coefficients. err defines the loop termination matching error threshold. P ,s R are subvectors of vector s, which are the elements in vector s corresponding to the element numbers in P and R.

[0101] S230: Using the non-zero values ​​and delays of the sparse linear filter coefficients, weightedly superimpose the forward wavefield delays to synthesize the return wavefield.

[0102] In this embodiment, the forward wavefield is the source wavelet field determined by the source wavelet data, that is, the forward wavefield can also be understood as the source wavelet field. In a specific embodiment, the forward wavefield can be determined by solving the wave equation by the source wavelet data, and the determination process of the forward wavefield can refer to the prior art and is not limited to that mentioned in this application. Specifically, the forward wavefield delay of the synthetic return wavefield is determined according to the delay of the non-zero coefficient of the sparse linear filter, and the superposition coefficient is determined according to the value of the filter coefficient.

[0103] In the present application, the numerical value and delay of the non-zero sparse linear filter coefficient are used to weight the time delay of the forward wave field and synthesize the return wave field, so as to reduce the calculation amount of the return wave data to solve the complex wave equation, avoid the strategy of solving the wave equation, and obtain the sparse linear filter coefficient by combining the seismic return data equal to the convolution of the source wavelet data and the sparse linear filter, and then obtain the return wave field by the convolution of the forward wave field of the source wavelet and the sparse linear filter coefficient. This method is based on the principle of energy superposition, and uses the time delay superposition of the forward wave field to synthesize the return wave field for the return seismic data near the source. This method does not need to solve the wave equation when synthesizing the return wave field, and can synthesize the return wave field at any time at any position in the model, and its computational complexity does not increase with the increase of the model dimension and the complexity of the wave equation. The return wave field synthesized by the present invention does not rely on solving the wave equation during the return, and can directly synthesize the return wave field value at any spatial position and any time. The least squares accuracy of the return wave field calculated by this method is high. The computational complexity of this method is only related to the number of non-zero coefficients of the sparse filter being solved, and has nothing to do with the model dimension, boundary conditions, grid order, and complexity of the wave equation. Therefore, when synthesizing the return wave field of a three-dimensional high-order complex wave equation, this method has a significant advantage of low computational cost over traditional methods.

[0104] In one embodiment, in the step of using the non-zero sparse linear filter coefficient value and delay to perform time-delay weighted superposition on the forward wavefield to synthesize the return wavefield, the return wavefield is determined using the following formula:

[0105]

[0106] in, is the earthquake source position x s Excitation, at detector position x r The wave field value at the model space position x at time t transmitted by the received source record; U f (x,x′ s =x r , t-τ(i)) is the earthquake source at x r The wave field value of the forward wave field excited at the spatial position x after the time delay t-τ(i). In this way, the return wave field can be determined better.

[0107] The above-mentioned return wave field synthesis processing method is based on the convolution of the source return data as the source wavelet data and the sparse linear filter. The sparse linear filter coefficients are determined by the source return data and the source wavelet data. Then, according to the superposition principle, the forward wave field is time-delayed and weighted superimposed according to the delay and value of the non-zero coefficient of the sparse linear filter to synthesize the return wave field. In the traditional return wave field synthesis method, it is necessary to rely on the return data to solve the wave equation to obtain the return wave field. Compared with the traditional return wavelength synthesis processing method, the present application does not rely on solving the wave equation during the return, and can directly synthesize the return wave field value at any spatial position and any time. The computational complexity of this method is only related to the number of non-zero coefficients of the sparse filter to be solved, and has nothing to do with the model dimension, boundary conditions, grid order, and complexity of the wave equation. Therefore, when synthesizing the return wave field of a three-dimensional high-order complex wave equation, this method has a significant advantage of low computational cost compared to traditional methods. The backhaul wave field synthesis processing method of the present application does not need to rely on the backhaul data to solve the wave equation to obtain the backhaul wave field. The method of determining the backhaul wavelength through the sparse linear filter coefficient and the forward wavelength can greatly reduce the computing power, improve the processing efficiency, reduce the processing time, and reduce the memory consumption, thereby reducing the computing power cost. Moreover, the computational complexity is only related to the number of non-zero coefficients of the sparse filter to be solved, but is not related to the model dimension, boundary conditions, grid order, and the complexity of the wave equation. Therefore, it is more universal and can adapt to more application scenarios.

[0108] Example 2

[0109] This embodiment provides a specific application process of the return wave field synthesis processing method of the present application. Based on a two-dimensional complex underground physical model, the method of the present invention is tested using shot gather data to compare the accuracy and computational cost of the return wave field synthesized by the present method with the return wave field calculated by the traditional finite difference (using second-order time and fourth-order space grid) method of solving the wave equation. Please refer to Figure 2a , Figure 2b and Figure 2c , Figure 2a To transmit seismic data back, Figure 2b is the sparse filter coefficient calculated by the effective set method in the return wave field synthesis processing method of this application, Figure 2c For this application, compare seismic data with source wavelets time-delayed and stacked using sparse filter coefficients, see Figure 3a to Figure 3f , on the left Figure 3a , Figure 3c , Figure 3e For this application, the wave field synthesis processing method returns the wave field, the left Figure 3b , Figure 3d , Figure 3f The wave field calculated by the traditional finite difference method (using second-order time and fourth-order space grid) is Figure 3a to Figure 3f Compare the wave field snapshots at different propagation times. The results show that the return wave field synthesis processing method provided by the present application has a 99.5% match between the results of the method of the present invention and the results of the traditional method, but the calculation cost of the method of the present invention is 50% of that of the traditional method. Therefore, the computing power cost can be reduced, and the computational complexity is only related to the number of non-zero coefficients of the sparse filter to be solved, and has nothing to do with the model dimension, boundary conditions, grid order, and complexity of the wave equation, so it is more universal and can adapt to more application scenarios. The return wave field synthesis processing method of the present application, because there is no need to rely on the return data to solve the wave equation to obtain the return wave field, the method of determining the return wavelength by the sparse linear filter coefficient and the forward wavelength, greatly reduces the computing power, can improve processing efficiency, reduce processing time, and can also reduce memory consumption, so it can reduce the computing power cost, and the computational complexity is only related to the number of non-zero coefficients of the sparse filter to be solved, and has nothing to do with the model dimension, boundary conditions, grid order, and complexity of the wave equation, so it is more universal and can adapt to more application scenarios.

[0110] Example 3

[0111] In a second aspect, the present application also provides a return wave field synthesis processing device, see Figure 4 , the return wave field synthesis processing device includes an acquisition module, a filter coefficient determination module and a synthesis module connected in sequence;

[0112] The acquisition module is used to respectively acquire the earthquake source return wave data and the earthquake source wavelet data;

[0113] The filter coefficient determination module is used to determine the sparse linear filter coefficient according to the earthquake source return wave data and the earthquake source wavelet data acquired by the acquisition module;

[0114] The synthesis module is used to use the values ​​and delays of the non-zero sparse linear filter coefficients given by the filter coefficient determination module to perform weighted superposition on the forward wave field delay to synthesize the return wave field.

[0115] The above-mentioned return wave field synthesis processing device adopts the above-mentioned acquisition module, filter coefficient determination module and synthesis module, based on the source return data as the convolution of the source wavelet data and the sparse linear filter, the sparse linear filter coefficient is determined by the source return data and the source wavelet data, and then according to the superposition principle, according to the delay and value of the non-zero coefficient of the sparse linear filter, the forward wave field delay is weighted and superimposed to synthesize the return wave field. In the traditional return wave field synthesis method, it is necessary to rely on the return data to solve the wave equation to obtain the return wave field. Compared with the traditional return wavelength synthesis processing method, the present application does not rely on solving the wave equation during the return, and can directly synthesize the return wave field value at any spatial position and any time. The computational complexity of this method is only related to the number of non-zero coefficients of the sparse filter to be solved, and has nothing to do with the model dimension, boundary conditions, grid order, and complexity of the wave equation. Therefore, when synthesizing the return wave field of a three-dimensional high-order complex wave equation, this method has a significant advantage of low computational cost compared with the traditional method. The return wave field synthesis processing device of the present application does not need to rely on the return data to solve the wave equation to obtain the return wave field. The method of determining the return wavelength through the sparse linear filter coefficient and the forward wavelength can greatly reduce the computing power, improve the processing efficiency, reduce the processing time, and reduce the memory consumption, thereby reducing the computing power cost. Moreover, the calculation complexity is only related to the number of non-zero coefficients of the sparse filter to be solved, and is not related to the model dimension, boundary conditions, grid order, and the complexity of the wave equation. Therefore, it is more universal and can adapt to more application scenarios.

[0116] In one embodiment, the acquisition module includes a return wave data acquisition unit and a source wavelet data acquisition unit, the return wave data acquisition unit is used to acquire source return wave data, and the source wavelet data acquisition unit is used to acquire source wavelet data.

[0117] In one of the embodiments, the filter coefficient determination module is used to determine the sparse linear filter coefficients based on the least squares matching according to the earthquake source return wave data and the earthquake source wavelet data.

[0118] In one embodiment, the filter coefficient determination module is used to determine the sparse linear filter coefficients based on the least squares matching according to the earthquake source return wave data and the earthquake source wavelet data, wherein the earthquake return data d and the earthquake source wavelet data f satisfy the following relationship:

[0119] d=wf

[0120] In the above formula, w is the sparse linear filter coefficient.

[0121] In one embodiment, the filter coefficient determination module determines the sparse linear filter coefficient using the following formula:

[0122]

[0123] Among them, φ(w) is the least square matching error, w(i) is the filter coefficient, f(τ(i)) is the source wavelet data after a time delay of τ(i), and d is the source return wave data.

[0124] Optionally, least squares matching has an iteration number. When using least squares matching to solve a system of linear equations, there is a convergence precision, which is a key factor in determining the number of iterations. When solving matrix decomposition using alternating least squares matching, the maximum number of iterations is determined by the dimension of the score matrix and the degree of the coefficients of the score matrix. Specifically, if the dimension of the original matrix is ​​low or the degree of its coefficients is low (more sparse), then fewer iterations are required; conversely, if the dimension of the original matrix is ​​high or the degree of its coefficients is high (more dense), then more iterations are required.

[0125] In one of the embodiments, the filter coefficient determination module includes a formula determination unit and a filter coefficient solving unit;

[0126] The formula determination unit is used to establish the following formula based on the least squares matching according to the earthquake source return wave data and the earthquake source wavelet data:

[0127]

[0128] Among them, φ(w) is the least square matching error, w(i) is the sparse linear filter coefficient, f(τ(i)) is the source wavelet data after time delay τ(i), and d is the source return wave data.

[0129] The filter coefficient solving unit is used for solving the sparse linear filter coefficients according to the formula by adopting the effective set method.

[0130] In a specific embodiment, the filter coefficient solving unit is used to solve the sparse linear filter coefficients according to the formula using the effective set method, wherein the solving process is as follows:

[0131] (1) Input the source return wave data d, source wavelet data f, number of iterations n, and least squares matching error φ(w);

[0132] (2) Define the search vector e, which satisfies the following relationship:

[0133] e=F T dF T FWf

[0134] Where F is the Toblitz matrix of the source wavelet data, T is the transposed matrix, and w is the filter coefficient;

[0135] (3) When the number of iterations does not reach the maximum number of iterations, and the maximum value of the e vector is greater than the least squares matching error φ(w), define:

[0136]

[0137] Move the element number of m in vector e out of R and add it to P,

[0138] s P =[(F P ) T F P ] -1 (F P ) T d

[0139] Where, R is the number of sampling points, R = {1, ..., ns}, ns is the number of seismic data sample points;

[0140] In this embodiment, when the maximum value of the e vector is greater than the least squares matching error φ(w), the element number of m in the vector e is moved out of R and added to P, so as to continuously adjust s. P In this embodiment, F P is the Toblitz matrix F of the source wavelet based on the specific sampling point data in the corresponding sampling point R.

[0141] (4) When s P When it is less than 0, it is defined as:

[0142] α=-min n∈P [w n / (w n -s n )]

[0143] w=w+α(sw)

[0144] When w is less than or equal to 0, move the element numbers less than 0 in w out of P and add them to R. Define:

[0145] s P =[(F P ) T F P ] -1 (F P ) T d

[0146] s R =0

[0147] In this embodiment, s P and R , respectively corresponding to the corresponding sampling point data.

[0148] definition:

[0149] w=s

[0150] e=F T dF T FWf

[0151] When the maximum value of the e vector is less than or equal to the least square matching error φ(w), the iteration is terminated, otherwise the loop continues and returns to execute step (3).

[0152] In one embodiment, the synthesis module is used to use the non-zero sparse linear filter coefficient value and delay to perform time-delay weighted superposition on the forward wave field to synthesize the return wave field, wherein the return wave field is determined by the following formula:

[0153]

[0154] in, is the earthquake source position x s Excitation, at detector position x r The wave field value at the model space position x at time t transmitted by the received source record; U f (x,x′ s =x r , t-τ(i)) is the earthquake source at x r The wave field value of the forward wave field excited at the spatial position x after the time delay t-τ(i). In this way, the return wave field can be determined better.

[0155] The above-mentioned return wave field synthesis processing device is based on the convolution of the source return data and the sparse linear filter, and the sparse linear filter coefficients are determined by the source return data and the source wavelet data. Then, according to the superposition principle, the forward wave field is time-delayed and weighted superimposed according to the delay and value of the non-zero coefficient of the sparse linear filter to synthesize the return wave field. In the traditional return wave field synthesis method, it is necessary to rely on the return data to solve the wave equation to obtain the return wave field. Compared with the traditional return wavelength synthesis processing method, the present application does not rely on solving the wave equation during the return, and can directly synthesize the return wave field value at any spatial position and any time. The computational complexity of the return wave field synthesis processing device provided by the present application is only related to the number of non-zero coefficients of the sparse filter to be solved, and has nothing to do with the model dimension, boundary conditions, grid order, and complexity of the wave equation. Therefore, when synthesizing the return wave field of a three-dimensional high-order complex wave equation, this method has a significant advantage of low computational cost compared with the traditional method. The return wave field synthesis processing device of the present application does not need to rely on the return data to solve the wave equation to obtain the return wave field. The method of determining the return wavelength through the sparse linear filter coefficient and the forward wavelength can greatly reduce the computing power, improve the processing efficiency, reduce the processing time, and reduce the memory consumption, thereby reducing the computing power cost. Moreover, the calculation complexity is only related to the number of non-zero coefficients of the sparse filter to be solved, and is not related to the model dimension, boundary conditions, grid order, and the complexity of the wave equation. Therefore, it is more universal and can adapt to more application scenarios.

[0156] Example 4

[0157] In a third aspect, the present application further provides a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein when the processor executes the program, the steps of the return wave field synthesis processing method described in any of the above embodiments are implemented.

[0158] In one embodiment, the computer device may be Figure 5As shown. The computer device includes a processor, a memory, a network interface, a display screen and an input device connected through a system bus. Among them, the processor of the computer device is used to provide computing and control capabilities. The memory of the computer device includes a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system and a computer program. The internal memory provides an environment for the operation of the operating system and the computer program in the non-volatile storage medium. The network interface of the computer device is used to communicate with the server through a network connection. When the computer program is executed by the processor, a return wave field synthesis processing method as described in any of the above embodiments is implemented. The display screen of the computer device can be a liquid crystal display screen or an electronic ink display screen, and the input device of the computer device can be a touch layer covered on the display screen, or a key, trackball or touchpad set on the computer device housing, or an external keyboard, touchpad or mouse, etc.

[0159] Those skilled in the art will understand that Figure 5 The structure shown in the figure is only a block diagram of a part of the structure related to the solution of the present application, and does not constitute a limitation on the computer device to which the solution of the present application is applied. The specific computer device may include more or fewer components than those shown in the figure, or combine certain components, or have a different arrangement of components.

[0160] In one embodiment, when the computer program is executed by a processor, the following steps are also implemented:

[0161] Acquire the earthquake source return wave data and earthquake source wavelet data respectively;

[0162] Determining sparse linear filter coefficients according to the earthquake source return wave data and the earthquake source wavelet data;

[0163] The forward wavefield is time-delayed and weightedly superimposed using the non-zero values ​​and delays of the sparse linear filter coefficients to synthesize the return wavefield.

[0164] In one embodiment, when the computer program is executed by a processor, the following steps are also implemented:

[0165] The step of determining the sparse linear filter coefficients according to the earthquake source return wave data and the earthquake source wavelet data comprises:

[0166] According to the earthquake source return wave data and the earthquake source wavelet data, sparse linear filter coefficients are determined based on least squares matching.

[0167] In one embodiment, when the computer program is executed by a processor, the following steps are also implemented:

[0168] In the step of determining the sparse linear filter coefficients based on least squares matching, the sparse linear filter coefficients are determined using the following formula:

[0169]

[0170] Among them, φ(w) is the least square matching error, w(i) is the filter coefficient, f(τ(i)) is the source wavelet data after a time delay of τ(i), and d is the source return wave data.

[0171] In one embodiment, when the computer program is executed by a processor, the following steps are also implemented:

[0172] In the least squares matching formula, the effective set method is used to solve the sparse linear filter coefficients.

[0173] In one embodiment, when the computer program is executed by a processor, the following steps are also implemented:

[0174] In the step of using the non-zero values ​​and delays of the sparse linear filter coefficients to perform time-delay weighted superposition on the forward wavefield to synthesize the return wavefield, the return wavefield is determined using the following formula:

[0175]

[0176] in, is the earthquake source position x s Excitation, at detector position x r The wave field value at the model space position x at time t transmitted by the received source record; U f (x,x′ s =x r , t-τ(i)) is the earthquake source at x r The wave field value of the forward wave field excited at the spatial position x after a time delay of t-τ(i).

[0177] In one embodiment, when the computer program is executed by a processor, the following steps are also implemented:

[0178] The process of solving the coefficients of sparse linear filter using the effective set method is as follows:

[0179] (1) Input the source return wave data d, source wavelet data f, number of iterations n, and least squares matching error φ(w);

[0180] (2) Define the search vector e, which satisfies the following relationship:

[0181] e=F T dF T FWf

[0182] Where F is the Toblitz matrix of the source wavelet data, T is the transposed matrix, and w is the filter coefficient;

[0183] (3) When the number of iterations does not reach the maximum number of iterations, and the maximum value of the e vector is greater than the least squares matching error φ(w), define:

[0184]

[0185] Move the element number of m in vector e out of R and add it to P,

[0186] s P =[(F P ) T F P ] -1 (F P ) T d

[0187] Where, R is the number of sampling points, R = {1, ..., ns}, ns is the number of seismic data sample points;

[0188] In this embodiment, when the maximum value of the e vector is greater than the least squares matching error φ(w), the element number of m in the vector e is moved out of R and added to P, so as to continuously adjust s. P In this embodiment, F P is the Toblitz matrix F of the source wavelet based on the specific sampling point data in the corresponding sampling point R.

[0189] (4) When s P When it is less than 0, it is defined as:

[0190] α=-min n∈P [w n / (w n -s n )]

[0191] w=w+α(sw)

[0192] When w is less than or equal to 0, move the element numbers less than 0 in w out of P and add them to R. Define:

[0193] s P =[(F P ) T F P ] -1 (F P ) T d

[0194] s R =0

[0195] In this embodiment, s P and R , respectively corresponding to the corresponding sampling point data.

[0196] definition:

[0197] w=s

[0198] e=F T dF T FWf

[0199] When the maximum value of the e vector is less than or equal to the least square matching error φ(w), the iteration is terminated, otherwise the loop continues and returns to execute step (3).

[0200] Example 5

[0201] In a fourth aspect, the present application further provides a computer-readable storage medium having a computer program stored thereon, and when the computer program is executed by a processor, the steps of the return wave field synthesis processing method described in any of the above embodiments are implemented.

[0202] In one embodiment, when the computer program is executed by a processor, the following steps are also implemented:

[0203] Acquire the earthquake source return wave data and earthquake source wavelet data respectively;

[0204] Determining sparse linear filter coefficients according to the earthquake source return wave data and the earthquake source wavelet data;

[0205] The forward wavefield is time-delayed and weightedly superimposed using the non-zero values ​​and delays of the sparse linear filter coefficients to synthesize the return wavefield.

[0206] In one embodiment, when the computer program is executed by a processor, the following steps are also implemented:

[0207] The step of determining the sparse linear filter coefficients according to the earthquake source return wave data and the earthquake source wavelet data comprises:

[0208] According to the earthquake source return wave data and the earthquake source wavelet data, sparse linear filter coefficients are determined based on least squares matching.

[0209] In one embodiment, when the computer program is executed by a processor, the following steps are also implemented:

[0210] In the step of determining the sparse linear filter coefficients based on least squares matching, the sparse linear filter coefficients are determined using the following formula:

[0211]

[0212] Among them, φ(w) is the least square matching error, w(i) is the filter coefficient, f(τ(i)) is the source wavelet data after a time delay of τ(i), and d is the source return wave data.

[0213] In one embodiment, when the computer program is executed by a processor, the following steps are also implemented:

[0214] In the least squares matching formula, the effective set method is used to solve the sparse linear filter coefficients.

[0215] In one embodiment, when the computer program is executed by a processor, the following steps are also implemented:

[0216] In the step of using the non-zero values ​​and delays of the sparse linear filter coefficients to perform time-delay weighted superposition on the forward wavefield to synthesize the return wavefield, the return wavefield is determined using the following formula:

[0217]

[0218] in, is the earthquake source position x s Excitation, at detector position x r The wave field value at the model space position x at time t transmitted by the received source record; U f (x,x′ s =x r , t-τ(i)) is the earthquake source at x r The wave field value of the forward wave field excited at the spatial position x after a time delay of t-τ(i).

[0219] In one embodiment, when the computer program is executed by a processor, the following steps are also implemented:

[0220] The process of solving the coefficients of sparse linear filter using the effective set method is as follows:

[0221] (1) Input the source return wave data d, source wavelet data f, number of iterations n, and least squares matching error φ(w);

[0222] (2) Define the search vector e, which satisfies the following relationship:

[0223] e=F T dF T FWf

[0224] Where F is the Toblitz matrix of the source wavelet data, T is the transposed matrix, and w is the filter coefficient;

[0225] (3) When the number of iterations does not reach the maximum number of iterations, and the maximum value of the e vector is greater than the least squares matching error φ(w), define:

[0226]

[0227] Move the element number of m in vector e out of R and add it to P,

[0228] s P =[(F P ) T F P ] -1 (F P ) T d

[0229] Where, R is the number of sampling points, R = {1, ..., ns}, ns is the number of seismic data sample points;

[0230] In this embodiment, when the maximum value of the e vector is greater than the least squares matching error φ(w), the element number of m in the vector e is moved out of R and added to P, so as to continuously adjust s. P In this embodiment, F P is the Toblitz matrix F of the source wavelet based on the specific sampling point data in the corresponding sampling point R.

[0231] (4) When s P When it is less than 0, it is defined as:

[0232] α=-min n∈P [w n / (w n -s n )]

[0233] w=w+α(sw)

[0234] When w is less than or equal to 0, move the element numbers less than 0 in w out of P and add them to R. Define:

[0235] s P =[(F P ) T F P ] -1 (F P ) T d

[0236] s R =0

[0237] In this embodiment, s P and R , respectively corresponding to the corresponding sampling point data.

[0238] definition:

[0239] w=s

[0240] e=F T dF T FWf

[0241] When the maximum value of the e vector is less than or equal to the least square matching error φ(w), the iteration is terminated, otherwise the loop continues and returns to execute step (3).

[0242] Those skilled in the art can understand that all or part of the processes in the above-mentioned embodiment methods can be completed by instructing the relevant hardware through a computer program, and the computer program can be stored in a non-volatile computer-readable storage medium. When the computer program is executed, it can include the processes of the embodiments of the above-mentioned methods. Among them, any reference to memory, storage, database or other media used in the embodiments provided in this application can include non-volatile and / or volatile memory. Non-volatile memory can include read-only memory (ROM), programmable ROM (PROM), electrically programmable ROM (EPROM), electrically erasable programmable ROM (EEPROM) or flash memory. Volatile memory can include random access memory (RAM) or external cache memory. As an illustration and not limitation, RAM is available in many forms, such as static RAM (SRAM), dynamic RAM (DRAM), synchronous DRAM (SDRAM), double data rate SDRAM (DDRSDRAM), enhanced SDRAM (ESDRAM), synchronous link (Synchlink) DRAM (SLDRAM), memory bus (Rambus) direct RAM (RDRAM), direct memory bus dynamic RAM (DRDRAM), and memory bus dynamic RAM (RDRAM).

[0243] The above-mentioned return wave field synthesis processing method and device are based on the convolution of the source return data as the source wavelet data and the sparse linear filter. The sparse linear filter coefficients are determined by the source return data and the source wavelet data. Then, according to the superposition principle, the forward wave field is time-delayed and weighted superimposed according to the delay and value of the non-zero coefficient of the sparse linear filter to synthesize the return wave field. In the traditional return wave field synthesis method, it is necessary to rely on the return data to solve the wave equation to obtain the return wave field. Compared with the traditional return wavelength synthesis processing method, the present application does not rely on solving the wave equation during the return, and can directly synthesize the return wave field value at any spatial position and any time. The computational complexity of this method is only related to the number of non-zero coefficients of the sparse filter to be solved, and has nothing to do with the model dimension, boundary conditions, grid order, and complexity of the wave equation. Therefore, when synthesizing the return wave field of a three-dimensional high-order complex wave equation, this method has a significant advantage of low computational cost compared with traditional methods. The backhaul wave field synthesis processing method of the present application does not need to rely on the backhaul data to solve the wave equation to obtain the backhaul wave field. The method of determining the backhaul wavelength through the sparse linear filter coefficient and the forward wavelength can greatly reduce the computing power, improve the processing efficiency, reduce the processing time, and reduce the memory consumption, thereby reducing the computing power cost. Moreover, the computational complexity is only related to the number of non-zero coefficients of the sparse filter to be solved, but is not related to the model dimension, boundary conditions, grid order, and the complexity of the wave equation. Therefore, it is more universal and can adapt to more application scenarios.

[0244] The technical features of the above-described embodiments can be combined arbitrarily. In order to make the description concise, all possible combinations of the technical features in the above-described embodiments are not described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification. It should be noted that "in one embodiment", "for example", "for example", etc. in this application are intended to illustrate the present application, rather than to limit the present application. The above-described embodiments only express several implementation methods of the present application, and the description is relatively specific and detailed, but it cannot be understood as a limitation on the scope of the patent application. It should be pointed out that for a person of ordinary skill in the art, without departing from the concept of the present application, several variations and improvements can also be made, which all belong to the scope of protection of the present application. Therefore, the scope of protection of the patent in this application shall be subject to the attached claims.

Claims

1. A return wave field synthesis processing method, characterized in that: The steps include: Acquire the earthquake source return wave data and earthquake source wavelet data respectively; Determining sparse linear filter coefficients according to the earthquake source return wave data and the earthquake source wavelet data; The forward wavefield is time-delayed and weightedly superimposed using the non-zero values ​​and delays of the sparse linear filter coefficients to synthesize the return wavefield.

2. The return wave field synthesis processing method according to claim 1, characterized in that: The step of determining the sparse linear filter coefficients according to the earthquake source return wave data and the earthquake source wavelet data comprises: According to the earthquake source return wave data and the earthquake source wavelet data, sparse linear filter coefficients are determined based on least squares matching.

3. The return wave field synthesis processing method according to claim 2, characterized in that: In the step of determining the sparse linear filter coefficients based on least squares matching, the sparse linear filter coefficients are determined using the following formula: Among them, φ(w) is the least square matching error, w(i) is the filter coefficient, f(τ(i)) is the source wavelet data after a time delay of τ(i), and d is the source return wave data.

4. The return wave field synthesis processing method according to claim 3, characterized in that: In the least squares matching formula, the effective set method is used to solve the sparse linear filter coefficients.

5. The return wave field synthesis processing method according to claim 3, characterized in that: In the step of using the non-zero values ​​and delays of the sparse linear filter coefficients to perform time-delay weighted superposition on the forward wavefield to synthesize the return wavefield, the return wavefield is determined using the following formula: in, is the earthquake source position x s Excitation, at detector position x r The wave field value at the model space position x at time t transmitted by the received source record; U f (x,x′ s =x r , t-τ(i)) is the earthquake source at x r The wave field value of the forward wave field excited at the spatial position x after a time delay of t-τ(i).

6. A return wave field synthesis processing device, characterized in that: include: An acquisition module, used to respectively acquire earthquake source return wave data and earthquake source wavelet data; A filter coefficient determination module, used to determine the sparse linear filter coefficient according to the earthquake source return wave data and the earthquake source wavelet data acquired by the acquisition module; The synthesis module is used to use the values ​​and delays of the non-zero sparse linear filter coefficients given by the filter coefficient determination module to perform weighted superposition on the forward wave field delay to synthesize the return wave field.

7. The return wave field synthesis processing device according to claim 6, characterized in that: The filter coefficient determination module is used to determine the sparse linear filter coefficients based on the least squares matching according to the earthquake source return wave data and the earthquake source wavelet data.

8. The return wave field synthesis processing device according to claim 7, characterized in that: The filter coefficient determination module uses the following formula to determine the sparse linear filter coefficient: Among them, φ(w) is the least square matching error, w(i) is the filter coefficient, f(τ(i)) is the source wavelet data after a time delay of τ(i), and d is the source return wave data.

9. A computer device comprising a memory, a processor and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the program, the steps of the return wave field synthesis processing method described in any one of claims 1 to 5 are implemented.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps of the return wave field synthesis processing method according to any one of claims 1 to 5 are implemented.

Images