System and method for interpolating seismic data

Abstract

A system and method of interpolating seismic data is provided. The system and method form a plurality of pairwise Hankel tensors from acquired seismic data, and a respective pairwise Hankel tensor for each of a plurality of originally collected frequency slices, perform tensor completion on each of said pairwise Hankel tensors to recover a plurality of interpolated frequency slices, and combine said plurality of interpolated frequency slices with said originally collected frequency slices to form a set of trace data of a geographical area of interest.

I claim: 
     
       1. A seismic survey method that improves an image of an explored subsurface geological formation by mitigating missing or faulty data, the method comprising:
 generating seismic waves with seismic sources; 
 recording, with seismic receivers, seismic data upon detecting the seismic waves including ones traveling through the explored subsurface geological formation; 
 determining that one or more traces of the seismic data are missing, corrupt or noisy; 
 forming, with a computing device, a plurality of pairwise Hankel tensors from the acquired seismic data, so that a pairwise Hankel tensor corresponds to each of a plurality of originally collected frequency slices and has two orders for each spatial dimension of the originally collected frequency slices, a number of spatial dimensions being maxim four; 
 performing a tensor completion for at least one of the pairwise Hankel tensors, to obtain interpolated seismic data suitable to replace the one or more traces of the seismic data that are missing, corrupt or noisy; 
 combining the seismic data with the interpolated seismic data to obtain a complete set of traces; and 
 generating the image based on the complete set of traces, to estimate an oil and gas yielding potential of the explored subsurface geological formation. 
 
     
     
       2. The method of  claim 1 , wherein the plurality of pairwise Hankel tensors are formed by:
 converting each of said seismic traces from a t-x domain to an f-x domain to form the plurality of frequency slices with frequencies ranging from F 0 , of about 0 Hertz, to F N , about a Nyquist frequency, 
 each of said plurality of originally collected frequency slices includes a first number of spatial dimensions. 
 
     
     
       3. The method of  claim 2 , wherein for a two spatial dimension frequency slice, S, with row dimension i and column dimension j, a fourth order pairwise Hankel tensor T H(i,j,m,n)  is formed according to
     T   H(i,j,m,n)   =S ( i+j −1, m+n −1),
 
 wherein the pairwise Hankel tensor T H(i,j,m,n) , i is a row position and j is a column position within each of a sub-matrix within said pairwise Hankel tensor, T H(i,j,m,n) , and m is a column position and n is a row position of the sub-matrices for the pairwise Hankel tensor T H(i,j,m,n) . 
 
     
     
       4. The method of  claim 1 , wherein said forming of the plurality of pairwise Hankel tensors comprises:
 determining a number of spatial dimensions D of each of said plurality of originally collected frequency slices S; 
 determining a length of each of said spatial dimensions D, wherein for each dimension i of said plurality of frequency slices, i=1 . . . , D a dimension is L(i) is length; 
 specifying a pairwise Hankel tensor number of modes as 2D; 
 specifying a length of each of said pairwise Hankel tensor modes as L(i)/2+1 for a first and any odd number of said tensor modes, (L(i)+1)/2 for a second and any even number of said pairwise Hankel tensor modes, and further 
 wherein if any tensor mode length is a fraction said fractional tensor mode length is reduced to a next lowest whole number; and a value for any particular pairwise Hankel tensor element is
     S ( j   1   +j   2 −1 ,j   3   +j   4 −1 , . . . ,j   m   +j   n −1),
 
 
 
       wherein m=2D−1, n=2D, and j's range over all possible values for said pairwise Hankel tensor. 
     
     
       5. The method of  1 , wherein said performing of tensor completion comprises:
 determining a low rank tensor R that minimizes ∥Z(T H −R)∥ F , which is a Frobenius norm of said pairwise Hankel tensor TH, wherein ∥·∥ F  is the Frobenius norm of the tensor elements, and Z(·) is an operator that zeroes out all elements that are unknown in T H . 
 
     
     
       6. The method of  claim 5 , wherein said tensor R provides an approximation to said unknown tensor elements of said Hankel tensor T H . 
     
     
       7. The method of  claim 6 , wherein said determining of the low rank tensor R comprises:
 evaluating 
 
       
         
           
             
               
                 R 
                 = 
                 
                   
                     ∑ 
                     
                       i 
                       = 
                       1 
                     
                     k 
                   
                   ⁢ 
                   
                     
                       u 
                       1 
                       i 
                     
                     ∘ 
                     
                         
                     
                     ⁢ 
                     
                       u 
                       2 
                       i 
                     
                     ∘ 
                     … 
                     ∘ 
                     
                       u 
                       p 
                       i 
                     
                   
                 
               
               , 
             
           
         
         wherein said vectors u p   i , are determined using a series of linear least squares approximations based on an initial estimate of said vectors u p   i , and further wherein i is valued from 1 to k, k being said rank of said pairwise Hankel tensor T, and p is said order of said pairwise Hankel tensor T, said vectors iteratively changed to minimize the Frobenius norm of elements of said pairwise Hankel tensor T. 
       
     
     
       8. The method of  claim 7 , further comprising:
 recovering said interpolated frequency slices from said determined low rank tensor R, by averaging over every tensor element in which each frequency slice value was originally placed. 
 
     
     
       9. The method of  claim 1 , wherein said combining of said plurality of interpolated frequency slices with said originally collected frequency slices to form a complete set of traces comprises:
 performing an inverse discrete Fourier transform for each of said plurality of traces in the f-x domain to obtain interpolated trace data in said t-x domain; and 
 combining said interpolated traces in said t-x domain with original traces in said t-x domain that formed said original collected frequency slice data. 
 
     
     
       10. A seismic survey system that improves an image of an explored subsurface geological formation by mitigating missing or faulty data, comprising:
 a receiver configured to acquire transmitted seismic data detecting seismic waves with seismic sources after the seismic waves travelled through the explored subsurface geological formation; and 
 a processor configured to,
 determine that one or more traces of the seismic data are missing, corrupt or noisy; 
 form a plurality of pairwise Hankel tensors from said acquired seismic data, a respective pairwise Hankel tensor being formed for each of a plurality of originally collected frequency slices and each of the pairwise Hankel tensors having two orders for each spatial dimension of the originally collected frequency slices, a number of spatial dimensions being maxim four, 
 perform tensor completion at least one of said pairwise Hankel tensors to obtain interpolated seismic data suitable to replace the one or more traces of the seismic data that are missing, corrupt or noisy, 
 combine the seismic data with the interpolated seismic data to obtain a complete set of traces, and 
 generate the image of the explored subsurface geological formation based on the complete set of traces, for estimating an oil and gas yielding potential of the explored subsurface geological formation. 
 
 
     
     
       11. The system of  claim 10 , wherein said processor is further configured, when forming a plurality of pairwise Hankel tensors, to
 convert each of said seismic trace data from a t-x domain to an f-x domain to form the plurality of frequency slices with frequencies ranging from F 0 , of about 0 Hertz, to F N , about a Nyquist frequency, and 
 each of said plurality of originally collected frequency slices including a first number of spatial dimensions. 
 
     
     
       12. The system of  claim 11 , wherein, for a two spatial dimension frequency slice, S, with row dimension i and column dimension j, a fourth order pairwise Hankel tensor T H(i,j,m,n)  is formed according to
     T   H(i,j,m,n)   =S ( i+j −1, m+n −1),
 
 wherein the pairwise Hankel tensor T H(i,j,m,n) , i is a row position and j is a column position within each of a sub-matrix within said pairwise Hankel tensor, T H(i,j,m,n) , and m is a column position and n is a row position of the sub-matrices for the pairwise Hankel tensor T H(i,j,m,n) . 
 
     
     
       13. The system of  claim 10 , wherein said processor is further configured to form the plurality of pairwise Hankel tensors by
 determining a number of spatial dimensions D of each of said plurality of originally collected frequency slices S; 
 determining a length of each of said spatial dimensions D, wherein for each dimension i of said plurality of frequency slices, i=1 . . . , D a dimension is L(i) is length; 
 specifying a pairwise Hankel tensor number of modes as 2D; 
 specifying a length of each of said pairwise Hankel tensor modes as L(i)/2+1 for a first and any odd number of said tensor modes, (L(i)+1)/2 for a second and any even number of said pairwise Hankel tensor modes, and further wherein if any tensor mode length is a fraction said fractional tensor mode length is reduced to a next lowest whole number; and a value for any particular pairwise Hankel tensor element is
     S ( j   1   +j   2 −1 ,j   3   +j   4 −1 , . . . ,j   m   +j   n −1),
 
 
 wherein m=2D−1, n=2D, and j's range over all possible values for said pairwise Hankel tensor. 
 
     
     
       14. The system of  claim 10 , wherein said processor is further configured, when performing tensor completion, to
 determine a low rank tensor R that minimizes ∥Z(T H −R)∥ F , which is a Frobenius norm of said pairwise Hankel tensor TH, wherein ∥·∥ F  is the Frobenius norm of the tensor elements, and Z(·) is an operator that zeroes out all elements that are unknown in T H . 
 
     
     
       15. The system of  claim 14 , wherein said tensor R provides an approximation to said unknown tensor elements of said Hankel tensor T H . 
     
     
       16. The system of  claim 14 , wherein said processor is further configured, when determining a low rank tensor R, to evaluate 
       
         
           
             
               
                 R 
                 = 
                 
                   
                     ∑ 
                     
                       i 
                       = 
                       1 
                     
                     k 
                   
                   ⁢ 
                   
                     
                       u 
                       1 
                       i 
                     
                     ∘ 
                     
                         
                     
                     ⁢ 
                     
                       u 
                       2 
                       i 
                     
                     ∘ 
                     … 
                     ∘ 
                     
                       u 
                       p 
                       i 
                     
                   
                 
               
               , 
             
           
         
         wherein said vectors u p   i , are determined using a series of linear least squares approximations based on an initial estimate of said vectors u p   i , and further wherein i is valued from 1 to k, k being said rank of said pairwise Hankel tensor T, and p is said order of said pairwise Hankel tensor T, said vectors iteratively changed to minimize the Frobenius norm of elements of said pairwise Hankel tensor T. 
       
     
     
       17. The system of  claim 16 , wherein said processor is further configured to recover said interpolated frequency slices from said determined low rank tensor R, by averaging over every tensor element in which each frequency slice value was originally placed. 
     
     
       18. The system of  claim 10 , wherein said processor is further configured, when combining said the seismic data with the interpolated seismic data, to perform an inverse discrete Fourier transform for each of said plurality of interpolated traces in the f-x domain to obtain interpolated trace data in said t-x domain, and
 combine said interpolated traces in said t-x domain with original traces in said t-x domain that formed said original collected frequency slice data. 
 
     
     
       19. The system of  claim 10 , wherein the system is a land seismic system. 
     
     
       20. The system of  claim 10 , wherein the system is a marine seismic system.