Seismic imaging point spread function calculation method based on space-time variable wavelets

Abstract

The invention provides a seismic imaging point spread function calculation method based on space-time variable wavelets, and the method comprises the steps: 1, inputting seismic data, and extracting seismic wavelets in different time periods; step 2, respectively calculating self-correlation for the seismic sub-bands, and calculating a fourth-order derivative; 3, point spread function calculation based on space-time variable wavelets is carried out, and ray tracing is carried out to form a shot point travel time field; 4, inputting new data, and circulating the steps 1-3 until all seismic data are calculated; and 5, carrying out deconvolution operation to obtain a high-resolution fidelity imaging profile. According to the seismic imaging point spread function calculation method based on the space-time variable wavelets, the wave number spectrum distribution of a final calculation point spread function result is improved, the consistency of the wave number spectrum and an actual imaging result is improved, the stability of final deconvolution can be improved, and the effect of least square migration of an imaging domain is improved.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic exploration, and in particular to a method for calculating a seismic imaging point spread function based on a spatio-temporal variable wavelet. Background Art

[0002] Migration imaging is one of the core technologies for seismic data processing in oil and gas exploration. The goal of traditional migration techniques is to locate subsurface structures, while inversion imaging represented by least-squares migration can better estimate the true reflection coefficients underground. The monograph "Seismic Inversion" by Schuster in 2017 details the method principle of least-squares migration and the challenges in method theory and implementation. In 2019, Schuster pointed out in the SEG annual meeting paper "Least-squares migration: Current and future directions" that the computational amount of least-squares migration in both the data domain and the imaging domain is very large, especially in the implementation of the data domain, where the computational amount is 10 to 30 times, and there are also convergence problems. Ru-Shan Wu et al. in the 2005 article "Directional illumination analysis using beamlet decomposition and propagation" in Geophysics and Xiao-Bi Xie et al. in the 2006 article "Wave-equation based seismic illumination analysis" in Geophysics pointed out that from the perspective of imaging correction of illumination, the true reflection coefficients underground can also be obtained, and it is relatively stable. The difficulty lies in the resolution function that characterizes the underground illumination and propagation factors, that is, the computational amount of the point spread function is relatively large and the consistency with the data is not strong enough. Subsequently, Chen and Xie proposed an algorithm for efficiently calculating the two-dimensional point spread function in the 2015 SEG annual meeting paper "An efficient method for broadband seismic illumination and resolution analyses". Wang Xuejun and Ren Haoran et al. continued to perform two-dimensional viscoacoustic medium inversion imaging by calculating the point spread function in the 2019 paper "Inversion imaging of viscoacoustic media based on point spread function" in Geophysical Prospecting for Petroleum.

[0003] At present, two-dimensional and three-dimensional migration imaging technologies are relatively mature and can be applied to industrial production. Calculating the point spread function and further performing deconvolution to improve the resolution of the imaging profile are areas that have been studied extensively in recent years. Almost all current theories for calculating the point spread function are based on known wavelets. However, in practical applications, wavelets exhibit spatio-temporal variations, which results in a mismatch between the calculated point spread function and the actual requirements, thereby affecting the final application effect. How to improve the matching of the point spread function with actual data and enhance the final imaging effect still requires in-depth research.

[0004] All of the above existing technologies are significantly different from the present invention and fail to solve the technical problems we aim to address. Therefore, we have invented a new method for calculating the seismic imaging point spread function based on spatio-temporal variable wavelets. Summary of the Invention

[0005] The objective of the present invention is to provide a method for calculating the seismic imaging point spread function based on spatio-temporal variable wavelets, which can improve the stability of the final deconvolution and enhance the effect of least-squares migration in the imaging domain.

[0006] The objective of the present invention can be achieved by the following technical measures: A method for calculating the seismic imaging point spread function based on spatio-temporal variable wavelets, which includes:

[0007] Step 1: Input a single seismic data trace and extract seismic wavelets at different time intervals.

[0008] Step 2: Calculate the autocorrelation for each seismic wavelet band and compute the fourth-order derivative.

[0009] Step 3: Perform the calculation of the point spread function based on spatio-temporal variable wavelets, and ray tracing to form the travel-time field between shot and receiver points.

[0010] Step 4: Input a new single data trace and loop through Steps 1 - 3 until all seismic data has been calculated.

[0011] Step 5: Perform deconvolution operations to obtain a high-resolution and fidelity imaging profile.

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

[0013] In Step 1, when inputting a single seismic data trace, extract seismic wavelets at different time intervals according to the parameters for calculating the point spread function, and store the seismic wavelets to form a set of wavelet trace groups for different time intervals.

[0014] In Step 1, change the conventional fixed given seismic wavelet to a spatio-temporal variable seismic wavelet to improve the matching of the subsequent calculated point spread function with real data and enhance the quality of the point spread function.

[0015] In step 3, taking advantage of the characteristic of Kirchhoff migration implemented trace by trace, seismic wavelets at different time intervals are directly extracted from the input single-trace seismic data to generate a time-varying wavelet trace combination. Using this combination to calculate the point spread function can minimize the difference between the wavenumber spectrum of the point spread function and the wavenumber spectrum of the actual data caused by the difference in the frequency bands of the wavelet and the actual data, and improve the stability and inversion effect of subsequent point spread function deconvolution.

[0016] In step 3, in reflection seismic exploration, the process of seismic data acquisition or forward modeling is as follows:

[0017] Lm = d obs (1)

[0018] where L represents the wave propagation operator, m represents the underground reflection parameter, d represents the acquired data, and the superscript obs refers to the observation, that is, the observed data. Mathematically, the exact solution of the underground reflection parameter m can be obtained by solving the inverse of the propagation operator. Usually, the inverse of L does not exist. To estimate the underground reflection parameter, migration is introduced in exploration seismology. The observed data is directly migrated and backpropagated to obtain the structural information of the underground reflection parameter. The process is as follows

[0019]

[0020] where the imaging result of the underground reflection parameter structure, L T represents the migration operator, that is, the transpose of the forward propagation operator. However, the above classical migration process can only obtain a relatively blurred reflection structure of the underground parameters. More accurate estimation of the underground parameters can be obtained through iterative inversion least-squares migration. The objective functional for inversion imaging in the least-squares sense can be written as

[0021]

[0022] The above objective function can be solved by the gradient-guided method. The normal equation of the above formula is as follows

[0023]

[0024] where H = L T L is the second derivative of the objective function in formula (3), that is, the Hessian matrix of the linear operator L. It can be seen from the above formula that the result of traditional imaging is equivalent to a blurred result of the Hessian matrix acting on the true underground reflection coefficient. The least-squares migration in the imaging domain attempts to directly solve the normal equation (4), estimate the inverse of the Hessian matrix, and act on the result of conventional imaging to obtain the inversion imaging result

[0025]

[0026] where The comparison after Hessian correction has a higher-precision least-squares migration result.

[0027] In step 3, under the Born approximation, the forward operator in equation (1) is written out specifically, and its form is as follows

[0028] d obs (x r , x s , ω) = ∫ω 2 G(x, x r , ω)G(x, x s , ω)s(ω)m(x)dV (6)

[0029] In the above formula, d obs (x r , x s , ω) represents the frequency domain, the surface source-receiver pair is located at (x r , x s ) for seismic acquisition (forward modeling) data, s(ω) represents the wavelet, and G(x, x s , ω) and G(x, x r , ω) are the Green's functions of the source point and the receiver point respectively; substituting the above formula into the Hessian matrix, H = L T L, we can obtain

[0030]

[0031] where ξ = (x s , x r ) represents the surface source-receiver pair, x i and x j both represent a certain underground grid point. Based on the asymptotic approximation theory, the Green's function under the WKBJ approximation can be expressed as

[0032]

[0033] where, T(x, x s ) represents the travel time term from the source (receiver) point to the underground imaging point, and A(x, x s ) represents the corresponding amplitude term; substituting equation (8) into equation (7) can obtain the expression of the Hessian matrix as follows

[0034]

[0035] where, R ss is the cross-correlation function of the wavelets; in the above formula, x i and x jThe dimensions are all in the model space. It can be seen that the size of the Hessian matrix is the square of the model. However, the Hessian matrix is sparse and banded, showing the characteristic that its elements have local support, that is, for a single underground point x i , only the surrounding x j is effective; according to its locality, for any underground point x i , only calculate the surrounding x j within a limited range as the point spread function of this point.

[0036] In step 3, the above derivation of equation (8) ignores the frequency attenuation caused by the inelastic effect of the medium and does not consider the wavelet differences caused by factors such as source excitation and geophone coupling during the actual observation process. Therefore, in practical applications, equation (9) needs to be modified as necessary to consider the wavelet differences at different time periods. Equation (9) is transformed into

[0037]

[0038] In equation (10), R ss (T(x r ; x j ; x s ) + T(x r ; x i ; x s )) is the cross-correlation of the wavelets corresponding to the moments of T(x r ; x j ; x s ) + T(x r ; x i ; x s ); in this way, the wavelets corresponding to different moments are different, and the time-varying characteristics of the point spread function are realized; in actual calculations, the wavelets at different moments need to be extracted separately for each seismic data trace, that is, the calculation of the point spread function based on the spatio-temporal variable wavelet is realized.

[0039] In step 3, ray tracing forms the travel-time field between shot points and receivers, and calculates the H corresponding to different points according to formula (10) to obtain the contributions of the point spread function corresponding to different traces, which are accumulated into the final point spread function volume; the conventional calculation formula is extended in a time-varying manner to obtain formula (10), so that a time-varying wavelet can be input, and combined with the time-varying wavelets extracted trace by trace in step 1, the point spread function based on the spatio-temporal variable wavelet is obtained.

[0040] In step 5, for the finally calculated point spread function, combined with the imaging profile obtained by conventional imaging, deconvolution operation can be performed to obtain a high-resolution and fidelity imaging profile.

[0041] The object of the present invention can also be achieved by the following technical measures: a seismic imaging point spread function calculation system based on space-time variant wavelets. This seismic imaging point spread function calculation system based on space-time variant wavelets uses a seismic imaging point spread function calculation method based on space-time variant wavelets to obtain a point spread function that better matches seismic data for deconvolution operations to obtain a high-resolution and fidelity imaging profile.

[0042] The seismic imaging point spread function calculation method based on space-time variant wavelets in the present invention can be applied to imaging amplitude correction and inversion in seismic exploration data processing. In the present invention, seismic wavelets corresponding to different time periods are extracted from each seismic data trace, fully considering the frequency changes of signals in different time periods caused by inelastic factors of underground media and frequency changes caused by seismic sources and geophones, etc., improving the wavenumber spectrum distribution of the final calculated point spread function result, enhancing the consistency between the wavenumber spectrum and the actual imaging result, thereby being able to improve the stability of the final deconvolution and the effect of least squares migration in the imaging domain. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] Figure 1 Schematic diagram of the basic concept of the point spread function in a specific embodiment of the present invention;

[0044] Figure 2 Schematic diagram of seismic wavelets in different time periods extracted from one of the data traces in a specific embodiment of the present invention;

[0045] Figure 3 Schematic diagram of the point spread function calculated using space-time variant wavelets in a specific embodiment of the present invention;

[0046] Figure 4 Schematic diagram of the conventional imaging profile of a certain data in a specific embodiment of the present invention;

[0047] Figure 5 Schematic diagram of the imaging profile after deconvolution of the point spread function calculated using this patent based on the conventional imaging profile of a certain data in a specific embodiment of the present invention;

[0048] Figure 6 In a specific embodiment of the present invention Figure 4 and Figure 5 Schematic diagram of the spectrum comparison corresponding to the imaging body;

[0049] Figure 7 Flowchart of a specific embodiment of the seismic imaging point spread function calculation method based on space-time variant wavelets of the present invention;

[0050] Figure 8 Schematic diagram of the point spread function calculation result in a specific embodiment of the present invention;

[0051] Figure 9 Schematic diagram of a conventional imaging profile in a specific embodiment of the present invention;

[0052] Figure 10 In a specific embodiment of the present invention, for Figure 9 Schematic diagram of the imaging profile after deconvolving the point spread function;

[0053] Figure 11 Schematic diagram of the sigsbee model velocity model for calculating the point spread function in a specific embodiment of the present invention;

[0054] Figure 12 Schematic diagram of the spatio-temporal variable point spread function calculated based on the sigsbee model in a specific embodiment of the present invention. Detailed implementation manners

[0055] It should be noted that the following detailed description is exemplary and is intended to provide further illustration of the present invention. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by those of ordinary skill in the technical field to which the present invention belongs.

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

[0057] The method for calculating the seismic imaging point spread function based on spatio-temporal variable wavelets in the present invention includes:

[0058] Step 1: Input a seismic data trace, extract seismic wavelets in different time periods according to the parameters for calculating the point spread function, and store the seismic wavelets to form a set of wavelet group traces in different time periods. Changing from the conventional fixed seismic wavelet to a spatio-temporal variable seismic wavelet can improve the matching between the subsequent calculation result of the point spread function and the real data, and improve the quality of the point spread function.

[0059] Step 2: Calculate the autocorrelation for each seismic wavelet band respectively, and calculate the fourth-order derivative.

[0060] Step 3: Perform the calculation of the point spread function based on spatio-temporal variable wavelets, and ray tracing forms the travel time field between shot points and receivers.

[0061] The seismic imaging point spread function is a locally supported Hessian matrix. The consistency between the calculated result of the point spread function and the vertical and horizontal wave numbers of the actual data is one of the key factors affecting the final inversion result. In view of the characteristics that the high-frequency information of seismic data gradually decays with the increase of propagation time, and there are usually frequency and amplitude differences in seismic traces at different spatial positions, the present invention proposes a method for calculating the point spread function based on a spatio-temporal variable wavelet. This method utilizes the characteristic of Kirchhoff migration implemented trace by trace, directly extracts seismic wavelets at different time periods from the input single-trace seismic data, generates a time-varying wavelet trace combination, and uses this combination to calculate the point spread function, minimizing the difference between the wave number spectrum of the point spread function and the wave number spectrum of the actual data caused by the frequency band difference between the wavelet and the actual data, and improving the stability and inversion effect of subsequent point spread function deconvolution.

[0062] In the classical seismic data acquisition and imaging system, the Hessian matrix is the system response including the acquisition (forward propagation of seismic waves) and imaging (backward propagation of seismic waves) processes. The point spread function is a locally supported Hessian matrix and can be regarded as the impulse response of the seismic data acquisition and imaging system. The core problem of seismic data inversion imaging is to correct the system effects caused by observation and propagation and image the underground parameters faithfully. In reflection seismic exploration, the process of seismic data acquisition or forward modeling is as follows:

[0063] Lm = d obs (1)

[0064] where L represents the wave propagation operator, m represents the underground reflection parameter, and d represents the acquired data. Mathematically, the exact solution of the underground reflection parameter m can be obtained by solving the inverse of the propagation operator. Usually, the inverse of L does not exist. To estimate the underground reflection parameter, migration is introduced in exploration seismology, and the observed data is directly migrated backward to obtain the structural information of the underground reflection parameter. The process is as follows

[0065]

[0066] where, The imaging result of the structure of the underground reflection parameter, L T represents the migration operator (the transpose of the forward propagation operator). However, the above classical migration process can only obtain a relatively blurred reflection structure of the underground parameters. A more accurate estimation of the underground parameters can be obtained through iterative inversion least-squares migration. The objective functional for inversion imaging in the least-squares sense can be written as

[0067]

[0068] Generally, the above objective function can be solved by a gradient-guided method. The normal equation of the above formula is as follows

[0069]

[0070] Among them, H = L T L is the second derivative of the objective function in Equation (3), that is, the Hessian matrix of the linear operator L. It can be seen from the above formula that the result of traditional imaging is equivalent to a blurred result of the Hessian matrix acting on the true underground reflection coefficient. The least-squares migration in the imaging domain attempts to directly solve the normal equation (4), estimate the inverse of the Hessian matrix, and act on the conventional imaging result to obtain the inversion imaging result

[0071]

[0072] Among them Compared with the Hessian correction There is a least-squares migration result with higher accuracy

[0073] Under the Born approximation, expand the forward operator in Equation (1) specifically, and its form is as follows

[0074] d obs (x r , x s , ω) = ∫ω 2 G(x, x r , ω)G(x, x s , ω)s(ω)m(x)dV (6)

[0075] In the above formula, d obs (x r , x s , ω) represents the frequency domain, the surface source-receiver pair is located at (x r , x s ) seismic acquisition (forward modeling) data, s(ω) represents the wavelet, and G(x, x s , ω) and G(x, x r , ω) are the Green's functions of the source point and the receiver point respectively. Substitute the above formula into the Hessian matrix, H = L T L, and we can get

[0076]

[0077] Among them, ξ = (x s , x r ) represents the surface source-receiver pair, x i and x j both represent a certain underground grid point. As shown in Figure 1 , where the box part is the travel time and amplitude calculation range of the corresponding stake point, and the center point of the point spread function and its neighborhood coordinates are x i and x j. Based on the asymptotic approximation theory, the Green's function under the WKBJ approximation can be expressed as

[0078]

[0079] where \(T(x,x s )\) represents the travel time term from the shot (geophone) point to the subsurface imaging point, and \(A(x,x s )\) represents the corresponding amplitude term. Substituting equation (8) into equation (7), the expression of the Hessian matrix can be obtained as follows

[0080]

[0081] where \(R ss \) is the cross-correlation function of the wavelets. In the above formula, the dimensions of \(x i \) and \(x j \) are both the model space. It can be seen that the size of the Hessian matrix is the square of the model. However, the Hessian matrix is sparse banded, which means that its elements have the property of local support. That is, for a single subsurface point \(x i \), only the surrounding \(x j \) is effective. According to its locality, for any subsurface point \(x i \), we only calculate the surrounding finite range of \(x j \) as the point spread function of this point.

[0082] The above derivation of equation (8) ignores the frequency attenuation caused by the inelastic effect of the medium and does not consider the wavelet differences caused by factors such as source excitation and geophone coupling in the actual observation process. Therefore, in practical applications, it is necessary to make necessary modifications to equation (9) to consider the wavelet differences at different time periods. Equation (9) is transformed into

[0083]

[0084] In equation (10), \(R ss (T(x r ; x j ; x s ) + T(x r ; x i ; x s ))\) is the cross-correlation of the wavelets corresponding to the time \(T(x r ; x j ; x s ) + T(x r ; x i ; x s ). In this way, the wavelets corresponding to different times are different, and the time-varying characteristics of the point spread function are realized. In actual calculations, it is necessary to extract the wavelets at different times for each seismic data trace, that is, the calculation of the point spread function based on the spatio-temporal variable wavelet is realized.

[0085] Ray tracing forms the travel time field between shot points and receivers, and according to formula (10), H corresponding to different points is calculated to obtain the contribution of the point spread function corresponding to different channels, which is accumulated into the final point spread function volume. This is also the core content of the present invention. The conventional calculation formula is extended in a time-varying manner to obtain formula (10), so that a time-varying wavelet can be input. By combining the time-varying wavelets extracted channel by channel in step S1, a point spread function based on spatio-temporal variable wavelets is obtained.

[0086] Step 4: Input a new channel of data, and loop through steps 1 - 3 until all seismic data is calculated.

[0087] Step 5: For the finally calculated point spread function, perform deconvolution operation in combination with the imaging profile obtained by conventional imaging to obtain a high-resolution and fidelity imaging profile.

[0088] The following are several specific embodiments of applying the present invention

[0089] Embodiment 1

[0090] As Figure 7 shown, Figure 7 is a flow chart of the method for calculating the point spread function of seismic imaging based on spatio-temporal variable wavelets of the present invention. The method for calculating the point spread function of seismic imaging based on spatio-temporal variable wavelets includes:

[0091] S1: Input a channel of seismic data. Since the final formula (10) requires extracting wavelets in different time periods, seismic wavelets in different time periods are extracted according to the parameters for calculating the point spread function, and the seismic wavelets are stored to form a set of wavelet group channels in different time periods; as Figure 2 shown, a series of sub-bands extracted from a certain channel of seismic data are presented, and the shape of the wavelet gradually changes from shallow to deep. This is also one of the supporting technical points of the present invention. Changing from the conventional fixed seismic wavelet to a spatio-temporal variable seismic wavelet can improve the matching between the subsequent calculation result of the point spread function and the real data, and improve the quality of the point spread function.

[0092] S2: Calculate the autocorrelation for each seismic sub-band respectively, and calculate the fourth-order derivative to implement the autocorrelation operation and the fourth-order derivative operation on the wavelet in formula (10).

[0093] S3: Ray tracing forms the travel time field between shot points and receivers, and according to formula (10), H corresponding to different points is calculated to obtain the contribution of the point spread function corresponding to different channels, which is accumulated into the final point spread function volume (as Figure 3 shown). This is also the core content of the present invention. The conventional calculation formula is extended in a time-varying manner to obtain formula (10), so that a time-varying wavelet can be input. By combining the time-varying wavelets extracted channel by channel in step S1, a point spread function based on spatio-temporal variable wavelets is obtained.

[0094] S4: Input a new set of data and loop through S1 - S3 until all seismic data calculations are completed.

[0095] For the finally calculated point spread function ( Figure 3 ), combined with the imaging profile obtained from conventional imaging ( Figure 4 ), performing deconvolution operation can obtain a high - resolution and fidelity imaging profile ( Figure 5 ). According to Figure 6 , where the black line is the spectrum of the conventional imaging profile and the gray curve is the spectrum after deconvolution of the point spread function of this patent, it can be seen that the high - resolution imaging profile has a wider frequency band than the conventional imaging.

[0096] Example 2

[0097] In the specific Example 2 of applying the present invention, Figure 8 is the result of the point spread function obtained by using the SEG three - dimensional research model data and the acquisition system with the technology of this patent. It can be seen that the energy of the point spread function gradually weakens from shallow to deep. In the area where seismic wave illumination is weak under the salt dome, the energy of the calculated point spread function is also relatively weak. Using the Figure 8 point spread function to perform deconvolution on the conventional imaging profile ( Figure 9 ) can obtain an imaging profile ( Figure 10 ) with improved resolution and fidelity. Comparing Figure 9 and Figure 10 , it can be seen that Figure 10 has more balanced energy and higher resolution.

[0098] Example 3

[0099] In the specific Example 3 of applying the present invention, Figure 11 uses the velocity model in the sedimentary area on the left side of the SEG sigsbee model as a test case, Figure 12 is the point spread function calculated by using the technology of this patent. It can be clearly seen that the total lateral resolution of the shallow layer of the point spread function calculated by using this patent is higher and the point spread function is more focused. The deeper it goes, the larger the points of the point spread function and the lower the resolution, which corresponds to the resolution of true migration imaging.

[0100] Finally, it should be noted that: The above - mentioned are only the preferred embodiments of the present invention and are not used to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, for those skilled in the art, they can still modify the technical solutions recorded in the foregoing embodiments, or perform equivalent replacements for some of the technical features. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

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

Claims

1. A method for calculating the seismic imaging point spread function based on spatio-temporal variable wavelets, characterized in that, The method for calculating the seismic imaging point spread function based on spatio-temporal variable wavelets includes: Step 1: Input a seismic data trace, and extract seismic wavelets at different time intervals; Step 2: Calculate the autocorrelation for each seismic wavelet band respectively, and calculate the fourth-order derivative; Step 3: Perform the calculation of the point spread function based on spatio-temporal variable wavelets, and ray tracing to form the travel-time field between shot and receiver; Step 4: Input a new seismic data trace, and loop through Steps 1 - 3 until all seismic data is calculated; Step 5: Perform deconvolution operation to obtain a high-resolution and fidelity imaging profile.

2. The method for calculating the seismic imaging point spread function based on spatio-temporal variable wavelets according to claim 1, wherein In Step 1, input a seismic data trace, extract seismic wavelets at different time intervals according to the parameters for calculating the point spread function, and store the seismic wavelets to form a set of wavelet trace groups at different time intervals.

3. The method for calculating the seismic imaging point spread function based on spatio-temporal variable wavelets according to claim 2, characterized in that In Step 1, change the conventional fixed seismic wavelet to a spatio-temporal variable seismic wavelet to improve the matching between the subsequent calculated result of the point spread function and the real data, and improve the quality of the point spread function.

4. The method for calculating the seismic imaging point spread function based on spatio-temporal wavelets according to claim 1, wherein In Step 3, taking advantage of the characteristic of Kirchhoff migration implemented trace by trace, directly extract seismic wavelets at different time intervals from the input single seismic data trace, generate a time-varying wavelet trace combination, and use this combination to calculate the point spread function, minimizing the difference between the wavenumber spectrum of the point spread function and the wavenumber spectrum of the real data caused by the frequency difference between the wavelet and the actual data, and improving the stability and inversion effect of the subsequent deconvolution of the point spread function.

5. The method for calculating the seismic imaging point spread function based on spatio-temporal variable wavelets according to claim 4, characterized in that, In Step 3, in reflection seismic exploration, the seismic data acquisition or forward modeling process is as follows: Lm = d obs (1) Where, L represents the wave propagation operator, m represents the underground reflection parameter, d represents the acquired data, and the superscript obs refers to the observation, i.e., the observed data; mathematically, the exact solution of the underground reflection parameter m can be obtained by solving the inverse of the propagation operator; usually, the inverse of L does not exist. To estimate the underground reflection parameter, migration is introduced in exploration seismology, and the observed data is directly migrated and back-propagated to obtain the structural information of the underground reflection parameter; the process is as follows Among them, Structural imaging result of the underground reflection parameter, L T denotes the migration operator, i.e., the transpose of the forward propagation operator; however, the above classical migration process can only obtain a relatively blurred reflection structure of the underground parameters, and a more accurate underground parameter estimation can be obtained through the least-squares migration of iterative inversion; the objective functional for inversion imaging in the least-squares sense can be written as The above objective function can be solved by the gradient-guided method; the normal equation of the above formula is as follows where H = L T L is the second derivative of the objective function in Equation (3), which is the Hessian matrix of the linear operator L. It can be seen from the above equation that the result of traditional imaging is equivalent to a blurred result of the Hessian matrix acting on the true underground reflection coefficient. The least squares migration in the imaging domain attempts to directly solve the normal equation (4), estimate the inverse of the Hessian matrix, and apply it to the conventional imaging result to obtain the inversion imaging result Among them, The comparison after Hessian correction has a higher-precision least-squares migration result.

6. The method for calculating the seismic imaging point spread function based on spatio-temporal variable wavelets according to claim 5, wherein, In Step 3, under the Born approximation, expand the forward operator in Equation (1) specifically, and its form is as follows d obs (x r ,x s ,ω)=∫ω 2 G(x,x r ,ω)G(x,x s ,ω)s(ω)m(x)dV (6) In the above formula, d obs (x r , x s , ω) represents the frequency domain, the surface shot-receiver points are located at the seismic acquisition forward modeling data of (x r , x s ), s(ω) represents the wavelet, G(x, x s , ω) and G(x, x r , ω) are the Green's functions of the shot point and the receiver point respectively; substituting the above formula into the Hessian matrix, H = L T L, we can obtain where ξ=(x s ,x r ) represents the surface shot-receiver pair, x i and x j both represent a certain underground grid point. Based on the asymptotic approximation theory, the Green's function under the WKBJ approximation can be expressed as Among them, T(x, x s ) represents the travel time term from the shot-receiver point to the subsurface imaging point, and A(x, x s ) represents the corresponding amplitude term; substituting equation (8) into equation (7) gives the expression of the Hessian matrix as follows where R ss is the cross-correlation function of the wavelets; in the above equation, the dimensions of x i and x j are both the model space. It can be seen that the size of the Hessian matrix is the square of the model. However, the Hessian matrix is sparse and banded, which means that its elements have the property of local support. That is, for a single point x i in the subsurface, only the x j around it is effective. According to its locality, for any point x i in the subsurface, only the x j in a finite range around it is calculated as the point spread function of this point.

7. The method for calculating the seismic imaging point spread function based on spatio-temporal variable wavelets according to claim 6, characterized in that, In Step 3, the above derivation of Equation (8) ignores the frequency attenuation caused by the inelastic effect of the medium, and does not consider the wavelet differences caused by factors such as the source excitation and the coupling of the geophones during the actual observation process. Therefore, in practical applications, it is necessary to make necessary modifications to Equation (9), considering the wavelet differences at different time intervals, and Equation (9) is transformed into In the formula (10), R ss (T(x r ; x j ; x s ) + T(x r ; x i ; x s )) is the cross-correlation of wavelets corresponding to the moment of T(x r ; x j ; x s ) + T(x r ; x i ; x s ); in this way, the wavelets corresponding to different moments are different, and the time-varying characteristics of the point spread function are realized; in actual calculation, it is necessary to extract the wavelets at different moments for each seismic data trace, that is, the calculation of the point spread function based on the spatio-temporal variable wavelet is realized.

8. The method for calculating the seismic imaging point spread function based on spatio-temporal variable wavelets according to claim 7, characterized in that In Step 3, ray tracing forms the travel-time field between shot and receiver, and calculate H corresponding to different points according to formula (10), obtain the contribution of the point spread function corresponding to different traces, and accumulate it into the final point spread function volume; the conventional calculation formula is extended in a time-varying manner to obtain formula (10), so that a time-varying wavelet can be input, and combined with the time-varying wavelets extracted trace by trace in Step 1, to obtain the point spread function based on spatio-temporal variable wavelets.

9. The method for calculating the seismic imaging point spread function based on spatio-temporal variable wavelets according to claim 1, characterized in that, In Step 5, for the finally calculated point spread function, combined with the imaging profile obtained by conventional imaging, perform deconvolution operation to obtain a high-resolution and fidelity imaging profile.

10. A seismic imaging point spread function calculation system based on spatio-temporal variable wavelets, characterized in that The seismic imaging point spread function calculation system based on spatio-temporal variable wavelets uses the method for calculating the seismic imaging point spread function based on spatio-temporal variable wavelets described in any one of claims 1-9 to obtain a point spread function that better matches seismic data, so as to perform deconvolution operation and obtain a high-resolution and fidelity imaging profile.

Images