Calculation method of seismic migration point spread function based on local time and space optimization

Abstract

The present invention discloses a method for calculating a seismic migration point spread function based on local time and space optimization, which relates to the field of exploration geophysical technology, including: acquiring data; calculating amplitude and complex travel time; Fourier transform and filtering; calculating a point scattering data volume; calculating a frequency spectrum; filtering and performing an inverse Fourier transform to obtain a point scattering response; calculating a Green's function and offsetting the point scattering data volume to a local space near the scattering point to calculate a point spread function; obtaining a simplified point spread function; calculating and weighting the point spread function; calculating a seismic image using a traditional adjoint migration method and decomposing it into local images; and deconvolution to obtain a reconstructed deconvolved image. The present invention calculates Gaussian beam-based born modeling and migration on a coarse grid based on local rays, interpolates the point spread function onto a fine image grid, and applies a high-dimensional Gaussian function to attenuate artifacts far from the center of the point spread function, which can significantly improve computational efficiency and obtain deep, high-quality imaging results.

Description

Technical Field

[0001] The present invention relates to the field of exploration geophysical technology, and in particular to a method for calculating a seismic offset point spread function based on local time and space optimization. Background Art

[0002] Seismic imaging is an important tool for detecting oil and gas resources and studying deep Earth structure. Traditional imaging methods based on rays and wave equations extend observational data and apply appropriate imaging conditions to construct subsurface impedance interfaces. However, due to band-limited source and receiver wavefields, incomplete seismic acquisition, and irregular subsurface illumination, neither Kirchoff migration nor one-way wave equation migration, nor even reverse-time migration, can produce high-quality images. Least-squares migration improves the resolution of seismic images and reduces artifacts caused by incomplete acquisition by continuously iteratively fitting observed data to predicted seismic data and calculating a generalized inverse for the subsurface reflectivity model. However, the enormous computational effort required for forward modeling and migration limits its widespread application in large-scale imaging problems.

[0003] Currently, to solve the problem of large least squares computational complexity, a common approach is to improve the computational efficiency of least squares migration in the data domain. However, this approach has two major problems: (1) it is difficult to converge the matching error between the predicted data and the actual data to an accurate migration imaging result through a short iteration time; (2) with the current level of computing power, the practical application of least squares migration is still in its infancy. An alternative approach is to promote the implementation of non-iterative least squares migration technology in the imaging domain. Its essence is to more efficiently calculate the point spread function (PSF). However, if the PSF sampling is too coarse, it cannot capture fine structural morphology, and if the sampling is too fine, there will be overlap between the PSFs, resulting in low imaging accuracy and reduced signal-to-noise ratio of the imaging profile.

[0004] Therefore, in order to obtain high-resolution underground structures with less computing power, it is urgent to develop an efficient and robust method for calculating the seismic offset point expansion function to provide technical support for the successful strategic succession of my country's deep oil and gas. Summary of the Invention

[0005] Compared with conventional migration methods, least squares migration can improve imaging resolution and deep imaging amplitude, and can eliminate the impact of poor acquisition lighting, but its computational complexity is too large, making it difficult to generalize to the imaging of actual seismic data. To address this problem, the present invention proposes a method for calculating the point spread function of seismic migration based on local time and space optimization. This method first calculates the Gaussian beam born modeling and migration based on local rays on a coarse grid, then interpolates the point spread function (PSF) onto a fine image grid, and finally applies a high-dimensional Gaussian function to attenuate artifacts far from the center of the point spread function (PSF). The method of the present invention can significantly improve computational efficiency and obtain deep, high-quality imaging results, with good practical application prospects.

[0006] To achieve the above object, the present invention adopts the following technical solutions:

[0007] The method for calculating the seismic migration point spread function based on local time and space optimization includes the following steps:

[0008] Step S1: Obtain input data, which includes the source wavelet f(t), background velocity field v0, and source position m s and detection point position m r ;

[0009] Step S2: Calculate the distance from the earthquake source m s After passing through scattering point m, it reaches detection point m r The amplitude A(m r ;m;m s ) and complex travel time T(m r ;m;m s );

[0010] Step S3: Perform Fourier transform on the source wavelet f(t) to obtain F(w), and filter the source wavelet in the frequency domain to obtain the filtered source wavelet.

[0011] Step S4: Calculate the point scattering data volume d generated from the scattering point m ps (m r ;m;m s ,t);

[0012] Step S5: Calculate d ps (m r ;m;m s ,t)’s spectrum D ps (m r ;m;m s ,ω);

[0013] Step S6: Filter in the frequency domain and then perform inverse Fourier transform to obtain the filtered point scattering response.

[0014] Step S7: Calculate Green's function G(m; m0, ω) using Gaussian beam operator and transform the point scattering data volume into Offset to the local space near the scattering point m and calculate the corresponding point spread function P sf (m,h);

[0015] Step S8, obtaining a simplified point spread function PSF;

[0016] Step S9: Use bilinear interpolation to calculate the point spread function PSF on the fine image grid, and use a multidimensional Gaussian function to weight the point spread function PSF to obtain a weighted PSF image.

[0017] Step S10: Calculate the seismic image I using the traditional adjoint migration method and decompose it into local images I l (m,h);

[0018] Step S11: l (m,h) and Deconvolution, get the reconstructed deconvolution image I psfdecon .

[0019] Optionally, in step S2, the amplitude A (m r ;m;m s ) and complex travel time T(m r ;m;m s ) is:

[0020] A(m r ;m;m s )=A(m;m s )A(m;m r )

[0021] T(m r ;m;m s )=T(m;m s )+T(m;m r )

[0022] Where m s represents the spatial location of the earthquake source, m represents the spatial location of the underground scattering point, and m r Represents the spatial position of the detection point; T(m; m s ) and T(m;m r ) are Gaussian beams from m s and m r The complex travel time to reach m; A(m; m s ) and A(m;m r ) are Gaussian beams from m s and mr The amplitude reaching m, T(m r ;m;m s ) represents the distance from the earthquake source position m s Passing through underground scattering point m, arriving at detection point m r The complex travel time, A(m r ;m;m s ) represents the distance from the earthquake source position m s Passing through underground scattering point m, arriving at detection point m r Amplitude.

[0023] Optionally, in step S3, the filtered source wavelet The expression is:

[0024]

[0025] Where t represents time; ∫ is the integral symbol; ω represents angular frequency, T i represents the imaginary part of the complex-valued traveltime T; F(ω) represents the Fourier transform of the input source wavelet f(t), i represents the imaginary unit, and exp[] represents the exponential function.

[0026] Optionally, in step S4, the point scattering data volume d ps (m r ;m;m s ,t) is:

[0027]

[0028] Where v0 represents the background velocity field; p x 、p y 、p z is the component of the ray parameter; s represents the source; r represents the detector; δ() is the Kronecker function; Ω represents the target imaging area underground; m′ is the integration variable; Re[] represents the real part; T r (m r ;m;m s ) and T i (m r ;m;m s ) represent the complex travel time T(m r ;m;m s )'s hour and imaginary parts.

[0029] Optionally, in step S5, the spectrum D ps (m r ;m;m s ,ω) is:

[0030]

[0031] Optionally, in step S6, the filtered point scattering response The expression is:

[0032]

[0033] Where * represents complex conjugation.

[0034] Optionally, in step S7, the corresponding point spread function P sf The expression of (m,h) is:

[0035]

[0036] Where h is the offset vector of the underground, ∑ represents the sum symbol; T1 is the record length;

[0037] The expression of the Green's function G(m; m0, ω) is:

[0038]

[0039] Where det stands for the determinant, P and Q are 2*2 complex matrices calculated by dynamic ray tracing, τ is the travel time along the central ray, q represents the two-dimensional orthogonal coordinate perpendicular to the central ray, the superscript T represents the transpose, and the superscript -1 represents the inverse matrix.

[0040] Optionally, in step S8, the step of obtaining a simplified point spread function PSF specifically includes:

[0041] Calculate and store ps (m r ;m;m s ,t) and In the time window of two cycles of a seismic wavelet, the grid spacing is set to 5 to 8 times the grid spacing of the traditional adjoint migration, in order to fully capture the changes in the PSF associated with the model;

[0042] The simplified point spread function PSF expression is:

[0043]

[0044] where t0 is the reference travel time that can be calculated in advance using ray tracing, and Δt is the half-time window, which is set to two periods of the wavelet.

[0045] Optionally, in step S9, the weighted PSF image The expression is:

[0046]

[0047] Where σ is the standard deviation of the weighting function.

[0048] Optionally, in step S10, the local image I l The expression of (m,h) is:

[0049]

[0050] Where Δx, Δy, and Δz are the spacings between the PSF centers along different directions, respectively; I(m+h) represents the traditional adjoint migration result.

[0051] Optionally, in step S11, the reconstructed deconvolution image I psfdecon The expression is:

[0052]

[0053] Where, represents the multidimensional Fourier transform, Represents the corresponding inverse transform, ∈(m) is a spatially varying function designed to prevent division by 0, and is set to one thousandth of the maximum value of the wavenumber spectrum.

[0054] The beneficial effects of the present invention are:

[0055] 1. This paper proposes a method for calculating the point spread function (PSF) for seismic migration based on local temporal and spatial optimization. This method primarily addresses the computational complexity of data-domain least-squares migration. This method optimizes the Green's function using a Gaussian beam operator. This optimized Green's function is then used to calculate the point spread function (PSF) within a local time window. These local waveforms are then migrated to the local subsurface space, effectively calculating the PSF and obtaining high-quality images. Statistically, this method's computational time is less than half that of traditional data-domain least-squares migration.

[0056] 2. To address the limited bandwidth and uneven amplitude of conventional companion migration, this paper uses a local point spread function (PSF) to estimate an approximate Hessian and applies it to PSF deconvolution. Compared to conventional companion migration imaging, this method significantly improves image resolution while producing fewer imaging artifacts. PSF deconvolution also incorporates the diagonal Hessian effect, which helps compensate for geometric diffusion and uneven illumination. This method, based on local temporal and spatial optimization, solves the computationally expensive problem of least-squares migration while ensuring accurate imaging of deep and ultra-deep oil and gas targets. Ultimately, it advances the practical application of high-resolution seismic migration in large-scale 3D oil and gas exploration. BRIEF DESCRIPTION OF THE DRAWINGS

[0057] Figure 1 This is a flow chart of the method for calculating the seismic migration point spread function based on local time and space optimization of the present invention;

[0058] Figure 2 A typical Zhongyuan Fault model in the central exploration area of ​​my country is shown as an embodiment of the present invention;

[0059] Figure 3 Common shot gathers calculated using a staggered grid finite difference forward modeling method according to an embodiment of the present invention;

[0060] Figure 4 Point spread functions (PSFs) calculated on a coarse grid as shown in one embodiment of the present invention;

[0061] Figure 5 This is the result of conventional accompanying migration imaging;

[0062] Figure 6 The PSF deconvolution imaging result shown in one embodiment of the present invention;

[0063] Figure 7 is the wavenumber spectrum of the true reflectivity model;

[0064] Figure 8 is the wavenumber spectrum corresponding to the conventional adjoint migration result;

[0065] Figure 9 FIG. 4 is a wavenumber spectrum corresponding to the PSF deconvolution imaging result shown in an embodiment of the present invention. DETAILED DESCRIPTION

[0066] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention without making any creative efforts shall fall within the scope of protection of the present invention.

[0067] Seismic migration point spread function calculation method based on local time and space optimization, such as Figure 1 As shown, the following steps are included:

[0068] Step S1: Obtain input data, which includes the source wavelet f(t), background velocity field v0, and source position m s and detection point position m r .

[0069] Step S2: Calculate the distance from the earthquake source m s After passing through scattering point m, it reaches detection point m r The amplitude A(m r ;m;m s ) and complex travel time T(m r ;m;m s );

[0070] The amplitude A(m r ;m;m s ) and complex travel time T(m r ;m;m s ) is:

[0071] A(m r ;m;m s )=A(m;m s )A(m;m r )

[0072] T(m r ;m;m s )=T(m;m s )+T(m;m r )

[0073] Where m s represents the spatial location of the earthquake source, m represents the spatial location of the underground scattering point, and m r Represents the spatial position of the detection point; T(m; m s ) and T(m;m r ) are Gaussian beams from m s and m r The complex travel time to reach m; A(m; m s ) and A(m;m r ) are Gaussian beams from m s and m r The amplitude reaching m, T(m r ;m;m s ) represents the distance from the earthquake source position m s Passing through underground scattering point m, arriving at detection point m r The complex travel time, A(m r ;m;m s ) represents the distance from the earthquake source position m s Passing through underground scattering point m, arriving at detection point m r Amplitude.

[0074] Step S3: Perform Fourier transform on the source wavelet f(t) to obtain F(w), and filter the source wavelet in the frequency domain to obtain the filtered source wavelet.

[0075] The filtered source wavelet The expression is:

[0076]

[0077] Where t represents time; ∫ is the integral symbol; ω represents angular frequency, T irepresents the imaginary part of the complex-valued traveltime T; F(ω) represents the Fourier transform of the input source wavelet f(t), i represents the imaginary unit, and exp[] represents the exponential function.

[0078] Step S4: Calculate the point scattering data volume d generated from the scattering point m ps (m r ;m;m s ,t);

[0079] The point scattering data volume d ps (m r ;m;m s ,t) is:

[0080]

[0081] Where v0 represents the background velocity field; p x 、p y 、p z is the component of the ray parameter; s represents the source; r represents the detector; δ() is the Kronecker function; Ω represents the target imaging area underground; m′ is the integration variable; Re[] represents the real part; T r (m r ;m;m s ) and T i (m r ;m;m s ) represent the complex travel time T(m r ;m;m s )'s hour and imaginary parts.

[0082] Step S5: Calculate d ps (m r ;m;m s ,t)’s spectrum D ps (m r ;m;m s ,ω);

[0083] The spectrum D ps (m r ;m;m s ,ω) is:

[0084]

[0085] Step S6: Filter in the frequency domain and then perform inverse Fourier transform to obtain the filtered point scattering response.

[0086] The filtered point scattering response The expression is:

[0087]

[0088] Where * represents complex conjugation.

[0089] Step S7: Calculate Green's function G(m; m0, ω) using Gaussian beam operator and transform the point scattering data volume into Offset to the local space near the scattering point m and calculate the corresponding point spread function P sf (m,h);

[0090] The corresponding point spread function P sf The expression of (m,h) is:

[0091]

[0092] Where h is the offset vector of the underground, ∑ represents the sum symbol; T1 is the record length;

[0093] The expression of the Green's function G(m; m0, ω) is:

[0094]

[0095] Where det stands for the determinant, P and Q are 2*2 complex matrices calculated by dynamic ray tracing, τ is the travel time along the central ray, q represents the two-dimensional orthogonal coordinate perpendicular to the central ray, the superscript T represents the transpose, and the superscript -1 represents the inverse matrix.

[0096] Step S8, obtaining a simplified point spread function PSF;

[0097] The steps for obtaining a simplified point spread function (PSF) include:

[0098] Calculate and store ps (m r ;m;m s ,t) and In the time window of two cycles of a seismic wavelet, the grid spacing is set to 5 to 8 times that of the conventional adjoint migration grid spacing to fully capture the changes in the PSF associated with the model;

[0099] The simplified point spread function PSF expression is:

[0100]

[0101] where t0 is the reference travel time that can be calculated in advance using ray tracing, and Δt is the half-time window, which is set to two periods of the wavelet.

[0102] Step S9: Use bilinear interpolation to calculate the point spread function PSF on the fine image grid, and use a multidimensional Gaussian function to weight the point spread function PSF to obtain a weighted PSF image.

[0103] The weighted PSF image The expression is:

[0104]

[0105] Where σ is the standard deviation of the weighting function.

[0106] Step S10: Calculate the seismic image I using the traditional adjoint migration method and decompose it into local images I l (m,h);

[0107] The local image I l The expression of (m,h) is:

[0108]

[0109] Where Δx, Δy, and Δz are the spacings between the PSF centers along different directions, respectively; I(m+h) represents the traditional adjoint migration result.

[0110] Step S11: l (m,h) and Deconvolution, get the reconstructed deconvolution image I psfdecon ;

[0111] The reconstructed deconvolution image I psfdecon The expression is:

[0112]

[0113] Where, represents the multidimensional Fourier transform, Represents the corresponding inverse transform, ∈(m) is a spatially varying function designed to prevent division by 0, and is set to one thousandth of the maximum value of the wavenumber spectrum.

[0114] The method provided by the present invention is applied to the typical Zhongyuan Fault Model in the central exploration area of ​​my country. The model has 475 vertical points and 1741 horizontal points. The horizontal and vertical spacing of the grid is set to 8 meters. Figure 2 As shown, use Figure 3 The forward modeling method based on staggered grid finite difference is shown in Figure 1 to calculate the common shot gathers. The point spread functions (PSFs) are calculated on a coarse grid with a spacing of 5-8 times the traditional offset grid spacing. Figure 4 shown.

[0115] Figure 5 This is the result of traditional adjoint migration. Although the basic structure has been imaged by adjoint migration, the resolution in the deep part is relatively low and the amplitude is weak due to reasons such as limited frequency band and irregular illumination. Figure 6 It is the result of PSF deconvolution imaging. It uses the local PSF to estimate an approximate Hessian and applies it to PSF deconvolution. Compared with the traditional accompanying migration image, it can significantly improve the image resolution and produce fewer imaging artifacts. Figure 7 Compared with the wavenumber spectrum of the true reflectivity model shown in Figure 2, the conventional adjoint migration produces a band-limited wavenumber spectrum due to the limited frequency band of the source wavelet, as shown in Figure 2. Figure 8 As shown. Figure 9 As shown, the PSF deconvolution wavenumber spectrum simultaneously expands the low- and high-wavenumber portions of the wavenumber spectrum, producing a wavenumber spectrum similar to the true reflectivity model, significantly improving spatial resolution. Furthermore, PSF deconvolution incorporates the diagonal Hessian effect, compensating for geometric diffusion and uneven illumination, resulting in an amplitude comparable to the true reflectivity model. In summary, compared to conventional adjoint migration methods, PSF deconvolution can achieve higher-quality imaging results while significantly reducing computational complexity, saving computing resources, and enabling rapid imaging, providing technical support for deep, high-precision oil and gas exploration.

[0116] Of course, the above description is not a limitation of the present invention, and the present invention is not limited to the above examples. Changes, modifications, additions or substitutions made by technicians in this technical field within the essential scope of the present invention should also fall within the scope of protection of the present invention.

Claims

1. A method for calculating the seismic migration point spread function based on local time and space optimization, characterized in that: The steps include: Step S1: Obtain input data, which includes the source wavelet f(t), background velocity field v0, and source position m s and detection point position m r ; Step S2: Calculate the distance from the earthquake source m s After passing through scattering point m, it reaches detection point m r The amplitude A(m r ;m;m s ) and complex travel time T(m r ;m;m s ); Step S3: Perform Fourier transform on the source wavelet f(t) to obtain F(w), and filter the source wavelet in the frequency domain to obtain the filtered source wavelet. Step S4: Calculate the point scattering data volume d generated from the scattering point m ps (m r ;m;m s ,t); Step S5: Calculate d ps (m r ;m;m s ,t)’s spectrum D ps (m r ;m;m s ,ω); Step S6: Filter in the frequency domain and then perform inverse Fourier transform to obtain the filtered point scattering response. Step S7: Calculate Green's function G(m; m0, ω) using Gaussian beam operator and transform the point scattering data volume into Offset to the local space near the scattering point m and calculate the corresponding point spread function P sf (m,h); Step S8, obtaining a simplified point spread function PSF; Step S9: Use bilinear interpolation to calculate the point spread function PSF on the fine image grid, and use a multidimensional Gaussian function to weight the point spread function PSF to obtain a weighted PSF image. Step S10: Calculate the seismic image I using the traditional adjoint migration method and decompose it into local images I l (m,h); Step S11: l (m,h) and Deconvolution, get the reconstructed deconvolution image I psfdecon .

2. The method for calculating the seismic migration point spread function based on local time and space optimization according to claim 1, characterized in that: In step S2, the amplitude A (m r ;m;m s ) and complex travel time T(m r ;m;m s ) is: A(m r ;m;m s )=A(m;m s )A(m;m r ) T(m r ;m;m s )=T(m;m s )+T(m;m r ) Where m s represents the spatial location of the earthquake source, m represents the spatial location of the underground scattering point, and m r Represents the spatial position of the detection point; T(m; m s ) and T(m;m r ) are Gaussian beams from m s and m r The complex travel time to reach m; A(m; m s ) and A(m;m r ) are Gaussian beams from m s and m r The amplitude reaching m, T(m r ;m; m s ) represents the distance from the earthquake source position m s Passing through underground scattering point m, arriving at detection point m r The complex travel time, A(m r ;m; m s ) represents the distance from the earthquake source position m s Passing through underground scattering point m, arriving at detection point m r Amplitude.

3. The method for calculating the seismic migration point spread function based on local time and space optimization according to claim 1, characterized in that: In step S3, the filtered source wavelet The expression is: Where, t represents time; ∫ is the integral symbol; ω represents the angular frequency, T i represents the imaginary part of the complex-valued traveltime T; F(ω) represents the Fourier transform of the input source wavelet f(t), i represents the imaginary unit, and exp[] represents the exponential function.

4. The method for calculating the seismic migration point spread function based on local time and space optimization according to claim 1, wherein: In step S4, the point scattering data volume d ps (m r ;m;m s ,t) is: Where v0 represents the background velocity field; p x 、p y 、p z is the component of the ray parameter; s represents the source; r represents the detector; δ() is the Kronecker function; Ω represents the target imaging area underground; m′ is the integral variable; Re[] represents the real part; T r (m r ;m;m s ) and T i (m r ;m;m s ) represent the complex travel time T(m r ;m;m s )’s real and imaginary parts; In step S5, the spectrum D ps (m r ;m;m s ,ω) is:

5. The method for calculating the seismic migration point spread function based on local time and space optimization according to claim 1, wherein: In step S6, the filtered point scattering response The expression is: Where ω represents the angular frequency, i represents the imaginary unit, and T i represents the imaginary part of the complex-valued traveltime T; F(ω) represents the Fourier transform of the input source wavelet f(t), and * represents the complex conjugate.

6. The method for calculating the seismic migration point spread function based on local time and space optimization according to claim 1, characterized in that: In step S7, the corresponding point spread function P sf The expression of (m,h) is: Where h is the offset vector of the underground, v0 represents the background velocity field; p x 、p y 、p z is the component of the ray parameter; s represents the source; r represents the detector; Re[] represents the real part; ∑ represents the sum symbol; T1 is the record length; The expression of the Green's function G(m; m0, ω) is: Where det stands for determinant, P and Q are 2*2 complex matrices calculated by dynamic ray tracing, τ is the travel time along the central ray, q represents the two-dimensional orthogonal coordinate perpendicular to the central ray, superscript T represents transpose, superscript -1 represents inverse matrix, and p x 、p y 、p z are the components of the ray parameters, i represents the imaginary unit, and ω represents the angular frequency.

7. The method for calculating the seismic migration point spread function based on local time and space optimization according to claim 1, characterized in that: In step S8, the step of obtaining a simplified point spread function PSF specifically includes: Calculate and store ps (m r ;m;m s ,t) and In the time window of two cycles of a seismic wavelet, the grid spacing is set to 5 to 8 times that of the traditional adjoint migration grid spacing; The simplified point spread function PSF expression is: Where t0 is the reference travel time that can be calculated in advance using ray tracing; Δt is the half time window, which is set to two periods of the wavelet, and p x 、p y 、p z are the components of the ray parameters; s represents the source; r represents the detector; v0 represents the background velocity field, h is the underground offset vector, and Re[] represents the real part.

8. The method for calculating the seismic migration point spread function based on local time and space optimization according to claim 1, characterized in that: In step S9, the weighted PSF image The expression is: Where σ is the standard deviation of the weighting function, and h is the offset vector of the underground.

9. The method for calculating the seismic migration point spread function based on local time and space optimization according to claim 1, wherein: In step S10, the local image I l The expression of (m,h) is: where Δx, Δy, and Δz are the spacing between PSF centers along different directions, I(m+h) represents the traditional adjoint migration result, and h is the offset vector in the underground.

10. The method for calculating the seismic migration point spread function based on local time and space optimization according to claim 1, wherein: In step S11, the reconstructed deconvolution image I psfdecon The expression is: Where, represents the multidimensional Fourier transform, represents the corresponding inverse transform, ∈(m) is a spatially varying function designed to prevent division by zero and is set to one thousandth of the maximum value of the wavenumber spectrum. Δx, Δy, and Δz are the spacing between the PSF centers along different directions, respectively; σ is the standard deviation of the weighting function, and h is the underground offset vector.

Images