Back-propagating wave field synthesis processing method and device, computer device and storage medium

Abstract

The present application relates to a kind of back wave field synthesis processing method, device, computer equipment and storage medium, method includes the following steps: respectively obtaining source back wave data and source wavelet data;According to source back wave data and source wavelet data, determine sparse linear filter coefficient;Using the numerical value of non-zero sparse linear filter coefficient and delay, time delay weighted superposition of forward wave field, synthesis back wave field.The above-mentioned back wave field synthesis processing method, based on source back data is the convolution of source wavelet data and sparse linear filter, determine sparse linear filter coefficient by source back data and source wavelet data, then according to the principle of superposition, time delay weighted superposition of forward wave field, synthesis back wave field.No need to rely on back data to obtain back wave field by solving wave equation, power is greatly reduced, can improve processing efficiency, reduce processing time, can also reduce memory consumption, so as to reduce power cost.

Description

Technical Field

[0001] This application relates to the field of earthquake source data backhaul data processing technology, and in particular to a backhaul wavefield synthesis processing method, apparatus, computer equipment, and storage medium. Background Technology

[0002] In the fields of velocity modeling and seismic imaging technology in oil and gas exploration and development, the analysis and processing of source return waves are commonly involved, serving as fundamental data for exploration analysis. Source return waves refer to the signals received by seismic detectors after the seismic waves generated by the earthquake source propagate underground and return to the surface. This signal contains information about the underground geological structure and is one of the basic data in seismic exploration. The characteristics of source return waves are closely related to factors such as the type of earthquake source, the arrangement of seismic detectors, and the underground geological structure. In seismic exploration, by acquiring and processing source return waves, the reflection coefficients of underground strata can be obtained, allowing for the analysis of strata properties and structure, providing crucial 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 require calculating the forward propagation wavefield of the seismic source and the return propagation wavefield of the seismic data. Conventional return propagation wavefield calculations use the returned seismic data as the source and solve the wave equation using differential equation solving methods such as the finite difference method or the finite element method within a discrete grid of the subsurface physical parameter model. In practical applications of geophysical methods based on pre-stack shot gather records, this conventional return propagation wavefield calculation method has two drawbacks: 1. It requires solving the wavefield at all model grid points within the computational scope simultaneously, thus necessitating a large amount of memory; 2. The computational complexity of the return propagation wavefield increases dramatically with the increase in model dimensionality and the complexity of the wave equation, thus hindering the practical application of more advanced geophysical methods. In other words, existing methods for synthesizing backpropagation waves suffer from several problems. They involve processing massive amounts of data, requiring the solution of wave fields at all model grid points within the computational scope, which consumes a large amount of memory. This results in low processing efficiency and long processing time. Furthermore, the processing is quite complex, and the complexity increases dramatically with the increase in model dimension and wave equation complexity. Therefore, these issues hinder their practical application in more advanced geophysical methods and limit their application scenarios. Summary of the Invention

[0004] Therefore, it is necessary to provide a method, apparatus, computer equipment, and storage medium for backhaul wave field synthesis processing that can improve processing efficiency, reduce processing time, reduce memory consumption, lower computing costs, and improve versatility.

[0005] In a first aspect, this application provides a method for synthesizing and processing a return wave field, comprising the following steps:

[0006] Separately acquire source return wave data and source wavelet data;

[0007] The coefficients of the sparse linear filter are determined based on the source return wave data and the source wavelet data.

[0008] By using the numerical values ​​of non-zero sparse linear filter coefficients and their delays, the forward propagation wavefield is weighted and superimposed to synthesize the return propagation wavefield.

[0009] In one embodiment, the step of determining the sparse linear filter coefficients based on the source return wave data and the source wavelet data includes:

[0010] Based on the source return wave data and the source wavelet data, the coefficients of the sparse linear filter are determined using least squares matching.

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

[0012]

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

[0014] In one embodiment, the least squares matching formula uses the effective set method to solve for the sparse linear filter coefficients.

[0015] In one embodiment, the step of synthesizing the return wavefield by weighted superposition of the forward propagation wavefield using the values ​​of non-zero sparse linear filter coefficients and the delay is determined by the following formula:

[0016]

[0017] in, x is the location of the earthquake source. s Excitation, at detector position x r The received source record transmits the wavefield value at the model's spatial location x at time t; U f (x, x′) s =x r ,t-τ(i)) is the earthquake source at x r The wave field value at spatial location x after a time delay t-τ(i) for the propagating wave field excited at point x.

[0018] Secondly, this application also provides a backhaul wavefield synthesis processing apparatus, comprising:

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

[0020] The filter coefficient determination module is used to determine the sparse linear filter coefficients based on the source return wave data and the source wavelet data acquired by the acquisition module.

[0021] The synthesis module is used to combine the values ​​and delays of the non-zero sparse linear filter coefficients given by the filter coefficient determination module, and to weight and superimpose the forward propagation wave field to synthesize the return propagation wave field.

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

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

[0024]

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

[0026] Thirdly, this application also provides a server, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the steps of the backhaul wavefield synthesis processing method as described in any of the above embodiments.

[0027] Fourthly, this application also provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the backhaul wave field synthesis processing method as described in any of the above embodiments.

[0028] The aforementioned method for synthesizing the returned wavefield is based on the assumption that the source-returned data is the convolution of source wavelet data and a sparse linear filter. The coefficients of the sparse linear filter are determined using the source-returned data and the source wavelet data. Then, according to the superposition principle, the forward-propagating wavefield is weighted and superimposed based on the delay and value of the non-zero coefficients of the sparse linear filter to synthesize the returned wavefield. In contrast, traditional methods for synthesizing returned wavefields require solving the wave equation using the returned data. This method, unlike traditional methods for synthesizing returned wavelengths, does not rely on solving the wave equation during the return process and can directly synthesize returned wavefield values ​​at any spatial location and time. The computational complexity of this method depends only on the number of non-zero coefficients of the solved sparse filter and is independent of the model dimension, boundary conditions, grid order, and the complexity of the wave equation. Therefore, when synthesizing returned wavefields with complex three-dimensional wave equations, this method has a significant advantage in terms of lower computational cost compared to traditional methods. The backhaul wavefield synthesis processing method of this application does not require solving the wave equation based on the backhaul data to obtain the backhaul wavefield. Instead, it determines the backhaul wavelength by using sparse linear filter coefficients and forward wavelength. This significantly reduces computational power, improves processing efficiency, reduces processing time, and reduces memory consumption, thus lowering computational costs. Moreover, the computational complexity is only related to the number of non-zero coefficients of the sparse filter being solved, and is independent of model dimension, boundary conditions, grid order, and the complexity of the wave equation. Therefore, it has greater universality and can be adapted to more application scenarios. Attached Figure Description

[0029] Figure 1 This is a flowchart of the backhaul wavefield synthesis processing method in one embodiment;

[0030] Figure 2a This is an example of a method for processing back-transmitted wavefield synthesis in one embodiment, showing the amplitude map of back-transmitted seismic data.

[0031] Figure 2b This is an amplitude diagram of the sparse filter coefficients in the backhaul wavefield synthesis processing method in one embodiment;

[0032] Figure 2c This is a comparison diagram of the amplitudes of the comparative seismic data and the source wavelet superimposed using a filter in one embodiment of the back-propagation wavefield synthesis processing method.

[0033] Figure 3a , Figure 3c , Figure 3e This is an example of a wavefield synthesis processing method applied to the test synthesis of shot gather data in a wavefield synthesis processing method.

[0034] Figure 3b , Figure 3d , Figure 3fThis is a composite wavefield map used in traditional return-propagation wavefield synthesis methods for testing shot gather data.

[0035] Figure 4 This is a schematic diagram of the processing device for the backhaul wavefield synthesis processing method in one embodiment;

[0036] Figure 5 This is a schematic diagram of the structure of a computer device in one embodiment. Detailed Implementation

[0037] To facilitate understanding of this application and to make the aforementioned objectives, features, and advantages of this application more apparent, a detailed description of specific embodiments of this application is provided below in conjunction with the accompanying drawings. Numerous specific details are set forth in the following description to provide a thorough understanding of this application, and preferred embodiments are shown in the accompanying drawings. However, this application can be implemented in many different forms and is not limited to the embodiments described herein. Rather, these embodiments are provided to provide a more thorough and complete understanding of the disclosure of this application. This application can be implemented in many other ways than those described herein, and those skilled in the art can make similar modifications without departing from the spirit of this application; therefore, this 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 commonly understood by one of ordinary skill in the art to which this application pertains. The terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit this application. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

[0038] Example 1

[0039] Firstly, this application provides a method for synthesizing and processing return wave fields; please refer to [link to relevant documentation]. Figure 1 The method for synthesizing and processing the return wave field includes the following steps:

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

[0041] In this application, both source return wave data and source wavelet data are instrument-received data. Optionally, the source wavelet refers to the seismic wave generated when an earthquake source is excited, which propagates in a viscoelastic medium. Its waveform changes over time, eventually forming a signal with a definite start time, finite energy, and a certain duration. This is the source wavelet. Return wave data refers to the data of seismic waves returning to the surface after reflection underground. When seismic waves encounter a concave interface with a large interface depth or a situation with a very large velocity gradient, a rotation phenomenon occurs, forming a rotating wave. The time-distance curve of the rotating wave is overlapping and intersecting. This is because the reflected waves from 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 the opposite of the order of their corresponding reflection points. The characteristics of the rotating wave are low apparent velocity, increased amplitude at the rotation point, and it always accompanies and interferes with normal reflected waves. Optionally, source wavelet data and return wave data can be received using a seismograph, a highly sensitive seismic wave receiving device that records seismic wave vibrations and converts them into analyzable data. Source wavelet data can be recorded by a seismograph to determine the seismic wave's onset time, energy, and waveform changes. Return wave data can also be recorded by a seismograph to obtain the signal after the seismic wave reflects underground, along with information such as its arrival time and amplitude. Of course, this is not the only option. It should be noted that for details on how to acquire source return wave data and source wavelet data using relevant instruments, please refer to existing technologies; this application will not elaborate further.

[0042] S220: Determine the sparse linear filter coefficients based on the source return wave data and the source wavelet data;

[0043] In this embodiment, since the earthquake backhaul data is equal to the convolution of the source wavelet data and the sparse linear filter, the coefficients of the sparse linear filter are determined based on the source backhaul 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 represents the coefficients of the sparse linear filter.

[0047] To better determine the sparse linear filter coefficients, in one embodiment, the step of determining the sparse linear filter coefficients based on the source return wave data and the source wavelet data includes:

[0048] Based on the source return wave data and the source wavelet data, the coefficients of the sparse linear filter are determined using least squares matching.

[0049] It should be noted that the traditional process of synthesizing the return wave requires solving the wave equation using seismic return data to obtain the seismic return wave, and then solving the wave equation using source wavelet data to obtain the source wavelet. The process of solving the wave equation using seismic return data to obtain the seismic return wave is particularly important, as geophysical methods based on the wave equation, such as reverse time migration and full waveform inversion, often require calculating both the forward propagation wavefield of the source and the return wavefield of the seismic data. The conventional calculation of the return wavefield uses the return seismic data as the source and solves the wave equation using differential equation solving methods such as the finite difference method and the finite element method within a discrete grid of the subsurface physical parameter model. In practical applications of geophysical methods based on pre-stack shot gather records, this conventional method of calculating the return wavefield has two drawbacks: 1. It requires solving the wavefield of all model grid points within the computational range simultaneously, thus requiring a large amount of memory; 2. The computational complexity of the return wavefield increases dramatically with the increase in model dimension and the complexity of the wave equation, thus hindering the practical application of more advanced geophysical methods. In this embodiment, to reduce the computational burden of solving complex wave equations from returned wave data, a strategy of avoiding solving the wave equations is adopted. Instead, the sparse linear filter coefficients are obtained by combining the seismic returned wave data with the source wavelet data (which is equivalent to the convolution of the source wavelet data and the sparse linear filter). Then, the returned wavefield is obtained by convolving the forward propagation wavefield of the source wavelet field with the sparse linear filter coefficients. This method is based on the principle of energy superposition, using the time delay superposition of the forward propagation wavefield to synthesize the returned wavefield from returned seismic data near the source. This method does not require solving the wave equation during returned wavefield synthesis and can synthesize returned wavefields at any time and location within the model. Its computational complexity does not increase with the increase of model dimension or the complexity of the wave equation. The returned wavefield synthesized by this invention does not rely on solving the wave equation during returned wave data and can directly synthesize returned wavefield values ​​at any spatial location and time. The least-squares accuracy of the returned wavefield calculated by this method is high. The computational complexity of this method depends only on the number of non-zero coefficients of the sparse filter being solved, and is independent of the model dimension, boundary conditions, grid order, and the complexity of the wave equation. Therefore, when synthesizing the backpropagation wavefield of a three-dimensional high-order complex wave equation, this method has a significant advantage in terms of low computational cost compared to traditional methods.

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

[0051]

[0052] Where φ(w) is the least squares matching error, w(i) is the coefficient of the sparse linear filter, 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 a set number of iterations. When solving linear equation systems using least squares matching, there exists a convergence accuracy, which is a key factor in determining the number of iterations. The maximum number of iterations required when solving matrix factorization using alternating least squares matching depends on the dimension and coefficient density of the scoring matrix. Specifically, if the original matrix has a low dimension or low coefficient density (sparse), fewer iterations are needed; conversely, if the original matrix has a high dimension or high coefficient density (dense), more iterations are required.

[0054] In this application, the least squares matching error and the number of iterations can be defined to obtain a more suitable sparse linear filter coefficient.

[0055] In one embodiment, the least squares matching formula uses the effective set method to solve for the sparse linear filter coefficients.

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

[0057] (1) Input source return wave data d, source wavelet data f, iteration number 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 Fw

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

[0061] (3) When the number of iterations has not reached the maximum number of iterations, and the maximum value of vector e is greater than the least squares matching error φ(w), the following definition applies:

[0062]

[0063] Move the element index 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}, and ns is the number of seismic data sampling points;

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

[0067] (4) When s P When less than 0, the definition is:

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

[0069] w = w + α(sw)

[0070] When w is less than or equal to 0, the elements in w with values ​​less than 0 are removed from P and added 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 s R Each corresponds to a specific sampling point data.

[0074] definition:

[0075] w = s

[0076] e = F T dF T Fw

[0077] When the maximum value of vector e is less than or equal to the least squares matching error φ(w), the iteration terminates; otherwise, the loop continues to return to step (3).

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

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

[0080]

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

[0082] w=0

[0083] s=0

[0084] e = F T dF T Fw

[0085] iter=1

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

[0087]

[0088] Move the element index 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] Remove the indices of elements in w that are less than 0 from 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 Fw

[0098] iter = iter + 1

[0099] return w

[0100] In this algorithm, F is the Tobleitz matrix of the source wavelet data, n is the maximum number of algorithm iterations, and w is defined as the upper limit of the number of zero coefficients. err defines the loop termination matching error threshold. P s R Let be a subvector of vector s, and let be the partial elements of vector s corresponding to the element indices of P and R.

[0101] S230: By using the numerical values ​​and delays of non-zero sparse linear filter coefficients, the forward propagation wavefield is weighted and superimposed to synthesize the return propagation wavefield.

[0102] In this embodiment, the propagating wavefield is the source wavelet field calculated from the source wavelet data. In other words, the propagating wavefield can also be understood as the source wavelet field. Specifically, the propagating wavefield can be determined by solving the wave equation using the source wavelet data. The process for determining the propagating wavefield can refer to existing technologies and is not limited to those mentioned in this application. Specifically, the propagating wavefield delay of the synthesized return wavefield is determined based on the delay of the non-zero coefficients of the sparse linear filter, and the superposition coefficient is determined based on the values ​​of the filter coefficients.

[0103] In this application, the propagating wavefield is synthesized by weighted superposition of the forward propagation wavefield using the numerical values ​​and delays of non-zero sparse linear filter coefficients. This reduces the computational burden of solving complex wave equations from the propagating wave data, avoiding the strategy of solving the wave equation by combining the seismic propagation data with the source wavelet data (which is equivalent to the convolution of the sparse linear filter) to obtain the sparse linear filter coefficients, and then convolving the forward propagation wavefield of the source wavelet with the sparse linear filter coefficients. This method is based on the principle of energy superposition, using the time delay superposition of the propagating wavefield to synthesize the propagating wavefield from the seismic data near the source. This method does not require solving the wave equation during the synthesis of the propagating wavefield, and can synthesize the propagating wavefield at any time and location in the model. Its computational complexity does not increase with the increase of model dimension or the complexity of the wave equation. The propagating wavefield synthesized by this invention does not depend on solving the wave equation during propagation and can directly synthesize the propagating wavefield value at any spatial location and time. The propagating wavefield calculated by this method has high least-squares accuracy. The computational complexity of this method depends only on the number of non-zero coefficients of the sparse filter being solved, and is independent of the model dimension, boundary conditions, grid order, and the complexity of the wave equation. Therefore, when synthesizing the backpropagation wavefield of a three-dimensional high-order complex wave equation, this method has a significant advantage in terms of low computational cost compared to traditional methods.

[0104] In one embodiment, the step of synthesizing the return wavefield by weighted superposition of the forward propagation wavefield using the values ​​of non-zero sparse linear filter coefficients and the delay is determined by the following formula:

[0105]

[0106] in, x is the location of the earthquake source. s Excitation, at detector position x r The received source record transmits the wavefield value at the model's spatial location x at time t; U f (x, x′) s =x r ,t-τ(i)) is the earthquake source at x r The wave field value at spatial location x after a time delay t-τ(i) for the forward propagation wave field excited at point x. In this way, the return propagation wave field can be determined relatively well.

[0107] The aforementioned method for synthesizing the returned wavefield is based on the assumption that the source-returned data is the convolution of source wavelet data and a sparse linear filter. The coefficients of the sparse linear filter are determined using the source-returned data and the source wavelet data. Then, according to the superposition principle, the forward-propagating wavefield is weighted and superimposed based on the delay and value of the non-zero coefficients of the sparse linear filter to synthesize the returned wavefield. In contrast, traditional methods for synthesizing returned wavefields require solving the wave equation using the returned data. This method, unlike traditional methods for synthesizing returned wavelengths, does not rely on solving the wave equation during the return process and can directly synthesize returned wavefield values ​​at any spatial location and time. The computational complexity of this method depends only on the number of non-zero coefficients of the solved sparse filter and is independent of the model dimension, boundary conditions, grid order, and the complexity of the wave equation. Therefore, when synthesizing returned wavefields with complex three-dimensional wave equations, this method has a significant advantage in terms of lower computational cost compared to traditional methods. The backhaul wavefield synthesis processing method of this application does not require solving the wave equation based on the backhaul data to obtain the backhaul wavefield. Instead, it determines the backhaul wavelength by using sparse linear filter coefficients and forward wavelength. This significantly reduces computational power, improves processing efficiency, reduces processing time, and reduces memory consumption, thus lowering computational costs. Moreover, the computational complexity is only related to the number of non-zero coefficients of the sparse filter being solved, and is independent of model dimension, boundary conditions, grid order, and the complexity of the wave equation. Therefore, it has greater universality and can be adapted to more application scenarios.

[0108] Example 2

[0109] This embodiment illustrates the specific application process of the proposed method for synthesizing and processing the returned wavefield. Based on a two-dimensional complex underground physical model, shot gather data is used to test the method of this invention. The accuracy and computational cost of the returned wavefield synthesized by this method are compared with those calculated by the traditional finite difference method (using second-order time and fourth-order spatial grids) for solving the wave equation. Please refer to [link to relevant documentation]. Figure 2a , Figure 2b and Figure 2c , Figure 2a To transmit earthquake data back, Figure 2b These are the sparse filter coefficients calculated using the effective set method in the backhaul wavefield synthesis processing method of this application. Figure 2c For a comparison of this application's seismic data with source wavelets obtained by time-delay superposition using sparse filter coefficients, please refer to [link to relevant documentation]. Figures 3a to 3f The left side Figure 3a , Figure 3c , Figure 3e The wavefield for the back-propagation wavefield synthesis processing method of this application is shown on the left. Figure 3b , Figure 3d , Figure 3f The wave field calculated using traditional finite difference methods (using second-order time and fourth-order spatial grids) is transformed by... Figures 3a to 3f Comparing wavefield snapshots at different propagation times, the results show that the wavefield synthesis processing method provided in this application achieves a 99.5% match with traditional methods, but its computational cost is only 50% of that of traditional methods. Therefore, it can reduce computational costs, and its computational complexity is only related to the number of non-zero coefficients of the solved sparse filter, and is independent of model dimension, boundary conditions, grid order, and the complexity of the wave equation. This makes it more universal and adaptable to more application scenarios. The wavefield synthesis processing method in this application, because it does not rely on solving the wave equation from the backhaul data to obtain the backhaul wavefield, but determines the backhaul wavelength through sparse linear filter coefficients and the forward propagation wavelength, significantly reduces computational power, improves processing efficiency, reduces processing time, and reduces memory consumption. This reduces computational costs, and its computational complexity is only related to the number of non-zero coefficients of the solved sparse filter, and is independent of model dimension, boundary conditions, grid order, and the complexity of the wave equation. Therefore, it is more universal and adaptable to more application scenarios.

[0110] Example 3

[0111] Secondly, this application also provides a backhaul wavefield synthesis processing apparatus, please refer to [link to relevant documentation]. Figure 4 The backhaul wavefield 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 acquire source return wave data and source wavelet data respectively;

[0113] The filter coefficient determination module is used to determine the sparse linear filter coefficients based on the source return wave data and the source wavelet data acquired by the acquisition module.

[0114] The synthesis module uses the values ​​and delays of the non-zero sparse linear filter coefficients given by the filter coefficient determination module to weight and superimpose the forward propagation wave field to synthesize the return propagation wave field.

[0115] The aforementioned echo wavefield synthesis processing device employs the acquisition module, filter coefficient determination module, and synthesis module. Based on the fact that the source echo data is a convolution of source wavelet data and a sparse linear filter, the sparse linear filter coefficients are determined using the source echo data and source wavelet data. Then, according to the superposition principle, the forward propagation wavefield is weighted and superimposed based on the delay and value of the non-zero coefficients of the sparse linear filter to synthesize the echo wavefield. In contrast, traditional echo wavefield synthesis methods require solving the wave equation using the echo data to obtain the echo wavefield. This application, compared to traditional echo wavelength synthesis methods, does not rely on solving the wave equation during echo propagation and can directly synthesize echo wavefield values ​​at any spatial location and any time. The computational complexity of this method is only related to the number of non-zero coefficients of the solved sparse filter and is independent of model dimension, boundary conditions, grid order, and the complexity of the wave equation. Therefore, when synthesizing echo wavefields with three-dimensional high-order complex wave equations, this method has a significant advantage in terms of low computational cost compared to traditional methods. The backhaul wavefield synthesis processing device of this application does not require solving the wave equation based on the backhaul data to obtain the backhaul wavefield. Instead, it determines the backhaul wavelength by using sparse linear filter coefficients and forward wavelength. This significantly reduces computational power, improves processing efficiency, reduces processing time, and reduces memory consumption, thus lowering computational costs. Moreover, the computational complexity is only related to the number of non-zero coefficients of the sparse filter being solved, and is independent of model dimension, boundary conditions, grid order, and the complexity of the wave equation. Therefore, it has greater versatility and can be adapted 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 embodiment, the filter coefficient determination module is used to determine sparse linear filter coefficients based on least squares matching, according to the source return wave data and the source wavelet data.

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

[0119] d = wf

[0120] In the above formula, w represents the coefficients of the sparse linear filter.

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

[0122]

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

[0124] Optionally, least squares matching has a set number of iterations. When solving linear equation systems using least squares matching, there exists a convergence accuracy, which is a key factor in determining the number of iterations. The maximum number of iterations required when solving matrix factorization using alternating least squares matching depends on the dimension and coefficient density of the scoring matrix. Specifically, if the original matrix has a low dimension or low coefficient density (sparse), fewer iterations are needed; conversely, if the original matrix has a high dimension or high coefficient density (dense), more iterations are required.

[0125] In one embodiment, 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 least squares matching, according to the source return wave data and the source wavelet data:

[0127]

[0128] Where φ(w) is the least squares matching error, w(i) is the coefficient of the sparse linear filter, 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 to solve the sparse linear filter coefficients for the formula using the effective set method.

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

[0131] (1) Input source return wave data d, source wavelet data f, iteration number 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 Fw

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

[0135] (3) When the number of iterations has not reached the maximum number of iterations, and the maximum value of vector e is greater than the least squares matching error φ(w), the following definition applies:

[0136]

[0137] Move the element index 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}, and ns is the number of seismic data sampling points;

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

[0141] (4) When s P When less than 0, the definition is:

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

[0143] w = w + α(sw)

[0144] When w is less than or equal to 0, the elements in w with values ​​less than 0 are removed from P and added 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 s R Each corresponds to a specific sampling point data.

[0148] definition:

[0149] w = s

[0150] e = F T dF T Fw

[0151] When the maximum value of vector e is less than or equal to the least squares matching error φ(w), the iteration terminates; otherwise, the loop continues to return to step (3).

[0152] In one embodiment, the synthesis module is used to synthesize the return wavefield by weighted superposition of the forward propagation wavefield using the values ​​of non-zero sparse linear filter coefficients and the delay, wherein the return wavefield is determined using the following formula:

[0153]

[0154] in, x is the location of the earthquake source. s Excitation, at detector position x r The received source record transmits the wavefield value at the model's spatial location x at time t; U f (x, x′) s =x r ,t-τ(i)) is the earthquake source at x r The wave field value at spatial location x after a time delay t-τ(i) for the forward propagation wave field excited at point x. In this way, the return propagation wave field can be determined relatively well.

[0155] The aforementioned echo wavefield synthesis processing device is based on the fact that the source echo data is the convolution of source wavelet data and a sparse linear filter. It determines the coefficients of the sparse linear filter using the source echo data and source wavelet data, and then, according to the superposition principle, weighted superposition of the forward propagation wavefield based on the delay and value of the non-zero coefficients of the sparse linear filter to synthesize the echo wavefield. In contrast, traditional echo wavefield synthesis methods require solving the wave equation using the echo data to obtain the echo wavefield. This application, compared to traditional echo wavelength synthesis methods, does not rely on solving the wave equation during echo propagation and can directly synthesize echo wavefield values ​​at any spatial location and any time. The computational complexity of the echo wavefield synthesis processing device provided in this application is only related to the number of non-zero coefficients of the solved sparse filter, and is independent of the model dimension, boundary conditions, grid order, and the complexity of the wave equation. Therefore, when synthesizing echo wavefields with three-dimensional high-order complex wave equations, this method has a significant advantage in terms of low computational cost compared to traditional methods. The backhaul wavefield synthesis processing device of this application does not require solving the wave equation based on the backhaul data to obtain the backhaul wavefield. Instead, it determines the backhaul wavelength by using sparse linear filter coefficients and forward wavelength. This significantly reduces computational power, improves processing efficiency, reduces processing time, and reduces memory consumption, thus lowering computational costs. Moreover, the computational complexity is only related to the number of non-zero coefficients of the sparse filter being solved, and is independent of model dimension, boundary conditions, grid order, and the complexity of the wave equation. Therefore, it has greater versatility and can be adapted to more application scenarios.

[0156] Example 4

[0157] Thirdly, this application also provides a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the steps of the backhaul wave field synthesis processing method as described in any of the above embodiments.

[0158] In one embodiment, the computer device can be as follows: Figure 5As shown. The computer device includes a processor, memory, network interface, display screen, and input devices connected via a system bus. The processor provides computing and control capabilities. The memory includes a non-volatile storage medium and internal memory. The non-volatile storage medium stores the operating system and computer programs. The internal memory provides an environment for the operation of the operating system and computer programs in the non-volatile storage medium. The network interface is used to communicate with a server via a network connection. When the computer program is executed by the processor, it implements a return wave field synthesis processing method as described in any of the above embodiments. The display screen can be a liquid crystal display (LCD) or an e-ink display. The input devices can be a touch layer covering the display screen, buttons, a trackball, or a touchpad mounted on the computer device casing, or an external keyboard, touchpad, or mouse.

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

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

[0161] Separately acquire source return wave data and source wavelet data;

[0162] The coefficients of the sparse linear filter are determined based on the source return wave data and the source wavelet data.

[0163] By using the numerical values ​​of non-zero sparse linear filter coefficients and their delays, the forward propagation wavefield is weighted and superimposed to synthesize the return propagation wavefield.

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

[0165] The step of determining the sparse linear filter coefficients based on the source return wave data and the source wavelet data includes:

[0166] Based on the source return wave data and the source wavelet data, the coefficients of the sparse linear filter are determined using least squares matching.

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

[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] Where φ(w) is the least squares matching error, w(i) is the filter coefficient, f(τ(i)) is the source wavelet data after time delay τ(i), and d is the source return wave data.

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

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

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

[0174] In the step of synthesizing the return wavefield by weighted superposition of the forward propagation wavefield using the values ​​of non-zero sparse linear filter coefficients and the delay, the return wavefield is determined using the following formula:

[0175]

[0176] in, x is the location of the earthquake source. s Excitation, at detector position x r The received source record transmits the wavefield value at the model's spatial location x at time t; U f (x, x′) s =x r ,t-τ(i)) is the earthquake source at x r The wave field value at spatial location x after a time delay t-τ(i) for the propagating wave field excited at point x.

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

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

[0179] (1) Input source return wave data d, source wavelet data f, iteration number 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 Fw

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

[0183] (3) When the number of iterations has not reached the maximum number of iterations, and the maximum value of vector e is greater than the least squares matching error φ(w), the following definition applies:

[0184]

[0185] Move the element index 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}, and ns is the number of seismic data sampling points;

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

[0189] (4) When s P When less than 0, the definition is:

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

[0191] w = w + α(sw)

[0192] When w is less than or equal to 0, the elements in w with values ​​less than 0 are removed from P and added 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 s R Each corresponds to a specific sampling point data.

[0196] definition:

[0197] w = s

[0198] e = F T dF T Fw

[0199] When the maximum value of vector e is less than or equal to the least squares matching error φ(w), the iteration terminates; otherwise, the loop continues to return to step (3).

[0200] Example 5

[0201] Fourthly, this application also provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the backhaul wave field synthesis processing method as described in any of the above embodiments.

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

[0203] Separately acquire source return wave data and source wavelet data;

[0204] The coefficients of the sparse linear filter are determined based on the source return wave data and the source wavelet data.

[0205] By using the numerical values ​​of non-zero sparse linear filter coefficients and their delays, the forward propagation wavefield is weighted and superimposed to synthesize the return propagation wavefield.

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

[0207] The step of determining the sparse linear filter coefficients based on the source return wave data and the source wavelet data includes:

[0208] Based on the source return wave data and the source wavelet data, the coefficients of the sparse linear filter are determined using least squares matching.

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

[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] Where φ(w) is the least squares matching error, w(i) is the filter coefficient, f(τ(i)) is the source wavelet data after time delay τ(i), and d is the source return wave data.

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

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

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

[0216] In the step of synthesizing the return wavefield by weighted superposition of the forward propagation wavefield using the values ​​of non-zero sparse linear filter coefficients and the delay, the return wavefield is determined using the following formula:

[0217]

[0218] in, x is the location of the earthquake source. s Excitation, at detector position x r The received source record transmits the wavefield value at the model's spatial location x at time t; U f (x, x′) s =x r ,t-τ(i)) is the earthquake source at x r The wave field value at spatial location x after a time delay t-τ(i) for the propagating wave field excited at point x.

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

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

[0221] (1) Input source return wave data d, source wavelet data f, iteration number 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 Fw

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

[0225] (3) When the number of iterations has not reached the maximum number of iterations, and the maximum value of vector e is greater than the least squares matching error φ(w), the following definition applies:

[0226]

[0227] Move the element index 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}, and ns is the number of seismic data sampling points;

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

[0231] (4) When s P When less than 0, the definition is:

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

[0233] w = w + α(sw)

[0234] When w is less than or equal to 0, the elements in w with values ​​less than 0 are removed from P and added 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 s R Each corresponds to a specific sampling point data.

[0238] definition:

[0239] w = s

[0240] e = F T dF T Fw

[0241] When the maximum value of vector e is less than or equal to the least squares matching error φ(w), the iteration terminates; otherwise, the loop continues to return to step (3).

[0242] Those skilled in the art will understand that all or part of the processes in the methods of the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium, and when executed, it can include the processes of the embodiments of the above methods. Any references to memory, storage, databases, 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. By way of illustration and not limitation, RAM is available in various forms, such as static RAM (SRAM), dynamic RAM (DRAM), synchronous DRAM (SDRAM), dual data rate SDRAM (DDRSDRAM), enhanced SDRAM (ESDRAM), synchronous link DRAM (SLDRAM), Rambus direct RAM (RDRAM), direct memory bus dynamic RAM (DRDRAM), and memory bus dynamic RAM (RDRAM), etc.

[0243] The aforementioned method and apparatus for synthesizing and processing return wavefields are based on the premise that the source return data is the convolution of source wavelet data and a sparse linear filter. The coefficients of the sparse linear filter are determined using the source return data and source wavelet data. Then, according to the superposition principle, the time delays of the forward propagation wavefield are weighted and superimposed based on the delays and values ​​of the non-zero coefficients of the sparse linear filter to synthesize the return wavefield. In contrast, traditional return wavefield synthesis methods require solving the wave equation using the return data to obtain the return wavefield. This application, compared to traditional methods for synthesizing return wavelengths, does not rely on solving the wave equation during return and can directly synthesize return wavefield values ​​at any spatial location and any time. The computational complexity of this method is only related to the number of non-zero coefficients of the solved sparse filter and is independent of the model dimension, boundary conditions, grid order, and the complexity of the wave equation. Therefore, when synthesizing return wavefields with three-dimensional high-order complex wave equations, this method has a significant advantage in terms of low computational cost compared to traditional methods. The backhaul wavefield synthesis processing method of this application does not require solving the wave equation based on the backhaul data to obtain the backhaul wavefield. Instead, it determines the backhaul wavelength by using sparse linear filter coefficients and forward wavelength. This significantly reduces computational power, improves processing efficiency, reduces processing time, and reduces memory consumption, thus lowering computational costs. Moreover, the computational complexity is only related to the number of non-zero coefficients of the sparse filter being solved, and is independent of model dimension, boundary conditions, grid order, and the complexity of the wave equation. Therefore, it has greater universality and can be adapted to more application scenarios.

[0244] The technical features of the embodiments described above can be combined arbitrarily. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as the combination of these technical features does not contradict each other, it should be considered within the scope of this specification. It should be noted that the terms "in one embodiment," "for example," and "again," etc., in this application are intended to illustrate the application and not to limit it. The embodiments described above only illustrate several implementation methods of this application, and their descriptions are relatively specific and detailed, but they should not be construed as limiting the scope of the patent application. It should be pointed out that for those skilled in the art, several modifications and improvements can be made without departing from the concept of this application, and these all fall within the protection scope of this application. Therefore, the protection scope of this patent application should be determined by the appended claims.

Claims

1. A method for back-propagating wavefield synthesis processing, characterized in that, The method comprises the following steps: obtaining source return wave data and source wavelet data respectively; determining sparse linear filter coefficients according to the source return wave data and the source wavelet data; performing time delay weighted stacking on a forward wave field using the values and delays of the non-zero sparse linear filter coefficients to synthesize a return wave field; The step of determining the sparse linear filter coefficients comprises: determining the sparse linear filter coefficients based on least square matching according to the source return wave data and the source wavelet data; In the step of determining the sparse linear filter coefficients based on least square matching, the sparse linear filter coefficients are determined using the following formula: wherein, is the least square matching error, is the filter coefficient, is the time delay the source wavelet data after, d is the source echo wave data.

2. The process of backpropagating wavefield synthesis according to claim 1, wherein, In the formula of the least square matching, the effective set method is used to solve the sparse linear filter coefficients.

3. The process of backpropagating wavefield synthesis according to claim 1, wherein, In the step of performing time delay weighted stacking on a forward wave field using the values and delays of the non-zero sparse linear filter coefficients to synthesize a return wave field, the return wave field is determined using the following formula: wherein, is the source position excitation, at the receiver position received source record backpropagates the wavefield value at the model space position at time t; is the source excitation at forward propagated wavefield at time delay at spatial position at time t.

4. A device for backpropagating wavefield synthesis processing, characterized by The method comprises the following steps: an obtaining module, configured to obtain source return wave data and source wavelet data respectively; a filter coefficient determining module, configured to determine sparse linear filter coefficients according to the source return wave data and the source wavelet data obtained by the obtaining module; The filter coefficient determining module is configured to determine the sparse linear filter coefficients based on least square matching according to the source return wave data and the source wavelet data. The filter coefficient determining module determines the sparse linear filter coefficients using the following formula: wherein, is the least square matching error, is the filter coefficient, is the time delay the source wavelet data after, d is the source echo wave data; a synthesizing module, configured to perform time delay weighted stacking on a forward wave field using the values and delays of the non-zero sparse linear filter coefficients given by the filter coefficient determining module to synthesize a return wave field.

5. A computer device comprising a memory, a processor, and a computer program stored on the memory and executable on the processor, characterized in that, The processor implements the steps of the return wave field synthesis processing method of any one of claims 1-3 when executing the program.

6. A computer-readable storage medium having stored thereon a computer program, characterized in that, The computer program implements the steps of the return wave field synthesis processing method of any one of claims 1-3 when executed by the processor.

Images