A high-resolution processing method for seismic data based on point spread function

Abstract

The present invention discloses a method for processing seismic data with high resolution based on a point spread function, comprising: inputting seismic records, conventional migration results and migration velocity, acquiring an observation system based on trace header information of the seismic records; calculating an initial point spread function based on source wavelet, travel time and amplitude; performing deconvolution processing on the migration results based on the point spread function to obtain an imaging result after deconvolution; updating the point spread function based on the imaging result and the conventional migration result, and outputting an image domain least squares migration imaging result and an updated point spread function; the present invention converts a strongly nonlinear parameter estimation inverse problem into an image domain least squares inversion imaging problem with better convexity; introducing an asymptotic Green's function under the WKBJ approximation to solve the computational problem of the Hessian matrix, and then calculating the point spread function of any underground point; introducing a priori knowledge of the point spread function and the reflection coefficient, constructing a strategy for iteratively updating the point spread function and the reflection coefficient, and obtaining a higher-resolution imaging result.

Description

Technical Field

[0001] The present invention relates to the field of geophysical exploration technology, and in particular to a high-resolution processing method for seismic data based on a point spread function. Background Art

[0002] High-resolution and high-fidelity seismic inversion imaging can be divided into two main routes. The first route is linear inversion based on the Bayesian estimation theory framework. This is the so-called least squares migration imaging, which is a linear sub-problem of full waveform inversion. The advantage of this method is that it has relatively low requirements on the observation system and the degree of data regularity, and the migration operator in the iterative solution process can flexibly select ray or wave types; the disadvantage of this method is that the computational efficiency is low and there are convergence problems for actual complex data. Another development route for high-fidelity imaging is the true amplitude imaging method under high-frequency approximation.

[0003] Least-squares migration imaging, based on the Bayesian estimation framework, can be implemented in either the data domain or the imaging domain. Generally, data-domain least-squares migration imaging typically requires 10 or more iterations, and its computational cost is more than an order of magnitude higher than conventional migration imaging. Furthermore, the inclusion of prior information is crucial, given factors such as inaccurate initial models, the inability of the Born approximation to accurately model wavefield perturbations, noise that does not conform to a Gaussian distribution, irregular observation systems, and unknown and space-varying wavelets. Without effective data preconditioning or prior information, iterative convergence of data-domain least-squares migration imaging can be slow or even non-convergent. Solving the normal equation for the error functional of data-domain least-squares migration imaging when the gradient is zero is the key to implementing least-squares migration imaging in the imaging domain. The key is to explicitly compute the Hessian matrix. Essentially, a column of the Hessian matrix corresponds to the point spread function at a point in the model space. This function depends on the propagating wavelet, the observation system, and the background velocity. It is a band-limited wave packet with a finite spatial distribution and exhibits locality and space-varying properties. Compared to the iterative least-squares migration method in the data domain, imaging-domain least-squares migration avoids the calculation of migration, demigration, and step size, resulting in a smaller overall computational load and more flexible operation. Therefore, imaging-domain least-squares migration is a more reasonable choice for practical applications.

[0004] However, due to the unknown source wavelet, the accuracy of background velocity and the precision of forward modeling / migration operators, the calculation results of point spread cannot be accurate. Therefore, it is necessary to construct a least squares migration imaging method in the imaging domain to correct the inaccurate point spread function. Summary of the Invention

[0005] The purpose of this section is to summarize some aspects of the embodiments of the present invention and briefly introduce some preferred embodiments. Some simplifications or omissions may be made in this section and the abstract and title of this application to avoid obscuring the purpose of this section, the abstract and the title of the invention, and such simplifications or omissions should not be used to limit the scope of the present invention.

[0006] In view of the above-mentioned problems, the present invention is proposed.

[0007] Therefore, the technical problem solved by the present invention is: the imaging artifact problem caused by the inaccuracy of the point spread function in the least squares migration imaging of seismic data imaging domain.

[0008] To solve the above technical problems, the present invention provides the following technical solution: a method for high-resolution processing of seismic data based on a point spread function, comprising:

[0009] Inputting seismic records, conventional migration results and migration velocity, and acquiring an observation system based on trace header information of the seismic records;

[0010] Analytically calculate the initial point spread function based on the travel time and amplitude of the source wavelet and spatial position;

[0011] performing deconvolution processing on the migration result based on the point spread function to obtain a deconvolution imaging result;

[0012] Based on the deconvolution imaging result and the conventional migration result, updating the point spread function, and finally outputting the image domain least squares migration imaging result and the updated point spread function;

[0013] As a preferred solution of the method for high-resolution processing of seismic data based on point spread function described in the present invention, wherein: the input seismic record includes

[0014] Input earthquake records in SEGY and SU formats.

[0015] As a preferred solution of the method for high-resolution processing of seismic data based on point spread function described in the present invention, the acquisition observation system includes:

[0016] Obtain the parameters of the observation system, namely the spatial position of the shot point and the spatial position of the receiver point.

[0017] As a preferred solution of the method for high-resolution processing of seismic data based on point spread function described in the present invention, the source wavelet includes:

[0018] It is assumed that the general Ricker wavelet is used as the source wavelet parameter and substituted into the analytical formula of the point spread function to calculate the initial point spread function.

[0019] As a preferred solution of the method for high-resolution processing of seismic data based on point spread function described in the present invention, the analytical formula of the point spread function is expressed as:

[0020]

[0021] in, is a spatial point x i The point spread function at R ss represents the autocorrelation of the source wavelet s(ω); T represents the travel time, A represents the amplitude value, and the amplitude of the scattered ray path is given by connecting the two incident rays: A(x r ;x j ;x s )=A(x j ,x r )A(x j ,x s ).

[0022] As a preferred solution of the method for high-resolution processing of seismic data based on point spread function described in the present invention, the analytical formula of the point spread function includes:

[0023] To solve the computational complexity problem of the Hessian matrix, the asymptotic Green's function under the WKBJ approximation is introduced to obtain an approximate expression of the Hessian matrix, namely the analytical formula of the point spread function, which is expressed as the integral of the ray amplitude and the autocorrelation of the source wavelet. In this way, only the travel time, amplitude and source wavelet are needed to calculate the point spread function of any point underground.

[0024] As a preferred embodiment of the method for high-resolution processing of seismic data based on point spread function described in the present invention, the deconvolution process includes:

[0025] By considering the sparsity and structural continuity of the reflection coefficient, L1 constraint and total variation constraint are introduced to construct the error functional J corresponding to the deconvolution process. r (r), and then realize deconvolution processing;

[0026] The error functional J corresponding to the deconvolution process r (r) is expressed as:

[0027]

[0028] Among them, * represents the convolution operator, F (i) is the point spread function of the i-th iteration, r represents the deconvolution in the current iteration, λ1, λ2 and λ3 represent hyperparameters, ‖r‖ p , ‖r‖ TVand Δr represent the Lp constraint, total variation constraint and second-order derivative constraint introduced on the deconvolution result, respectively, p∈[0,2].

[0029] As a preferred embodiment of the method for high-resolution processing of seismic data based on point spread function described in the present invention, the conventional migration imaging result includes:

[0030] The conventional migration imaging result is regarded as the convolution of the point spread function and the reflection coefficient. For the one-dimensional case, the convolution operation is written as the Toeplitz matrix product form. For the two-dimensional or three-dimensional case, the convolution operation and the correlation operation are regarded as a pair of conjugate operators to perform gradient updates. Among them, the correlation operation specifically calculates the correlation coefficient between two variables.

[0031] As a preferred solution of the method for high-resolution processing of seismic data based on point spread function described in the present invention, the updating of point spread function includes:

[0032] Introducing the smooth constraint μ1‖F‖ of the point spread function q , construct the error functional J F (F) to update the point spread function;

[0033] The error functional J at this time F (F) is expressed as:

[0034]

[0035] in, represents the correlation operator, r (i) is the deconvolution result of the i-th iteration, ‖F‖ q is the Lq constraint introduced for the point spread function, μ1 is its corresponding hyperparameter, q∈(1,2].

[0036] As a preferred embodiment of the method for high-resolution processing of seismic data based on point spread function described in the present invention, the output image domain least squares migration imaging result and the updated point spread function include:

[0037] Determine whether the error between the convolution of the point spread function and the deconvolution result and the conventional migration imaging result is less than a given threshold; when the error between the convolution of the point spread function and the deconvolution result and the conventional migration imaging result is less than the given threshold, output the image domain least squares migration imaging result and the updated point spread function.

[0038] The beneficial effects of the present invention are as follows: starting from probabilistic statistical inversion, the present invention transforms the strongly nonlinear parameter estimation inverse problem into a least squares inversion imaging problem with better convexity in the image domain; for the problem of the computational complexity of the Hessian matrix, the asymptotic Green's function under the WKBJ approximation is introduced to obtain an approximate expression of the Hessian matrix, thereby realizing the calculation of the point spread function of any point underground; in addition, the present invention introduces a priori knowledge of the point spread function and the reflection coefficient, constructs a strategy for iteratively updating the point spread function and the reflection coefficient, and obtains higher-resolution imaging results. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following briefly introduces the drawings required for describing the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. Those skilled in the art can also derive other drawings based on these drawings without inventive effort. Among them:

[0040] Figure 1 An overall flow chart of a method for high-resolution processing of seismic data based on a point spread function provided by one embodiment of the present invention;

[0041] Figure 2 A reflection coefficient profile provided for one embodiment of the present invention;

[0042] Figure 3 A conventional migration imaging result provided by one embodiment of the present invention;

[0043] Figure 4 The deconvolution result of the initial point spread function provided by one embodiment of the present invention;

[0044] Figure 5 A high-resolution imaging result provided by one embodiment of the present invention;

[0045] Figure 6 A comparison diagram of the initial point spread function, the actual point spread function, and the last updated point spread function measurement channel provided by one embodiment of the present invention;

[0046] Figure 7 A two-dimensional schematic diagram of a real point spread function provided by one embodiment of the present invention;

[0047] Figure 8 A two-dimensional schematic diagram of the last updated point spread function provided by one embodiment of the present invention;

[0048] Figure 9 A two-dimensional schematic diagram of an initial point spread function provided by one embodiment of the present invention. DETAILED DESCRIPTION

[0049] To make the above-mentioned objects, features, and advantages of the present invention more clearly understood, the following detailed description of the specific embodiments of the present invention is given in conjunction with the accompanying drawings. It is obvious that the described embodiments are only part of the embodiments of the present invention, not all of them. Based on the embodiments of the present invention, all other embodiments obtained by ordinary persons in this field without creative work should fall within the scope of protection of the present invention.

[0050] In the following description, many specific details are set forth to facilitate a full understanding of the present invention. However, the present invention may also be implemented in other ways different from those described herein. Those skilled in the art may make similar generalizations without violating the connotation of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0051] Secondly, the term "one embodiment" or "embodiment" herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The phrase "in one embodiment" appearing in various places throughout this specification does not necessarily refer to the same embodiment, nor does it refer to a separate or selective embodiment that is mutually exclusive of other embodiments.

[0052] The present invention is described in detail with reference to schematic diagrams. For ease of illustration, cross-sectional views of device structures may be partially enlarged and not to scale when describing embodiments of the present invention. Furthermore, the schematic diagrams are merely illustrative and should not limit the scope of the present invention. Furthermore, in actual production, the three-dimensional dimensions of length, width, and depth should be included.

[0053] In the description of the present invention, it should be noted that the terms "upper, lower, inner, and outer" and other references to orientations or positional relationships are based on the orientations or positional relationships shown in the accompanying drawings and are intended solely to facilitate and simplify the description of the present invention. They are not intended to indicate or imply that the devices or components referred to must have, be constructed, or operate in a specific orientation, and therefore should not be construed as limitations on the present invention. Furthermore, the terms "first, second, or third" are used for descriptive purposes only and should not be construed as indicating or implying relative importance.

[0054] In this disclosure, unless otherwise specified or limited, the terms "mounted," "connected," and "connected" should be interpreted broadly. For example, they may refer to fixed, removable, or integral connections. They may also refer to mechanical, electrical, or direct connections, indirect connections through an intermediary, or internal communication between two components. Those skilled in the art will understand the specific meanings of these terms in this disclosure.

[0055] Example 1

[0056] Reference Figure 1, as one embodiment of the present invention, provides a method for high-resolution processing of seismic data based on a point spread function, comprising:

[0057] S1: Input seismic records, conventional migration results and migration velocity, and obtain an observation system based on the trace header information of the seismic records;

[0058] Furthermore, input earthquake records in SEGY format and SU format;

[0059] It should be noted that for seismic data in SEGY and SU formats, each trace has a trace header that records the relevant information of the detector and the source.

[0060] Furthermore, the observation system parameters are obtained by scanning the trace header information of the seismic data;

[0061] It should be noted that the observation system parameters are specifically the spatial position of the shot point and the spatial position of the receiver point.

[0062] S2: Analytical calculation of the initial point spread function based on the travel time and amplitude of the source wavelet and spatial position;

[0063] Furthermore, it is assumed that the Ricker wavelet with general significance is substituted into the analytical formula of the point spread function as the source wavelet parameter to calculate the initial point spread function;

[0064] It should be noted that if the earthquake source in the work area is known, the known earthquake source wavelet can also be input.

[0065] Furthermore, to solve the computational problem of the Hessian matrix, the asymptotic Green's function under the WKBJ approximation is introduced, and the approximate expression of the Hessian matrix, namely the analytical formula of the point spread function, is obtained, which is expressed as the integral of the ray amplitude and the source wavelet autocorrelation;

[0066] The analytical formula of the point spread function is expressed as:

[0067]

[0068] in, is a spatial point x i The point spread function at R ss represents the autocorrelation of the source wavelet s(ω); T represents the travel time, A represents the amplitude value, and the amplitude of the scattered ray path is given by connecting the two incident rays: A(x r ;x j ;x s )=A(x j ,x r )A(x j ,x s ).

[0069] It should be noted that after introducing the asymptotic Green's function under the WKBJ approximation and obtaining the approximate expression of the Hessian matrix, it is only necessary to calculate the propagation time and amplitude of the ray and the source wavelet to calculate the point spread function of any point underground; the propagation time and amplitude of the ray can be obtained by ray tracing.

[0070] S3: performing deconvolution processing on the migration result based on the point spread function to obtain a deconvolution imaging result;

[0071] Furthermore, by considering the sparsity and structural continuity of the reflection coefficient, L1 constraint and total variation constraint are introduced to construct the error functional J corresponding to the deconvolution process. r (r), and then realize deconvolution processing;

[0072] The error functional J corresponding to the deconvolution process r (r) is expressed as:

[0073]

[0074] Among them, * represents the convolution operator, F (i) is the point spread function of the i-th iteration, r represents the deconvolution in the current iteration, λ1, λ2 and λ3 represent hyperparameters, ‖r‖ p , ‖r‖ TV and Δr represent the Lp constraint, total variation constraint and second-order derivative constraint introduced on the deconvolution result, respectively, p∈[0,2].

[0075] It should be known that deconvolution processing is to remove the blur from a blurred image to obtain a clearer image. The imaging result after deconvolution refers to the result after deconvolution processing.

[0076] It should be noted that the sparsity of the reflection coefficient corresponds to the long-tail distribution in probability statistics; the continuity and smoothness of the geological structure correspond to total variation regularization and second-order gradient regularization, respectively; based on the above two prior cognitions and considering the maximization of the posterior probability density, the L1 constraint and total variation constraint are introduced to construct the error functional J r (r), thereby realizing deconvolution processing.

[0077] Furthermore, it is determined whether the error between the convolution of the point spread function and the deconvolution result and the conventional migration imaging result is less than a given threshold; if it is less than the given threshold, the convergence state is reached, and the image domain least squares migration imaging result and the updated point spread function are output; otherwise, the point spread function update step is entered.

[0078] S4: Based on the deconvolution imaging result and the conventional migration result, the point spread function is updated, and finally the image domain least squares migration imaging result and the updated point spread function are output.

[0079] Furthermore, the conventional migration imaging result is regarded as the convolution result of the point spread function and the reflection coefficient. For the one-dimensional case, the convolution operation is written as the Toeplitz matrix product form. For the two-dimensional or three-dimensional case, the convolution operation and the correlation operation are regarded as a pair of conjugate operators to perform gradient update.

[0080] It should be noted that the correlation operation is to calculate the correlation coefficient between two variables to measure the degree of similarity.

[0081] Furthermore, we introduce the smooth constraint μ1‖F‖ of the point spread function q , construct the error functional J F (F) to update the point spread function;

[0082] Error functional J F (F) is expressed as:

[0083]

[0084] in, represents the correlation operator, r (i) is the deconvolution result of the i-th iteration, ‖F‖ q is the Lq constraint introduced for the point spread function, μ1 is its corresponding hyperparameter, q∈(1,2].

[0085] It should be noted that the smoothness of the point spread function corresponds to the Gaussian distribution in probability statistics. Based on the prior knowledge of the point spread function and the maximization of the posterior probability density, the smoothness constraint μ1‖F‖ of the point spread function is introduced. q , to construct the error functional J F (F), thereby updating the point spread function.

[0086] Furthermore, it is determined whether the error between the convolution of the point spread function and the deconvolution result and the conventional migration imaging result is less than a given threshold; if it is less than the given threshold, the convergence state is reached, and the image domain least squares migration imaging result and the updated point spread function are output; otherwise, the deconvolution processing step will be performed until convergence is completed.

[0087] It should be noted that the algorithm can converge when the number of iterations is less than 20.

[0088] Example 2

[0089] Reference Figures 2 to 9, which is an embodiment of the present invention, provides a high-resolution processing method for seismic data based on point spread function. In order to verify the beneficial effects of the present invention, scientific demonstration is carried out through economic benefit calculation and simulation experiments.

[0090] First, considering the structure of caves, thin interbeds and inclined strata, a velocity model is established and the corresponding reflection coefficient profile is calculated, such as Figure 2 As shown, the model length and width are 101 by 101 sampling points, and the sampling interval is 5 meters;

[0091] Secondly, the result obtained by two-dimensional convolution using the point spread function is used as the conventional migration imaging result, such as Figure 3 As shown;

[0092] Next, use Figure 9 The two-dimensional Gaussian function shown is used as the initial point spread function. Deconvolution processing is performed based on the initial point spread function and the conventional migration result to obtain the deconvolution imaging result, as shown in FIG. Figure 4 As shown; and based on the deconvolution imaging results and the conventional migration imaging results, the point spread function is updated;

[0093] Finally, after repeating the above steps 10 times, the convergence state is reached and the high-resolution imaging result is obtained, as shown in Figure 5.

[0094] pass Figure 6 The comparison chart of the initial point spread function, the actual point spread function and the last updated point spread function is shown in the figure. Figures 7-9 As can be seen from the schematic diagram, the method of the present invention can correct the errors in the point spread calculation results caused by unknown source wavelets, inaccurate background velocity and forward / migration operators, avoid erroneous structural imaging caused by erroneous point spread functions, and improve the resolution of imaging results.

[0095] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit the present invention. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solutions of the present invention may be modified or replaced by equivalents without departing from the spirit and scope of the technical solutions of the present invention, which should all be included in the scope of the claims of the present invention.

Claims

1. A method for high-resolution processing of seismic data based on point spread function, characterized in that: include: Inputting seismic records, conventional migration results and migration velocity, and acquiring an observation system based on trace header information of the seismic records; Analytically calculate the initial point spread function based on the travel time and amplitude of the source wavelet and spatial position; performing deconvolution processing on the migration result based on the point spread function to obtain a deconvolved imaging result; Based on the deconvolution imaging result and the conventional migration result, updating the point spread function, and finally outputting the image domain least squares migration imaging result and the updated point spread function; The analytical formula of the point spread function is expressed as: in, is a spatial point x i The point spread function at R ss represents the autocorrelation of the source wavelet s(ω); T represents the travel time, A represents the amplitude value, and the amplitude of the scattered ray path is given by connecting the two incident rays: A(x r ;x j ;x s )=A(x j ,x r )A(x j ,x s ); The analytical formula of the point spread function includes: To address the computational complexity of the Hessian matrix, we introduce the asymptotic Green's function under the WKBJ approximation, and derive an approximate expression for the Hessian matrix, namely the analytical formula for the point spread function, which is expressed as the integral of the ray amplitude and the autocorrelation of the source wavelet. This allows us to calculate the point spread function for any point underground, simply by using the travel time, amplitude, and source wavelet. The deconvolution process comprises: By considering the sparsity and structural continuity of the reflection coefficient, L1 constraint and total variation constraint are introduced to construct the error functional J corresponding to the deconvolution process. r (r), and then realize deconvolution processing; The error functional J corresponding to the deconvolution process r (r) is expressed as: Among them, * represents the convolution operator, F (i) is the point spread function of the i-th iteration, r represents the deconvolution in the current iteration, λ1, λ2 and λ3 represent hyperparameters, ||r|| p 、||r|| TV and Δr represent the Lp constraint, total variation constraint and second-order derivative constraint introduced on the deconvolution result, respectively, p∈[0,2].

2. The method for high-resolution processing of seismic data based on point spread function according to claim 1, wherein: The input earthquake record includes Input earthquake records in SEGY and SU formats.

3. The method for high-resolution processing of seismic data based on point spread function according to claim 1 or 2, characterized in that: The acquisition observation system includes: Obtain the parameters of the observation system, namely the spatial position of the shot point and the spatial position of the receiver point.

4. The method for high-resolution processing of seismic data based on point spread function according to claim 3, characterized in that: The source wavelet includes: It is assumed that the general Ricker wavelet is used as the source wavelet parameter and substituted into the analytical formula of the point spread function to calculate the initial point spread function.

5. The method for high-resolution processing of seismic data based on point spread function according to claim 4, characterized in that: The conventional migration imaging results include: The conventional migration imaging result is regarded as the convolution of the point spread function and the reflection coefficient. For the one-dimensional case, the convolution operation is written as the Toeplitz matrix product form. For the two-dimensional or three-dimensional case, the convolution operation and the correlation operation are regarded as a pair of conjugate operators to perform gradient updates. Among them, the correlation operation specifically calculates the correlation coefficient between two variables.

6. The method for high-resolution processing of seismic data based on point spread function according to claim 5, characterized in that: The updating point spread function includes: Introducing the smoothness constraint μ1||F|| of the point spread function q , construct the error functional J F (F) to update the point spread function; The error functional J at this time F (F) is expressed as: in, represents the correlation operator, r (i) is the deconvolution result of the i-th iteration, ||F|| q is the Lq constraint introduced for the point spread function, μ1 is its corresponding hyperparameter, q∈(1,2].

7. The method for high-resolution processing of seismic data based on point spread function according to claim 6, characterized in that: The output image domain least squares migration imaging result and the updated point spread function include: Determine whether the error between the convolution of the point spread function and the deconvolution result and the conventional migration imaging result is less than a given threshold; when the error between the convolution of the point spread function and the deconvolution result and the conventional migration imaging result is less than the given threshold, output the image domain least squares migration imaging result and the updated point spread function.

Images