Abstract:
A full wave inversion (FWI) may utilize an Amplitude-Frequency-Differentiation (AFD) or a Phase-Frequency-Differentiation (PFD) operation to form a velocity model of a subterranean formation utilizing recovered low wavenumber data. Received seismic data is processes to isolate two data signals at different frequencies. In an AFD operation, the two data signals are summed and the data of the envelope of the summed function is used for the FWI. In a PFD operation, the phase data of the quotient of the two data signals is used for the FWI. The FWI proceeds iteratively utilizing either the AFD or PFD data or with single frequency data until the cost function of the AFD or PFD is satisfied.

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS 
       [0001]    This application is a nonprovisional application which claims priority from U.S. provisional application No. 62/107,025, filed Jan. 23, 2015, the entirety of which is hereby incorporated by reference. 
     
    
     TECHNICAL FIELD/FIELD OF THE DISCLOSURE 
       [0002]    The present disclosure relates to processing of seismic data. 
       BACKGROUND OF THE DISCLOSURE 
       [0003]    In order to determine the composition and structure of an underground formation, seismic data may be utilized. As understood in the art, when a seismic wave propagates through the earth from a source to a receiver, the received seismic wave may bring information about geophysical properties of the subsurface volume. The properties may include, without limitation, seismic wave velocity, mass density, and anisotropy properties. In exploration geophysics, by recording the seismic waves at the surface by one or more geophones or hydrophones, a model of the subsurface strata may be determined through various seismic processing techniques. 
         [0004]    The time between transmission of the seismic wave and reception of the reflected seismic waves may be measured to determine the distance to an interface. More advanced techniques, such as reflection velocity tomography, are capable of using kinematic information contained in the seismic reflection data to transform the seismic data gathered from the reflected waves to determine a more detailed and complete model of the underground seismic velocity model. Another technique, known as full waveform inversion (FWI), utilizes full waveform information, i.e. both phase and amplitude data, both transmission seismic data and reflection seismic data, both primaries and multiples, contained in the seismic data gathered from the recorded seismic waves to generate a high resolution velocity model or other geophysical property model of the underground formation through optimal data fitting. 
         [0005]    Because of hardware limitations of geophones and hydrophones, seismic data is typically band limited, rendering low frequency seismic data unavailable. Due to the lack of low frequency data, large underground structures are less able to be identified due to issues including without limitation nonlinearity of the FWI and the cycle-skipping phenomenon. Nonlinearity and cycle-skipping can result in incorrect data fitting as, for example, the inversion may be trapped in local minima, preventing the updated models from progressing properly. Nonlinearity and cycle-skipping each increase in effect and likelihood with increasing frequency during an inversion operation. 
       SUMMARY 
       [0006]    The present disclosure provides for a method for generating a model of an underground formation from seismic data. The method may include measuring seismic data from one or more receivers in the time domain. The one or more receivers may be positioned to receive seismic waves generated by a source and emitted into the underground formation, the seismic waves forming the seismic data. The method may further include converting the seismic data from the time domain to the frequency domain to form frequency domain seismic data. The method may further include extracting, from the frequency domain seismic data, first measured frequency domain data from corresponding to a first signal at a first frequency. The method may further include extracting, from the frequency domain seismic data, second measured frequency domain data corresponding to a second signal at a second frequency. The second frequency is different from the first frequency. The method may further include generating, from a first velocity model, first simulated frequency domain data at the first frequency and second simulated frequency domain data at the second frequency. The method may further include calculating a gradient of a cost function utilizing the first and second measured frequency domain data and first and second simulated frequency domain data. The method may further include generating a second velocity model based at least in part on the calculated gradient, the second velocity model corresponding to at least one feature of the underground formation. 
     
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         [0007]    The present disclosure is best understood from the following detailed description when read with the accompanying figures. It is emphasized that, in accordance with the standard practice in the industry, various features are not drawn to scale. In fact, the dimensions of the various features may be arbitrarily increased or reduced for clarity of discussion. 
           [0008]      FIGS. 1A, 1B  depict a beat signal of an Amplitude-Frequency-Differentiation (AFD) method FWI consistent with embodiments of the present disclosure as created by two higher frequency signals. 
           [0009]      FIGS. 2A, 2B  depict a beat signal of an Amplitude-Frequency-Differentiation (AFD) method FWI consistent with embodiments of the present disclosure as received. 
           [0010]      FIGS. 3B, 3C  depict a beat signal of an Amplitude-Frequency-Differentiation (AFD) method FWI consistent with embodiments of the present disclosure as received showing velocity anomalies as depicted in  FIG. 3A . 
           [0011]      FIGS. 4A, 4B  depict phase data of a beat signal of a Phase-Frequency-Differentiation (PFD) method FWI consistent with embodiments of the present disclosure as received. 
           [0012]      FIGS. 5B, 5C  depict phase data of a beat signal of a Phase-Frequency-Differentiation (PFD) method FWI consistent with embodiments of the present disclosure as received showing velocity anomalies as depicted in  FIG. 5A . 
           [0013]      FIGS. 6B, 6C  depict phase data of a beat signal of a Phase-Frequency-Differentiation (PFD) method FWI consistent with embodiments of the present disclosure as received showing velocity anomalies as depicted in  FIG. 6A . 
           [0014]      FIGS. 7C, 7D  depict iterations of a model generated by a Phase-Frequency-Differentiation FWI consistent with embodiments of the present disclosure from the model depicted in  FIG. 7A  using an initial model as shown in  FIG. 7B . 
           [0015]      FIG. 8  is a flow chart depicting an FWI consistent with embodiments of the present disclosure. 
       
    
    
     DETAILED DESCRIPTION 
       [0016]    It is to be understood that the following disclosure provides many different embodiments, or examples, for implementing different features of various embodiments. Specific examples of components and arrangements are described below to simplify the present disclosure. These are, of course, merely examples and are not intended to be limiting. In addition, the present disclosure may repeat reference numerals and/or letters in the various examples. This repetition is for the purpose of simplicity and clarity and does not in itself dictate a relationship between the various embodiments and/or configurations discussed. 
         [0017]    In the following disclosure, it is to be understood that frequency of a signal may be expressed as either ordinary frequency, commonly measured in cycles per second or Hz, or angular frequency, measured in radians per second. As used herein, f refers to ordinary frequency and co refers to angular frequency. One having ordinary skill in the art with the benefit of this disclosure will understand that the frequencies are related according to: 
         [0000]      ω=2π f  
 
         [0000]    and thus unless explicitly stated otherwise, the use of one or the other is not intended to be limiting. 
         [0018]    In some embodiments of the present disclosure, a seismic geologic survey is carried out utilizing a seismic source commonly referred to as a “shot” to generate seismic waves. As understood in the art, the shot may be band limited, but include seismic waves in a range of frequencies. In some embodiments, the shot may create seismic waves from, for example and without limitation, 0.2 to 100 Hz. In analyzing the seismic waves received by the one or more geophones or hydrophones, frequency domain data may be extracted from the received time domain seismic waves corresponding to two frequencies of seismic waves using various signal processing techniques. These signals are referred to herein as S 1  and S 2 . The frequencies of S 1  and S 2  are selected to be different but close to each other. In some embodiments, S 1  and S 2  may each be between 4 and 100 Hz, between 5 and 30 Hz, or between 5 and 15 Hz. In some embodiments, for example and without limitation, S 1  and S 2  may be separated by 5 Hz or less, 2 Hz or less, or 0.5 Hz or less. 
         [0019]    As understood in the art, when two time domain single frequency acoustic waves are close in frequency, the interference between the two signals may result in an apparently single tone of temporally varying volume or amplitude. The rate of variation in volume, known as a beat frequency, is determined by the difference between the two sound frequencies. For a first signal, S 1 (t)=cos(2πf 1 t) and a second signal S 2 (t)=cos(2πf 2 t), where |f 1 −f 2 |&lt;&lt;f 1  and f 2 , the summation of the two signals is given by: 
         [0000]    
       
         
           
             
               
                 S 
                 3 
               
                
               
                 ( 
                 t 
                 ) 
               
             
             = 
             
               
                 
                   
                     S 
                     1 
                   
                    
                   
                     ( 
                     t 
                     ) 
                   
                 
                 + 
                 
                   
                     S 
                     2 
                   
                    
                   
                     ( 
                     t 
                     ) 
                   
                 
               
               = 
               
                 2 
                  
                 
                   cos 
                    
                   
                     ( 
                     
                       2 
                        
                       π 
                        
                       
                         
                           
                             f 
                             1 
                           
                           + 
                           
                             f 
                             2 
                           
                         
                         2 
                       
                        
                       t 
                     
                     ) 
                   
                 
                  
                 
                   cos 
                    
                   
                     ( 
                     
                       2 
                        
                       π 
                        
                       
                         
                           
                             f 
                             2 
                           
                           - 
                           
                             f 
                             1 
                           
                         
                         2 
                       
                        
                       t 
                     
                     ) 
                   
                 
               
             
           
         
       
     
         [0020]    The addition of S 1 (t) and S 2 (t) creates a perceived carrier signal (defined by 
         [0000]    
       
         
           
             
               
                 
                   
                     S 
                     c 
                   
                    
                   
                     ( 
                     t 
                     ) 
                   
                 
                 = 
                 
                   2 
                    
                   
                     cos 
                      
                     
                       ( 
                       
                         2 
                          
                         π 
                          
                         
                           
                             
                               f 
                               1 
                             
                             + 
                             
                               f 
                               2 
                             
                           
                           2 
                         
                          
                         t 
                       
                       ) 
                     
                   
                 
               
               ) 
             
             , 
           
         
       
     
         [0000]    whose amplitude varies slowly with the beat frequency 
         [0000]    
       
         
           
             
               f 
               b 
             
             = 
             
               
                 
                   f 
                   2 
                 
                 - 
                 
                   f 
                   1 
                 
               
               2 
             
           
         
       
     
         [0000]    (according to the 
         [0000]    
       
         
           
             cos 
              
             
               ( 
               
                 2 
                  
                 π 
                  
                 
                   
                     
                       f 
                       2 
                     
                     - 
                     
                       f 
                       1 
                     
                   
                   2 
                 
                  
                 t 
               
               ) 
             
           
         
       
     
         [0000]    term). Thus, by selecting the difference in frequency between S 1 (t) and S 2 (t), the high frequency carrier signal S c (t), is modulated with the low frequency beat signal. In some embodiments, as an example and without limitation, for time domain single frequency signals, by taking the envelope of S 3 ( t ), referred to herein as beat signal S B (t), low frequency data may be extracted from the high-frequency signals S 1 (t) and S 2 (t). As depicted in  FIG. 1A , signals S 1 (t) and S 2 (t) when superimposed create the amplitude variations of beat signal S B (t) as shown in  FIG. 1B , and the frequency of beat signal S B (t) is determined by the difference in frequency between signals S 1 (t), S 2 (t). By selecting the frequencies of S 1 (t) and S 2 (t) close to one another, beat signal S B (t) may have a low frequency. 
         [0021]    As understood in the art, when two frequency domain seismic waves recorded at multiple spatial locations are close in frequency, the interference between the two signals results in an apparently single tone of spatially varying amplitude. The rate of variation in amplitude is determined by the difference between the two seismic wave frequencies or wavenumbers. As understood in the art, wavenumber refers to the spatial frequency of a wave, and may be calculated according to: 
         [0000]    
       
         
           
             k 
             = 
             
               ω 
               v 
             
           
         
       
     
         [0000]    where k is the wavenumber and v is the phase velocity of the wave. For a first signal, S 1 (x)=exp(ik 1 x) and a second signal S 2 (x)=exp(ik 2 x), where x is the spatial location, |f 1 −f 2 |&gt;&gt;f 1  and f 2 , 
         [0000]    
       
         
           
             
               
                 k 
                 1 
               
               = 
               
                 
                   ω 
                   1 
                 
                 v 
               
             
             , 
             
               
                 k 
                 2 
               
               = 
               
                 
                   ω 
                   2 
                 
                 v 
               
             
             , 
           
         
       
     
         [0000]    the summation of the two signals is given by: 
         [0000]    
       
         
           
             
               
                 S 
                 3 
               
                
               
                 ( 
                 x 
                 ) 
               
             
             = 
             
               
                 
                   
                     S 
                     1 
                   
                    
                   
                     ( 
                     x 
                     ) 
                   
                 
                 + 
                 
                   
                     S 
                     2 
                   
                    
                   
                     ( 
                     x 
                     ) 
                   
                 
               
               = 
               
                 
                   2 
                   
                     | 
                     x 
                     | 
                   
                 
                  
                 
                   exp 
                    
                   
                     ( 
                     
                       
                         
                           
                             k 
                             1 
                           
                           + 
                           
                             k 
                             2 
                           
                         
                         2 
                       
                        
                       x 
                     
                     ) 
                   
                 
                  
                 
                   cos 
                    
                   
                     ( 
                     
                       
                         
                           
                             k 
                             2 
                           
                           - 
                           
                             k 
                             1 
                           
                         
                         2 
                       
                        
                       x 
                     
                     ) 
                   
                 
               
             
           
         
       
     
         [0022]    The addition of S 1 (x) and S 2 (x) creates a perceived carrier signal (defined by 
         [0000]    
       
         
           
             
               
                 
                   
                     S 
                     c 
                   
                    
                   
                     ( 
                     x 
                     ) 
                   
                 
                 = 
                 
                   exp 
                    
                   
                     ( 
                     
                       
                         
                           
                             k 
                             1 
                           
                           + 
                           
                             k 
                             2 
                           
                         
                         2 
                       
                        
                       x 
                     
                     ) 
                   
                 
               
               ) 
             
             , 
           
         
       
     
         [0000]    whose amplitude varies slowly with the beat wavenumber 
         [0000]    
       
         
           
             
               k 
               b 
             
             = 
             
               
                 
                   k 
                   2 
                 
                 - 
                 
                   k 
                   1 
                 
               
               2 
             
           
         
       
     
         [0000]    (according to the 
         [0000]    
       
         
           
             
               
                 
                   cos 
                    
                   
                     ( 
                     
                       
                         
                           
                             k 
                             2 
                           
                           - 
                           
                             k 
                             1 
                           
                         
                         2 
                       
                        
                       x 
                     
                     ) 
                   
                 
                  
                 term 
               
               ) 
             
             . 
           
         
       
     
         [0000]    Thus, by selecting the difference in frequency between S 1 (x) and S 2 (x), the high wavenumber carrier signal S c (x) is modulated with the low wavenumber beat signal. In some embodiments, as an example and without limitation, for frequency domain signals, by taking the envelope of S 3 (x), referred to herein as beat signal S B (x), low wavenumber data may be extracted from the high-wavenumber signals S 1 (x) and S 2 (x). As depicted in  FIG. 2A , signals S 1 (x) and S 2 (x) when superimposed create the amplitude variations of beat signal S B (x) as shown in  FIG. 2B , and the wavenumber of beat signal S B (t) is determined by the difference in wavenumber or frequency between signals S 1 (x), S 2 (x). By selecting the frequencies of S 1 (x) and S 2 (x) close to one another, beat signal S B (x) may have a low wavenumber. 
         [0023]    By utilizing the data in a low wavenumber beat signal S B  for a FWI operation, otherwise unavailable low wavenumber data may be obtained. 
         [0024]    An FWI inversion begins with an initial model, referred to herein as KO. After iterative inversions of seismic data, the model is gradually updated by using determined gradient information. Acoustic waves propagate according to: 
         [0000]      ∇ 2   p+ω   2   ρκρ=−iωQ  
 
         [0000]    where p is the pressure, Q is the injected volume acting as the source, ρ is the mass density, and κ is the compressibility of the medium. The seismic P-wave velocity, v p  is determined according to: 
         [0000]        v   p =1/√{square root over (ρκ)}
 
         [0025]    The model is updated until a model κ n =κ 0 +Δκ which is close to the true model (referred to herein as κ true ) is determined, such that the simulated data closely fits the measured data. The gradient information is calculated from the sensitivity kernel, i.e. the response of the simulated data to the perturbation of the medium property, according to: 
         [0000]    
       
         
           
             
               
                 ∂ 
                 p 
               
               
                 ∂ 
                 κ 
               
             
             = 
             
               
                 - 
                 
                   ρω 
                   2 
                 
               
                
               
                 
                   ∫ 
                   τ 
                 
                  
                 
                   
                     p 
                      
                     
                       ( 
                       
                         ω 
                         , 
                         
                           κ 
                           ; 
                           x 
                         
                         , 
                         
                           x 
                           r 
                         
                       
                       ) 
                     
                   
                    
                   
                     G 
                      
                     
                       ( 
                       
                         ω 
                         , 
                         
                           κ 
                           ; 
                           x 
                         
                         , 
                         
                           x 
                           s 
                         
                       
                       ) 
                     
                   
                    
                   
                      
                     x 
                   
                 
               
             
           
         
       
     
         [0000]    where G(ω, κ; x, x s ) is the Green&#39;s function, i.e., the solution to 
         [0000]      ∇ 2   p+ω   2 ρκρ=δ( x−x   s )
 
         [0026]    Because 
         [0000]    
       
         
           
             
               ∂ 
               p 
             
             
               ∂ 
               κ 
             
           
         
       
     
         [0000]    is a function of κ and ω, the inversion is nonlinear, and the level of nonlinearity increases with increasing frequency. Additionally, if the simulated data at κ n =κ 0  and the measurement data are out of phase for more than a half period, then the gradient-based search direction will be incorrect, a phenomenon known as cycle-skipping, causing errors to propagate in subsequent models. 
         [0027]    In some embodiments of the present disclosure, the low wavenumber information buried in beat signal S B  may be extracted to carry out the FWI operation and reconstruct a model for the underground formation without suffering from cycle-skipping. Rather than inverting data from signal S 1  at frequency f 1 , or data from signal S 2  at frequency f 2 , or data from signal S 3  at frequencies f 1  and f 2  using the conventional FWI method, for example and without limitation, the inversion may be carried out utilizing the low wavenumber data extracted from beat signal S B , which is constructed from high frequency seismic data at two different frequencies through mathematical operations and signal processing techniques. 
         [0028]    In some embodiments, FWI may be carried out utilizing the amplitude data of beat signal S B . In such an amplitude-frequency differentiation (AFD) inversion, beat signal S B  may be defined as the envelope of S 3 , which may then be utilized as the input data of at least one iteration of the FWI operation. As an example,  FIGS. 2A and 2B  depict simulated data received by 313 receivers located evenly from 100 m to 4000 m from a signal source at 0 m with constant seismic velocity v.  FIG. 2A  depicts real parts of frequency domain data S 1  and S 2  overlaid as received by the receivers.  FIG. 2B  depicts S B , calculated by computing the absolute value of the amplitude or the envelope of the summation of frequency domain data from signals S 1  at f 1  and S 2  at f 2 .  FIG. 2B  also depicts S LF . S LF , shown only for comparison purposes, corresponds to a true low frequency signal captured at the frequency 
         [0000]    
       
         
           
             
               
                 
                   f 
                   1 
                 
                 - 
                 
                   f 
                   2 
                 
               
               2 
             
             . 
           
         
       
     
         [0000]    In this nonlimiting example, S 1  and S 2  are 15 and 12 Hz respectively, with S LF  thus being 1.5 Hz. As can be seen, although S B  is constructed from high frequency seismic data, it includes all the same low wavenumber information as the true low frequency S LF , also shown by: 
         [0000]    
       
         
           
             
               
                 S 
                 1 
               
               + 
               
                 S 
                 2 
               
             
             = 
             
               
                 2 
                 
                   | 
                   x 
                   | 
                 
               
                
               
                 cos 
                  
                 
                   ( 
                   
                     
                       
                         
                           k 
                           1 
                         
                         - 
                         
                           k 
                           2 
                         
                       
                       2 
                     
                      
                     x 
                   
                   ) 
                 
               
                
               
                 exp 
                  
                 
                   ( 
                   
                     
                       - 
                       i 
                     
                      
                     
                       
                         
                           k 
                           1 
                         
                         + 
                         
                           k 
                           2 
                         
                       
                       2 
                     
                      
                     x 
                   
                   ) 
                 
               
             
           
         
       
     
         [0000]    and thus: 
         [0000]    
       
         
           
             
               | 
               
                 
                   S 
                   1 
                 
                 + 
                 
                   S 
                   2 
                 
               
               | 
             
             = 
             
               2 
               | 
               
                 
                   cos 
                    
                   
                     ( 
                     
                       
                         
                           
                             k 
                             1 
                           
                           - 
                           
                             k 
                             2 
                           
                         
                         2 
                       
                        
                       x 
                     
                     ) 
                   
                 
                  
                 
                   / 
                 
                  
                 x 
               
               | 
             
           
         
       
     
         [0029]    By inverting |S 1 +S 2 | instead of S 1  or S 2  alone, for example and without limitation, during FWI processes, low wavenumber information contained in beat signal may contribute to the recovery of large scale structures while mitigating issues with cycle-skipping and local minima without using true low frequency seismic data. Because S B  includes all the same low wavenumber information as the true low frequency S LF , the FWI result obtained by inverting data from beat signal S B  is expected to be identical or close to that obtained by inverting data S LF . 
         [0030]    As an example,  FIG. 3A  depicts velocity anomalies as ideally detectable by the receivers.  FIG. 3B  depicts real parts of frequency domain data from signals S 1  and S 2  as computed.  FIG. 3C  depicts beat signal S B  (the envelope of S 1 +S 2 ) and real part of true low frequency data S LF , showing the effect of the anomalies on S B . 
         [0031]    In order to carry out FWI utilizing the AFD data, beat signal S B  may be used to determine an initial model after which single-frequency FWI (for example using only S 1 , S 2 , or another seismic signal at another frequency) may be utilized to generate subsequent, updated models. In other embodiments, beat signal S B  may be used to determine the initial model κ 0  as well as each subsequent model κ n . In some embodiments, for each iteration of the FWI, beat signal S B  may be generated by S 1  and S 2  at constant frequencies. In other embodiments, for each iteration of the FWI, beat signal S B  may be generated by S 1  and S 2  at different frequencies selected to generate a different beat signal S B . In some embodiments, for each iteration, beat signal S B  may be generated by a variation in one or both of S 1  and S 2 . 
         [0032]    At each iteration, the quality of the fit may be calculated according to the cost function for AFD, given by: 
         [0000]    
       
         
           
             
               C 
                
               
                 ( 
                 
                   m 
                   n 
                 
                 ) 
               
             
             = 
             
               
                 
                   1 
                   2 
                 
                 || 
               
               | 
               
                 
                   S 
                    
                   
                     ( 
                     
                       ω 
                       1 
                     
                     ) 
                   
                 
                 + 
                 
                   S 
                    
                   
                     ( 
                     
                       ω 
                       2 
                     
                     ) 
                   
                 
               
                
               
                 | 
                 2 
               
                
               
                 - 
                 
                   | 
                   
                     
                       M 
                        
                       
                         ( 
                         
                           ω 
                           1 
                         
                         ) 
                       
                     
                     + 
                     
                       M 
                        
                       
                         ( 
                         
                           ω 
                           2 
                         
                         ) 
                       
                     
                   
                    
                   
                     | 
                     2 
                   
                    
                   
                     
                       || 
                       2 
                     
                      
                     
                       
                         + 
                         
                           λ 
                           n 
                         
                       
                        
                       
                         
                           R 
                           n 
                         
                          
                         
                           ( 
                           
                             m 
                             n 
                           
                           ) 
                         
                       
                     
                   
                 
               
             
           
         
       
     
         [0000]    where S is the simulated data from the model, M is the measurement data, m is the unknown vector, a function of seismic velocity model v or compressibility model κ, to be inverted, λ is the regularization parameter, R is the regularization term, w is the angular frequency, and n denotes the iteration index. By minimizing C(m n ) or iterating until it is within an acceptable error tolerance, the fit of the generated model may be improved. As understood by one having ordinary skill in the art with the benefit of this disclosure, as the simulated data becomes closer to the measurement data, the value of the cost function decreases. 
         [0033]    In some embodiments, FWI may be carried out utilizing phase data of S 1 /S 2  to render low wavenumber information to inversion processes. In such a phase-frequency differentiation (PFD) inversion, beat signal data is defined as Φ(S 1 /S 2 ), i.e., phase of S 1 /S 2 , rather than envelope of S 1 +S 2  as in AFD. In some embodiments, dividing S 1  by S 2  may, for example and without limitation, reduce the effects of source estimation uncertainty, measurement error, or unnecessary cross talk between reflection energy. 
         [0034]    Utilizing the phase-frequency-differentiation data (i.e., beat signal data in PFD method) may mitigate or eliminate errors caused by, for example and without limitation, sensitivity to amplitude measurement error, unknown attenuation parameters Q, and uncertainty of mass density. 
         [0035]    As an example,  FIGS. 4A and 4B  depict simulated data received by 313 receivers located evenly from 100 m to 4000 m from a signal source at 0 m with constant seismic velocity v.  FIG. 4A  depicts the phases of frequency domain data from signals S 1  at frequency f 1  and S 2  at frequency f 2  overlaid as received by the receivers.  FIG. 4B  depicts the phase of S B  (i.e., Φ(S 1 /S 2 )).  FIG. 4B  also depicts the phase of S LF . The phase of S LF , shown only for comparison purposes, corresponds to a true low frequency signal captured at the frequency of f 1 -f 2 . As can be seen, the phase of S 1 /S 2  and S LF  are identical. 
         [0036]    As an example,  FIG. 5A  depicts velocity anomalies as ideally detectable by the receivers.  FIG. 5B  depicts the phases of S 1  and S 2  as received.  FIG. 5C  depicts S B  (i.e., Φ(S 1 /S 2 )) and the phase of S LF , showing the effect of the anomalies on the phase of S B . Because S B  and the phase of data S LF  are identical, the FWI result obtained by inverting data S B  is expected to be identical or close to that obtained by inverting the data S LF . 
         [0037]    In order to carry out FWI utilizing the PFD data, the beat signal S B  may be used to determine an initial model after which single-frequency FWI (for example using only S 1 , S 2 , or another seismic signal at another frequency) may be utilized to generate subsequent, updated models. In other embodiments, beat signal S B  may be used to determine the initial model KO as well as each subsequent model κ n . In some embodiments, for each iteration of the FWI, beat signal S B  may be generated by S 1  and S 2  at constant frequencies. In other embodiments, for each iteration of the FWI, beat signal S B  may be generated by S 1  and S 2  at different frequencies selected to generate a different beat signal S B . In some embodiments, for each iteration, beat signal S B  may be generated by a variation in one or both of S 1  and S 2 . 
         [0038]    At each iteration, the quality of the fit may be calculated according to the cost function for PFD, given by: 
         [0000]    
       
         
           
             
               C 
                
               
                 ( 
                 
                   m 
                   n 
                 
                 ) 
               
             
             = 
             
               
                 1 
                 2 
               
               || 
               
                 
                   Φ 
                    
                   
                     [ 
                     
                       
                         S 
                          
                         
                           ( 
                           
                             ω 
                             2 
                           
                           ) 
                         
                       
                        
                       
                         / 
                       
                        
                       
                         S 
                          
                         
                           ( 
                           
                             ω 
                             1 
                           
                           ) 
                         
                       
                     
                     ] 
                   
                 
                 - 
                 
                   Φ 
                    
                   
                     [ 
                     
                       
                         M 
                          
                         
                           ( 
                           
                             ω 
                             2 
                           
                           ) 
                         
                       
                        
                       
                         / 
                       
                        
                       
                         M 
                          
                         
                           ( 
                           
                             ω 
                             1 
                           
                           ) 
                         
                       
                     
                     ] 
                   
                 
               
                
               
                 || 
                 2 
               
                
               
                 
                   + 
                   
                     λ 
                     n 
                   
                 
                  
                 
                   
                     R 
                     n 
                   
                    
                   
                     ( 
                     
                       m 
                       n 
                     
                     ) 
                   
                 
               
             
           
         
       
     
         [0000]    where S is the simulated data from the model, M is the measurement data, m is the unknown vector to be inverted, λ is the regularization parameter, R is the regularization term, ω is the angular frequency, Φ is the phase extraction operator, and n denotes the iteration index. By minimizing C(m n ) or iterating until it is within an acceptable error tolerance, the fit of the generated model may be improved. The gradient vector of the cost function with respect to the unknown vector m, a function of seismic velocity model v or compressibility model κ, is calculated. The gradient vector alone, or combined with other information, may be utilized in the FWI inversion processes for velocity model update direction searching to minimize the cost function. 
         [0039]    As understood by one having ordinary skill in the art with the benefit of this disclosure, as the simulated data becomes closer to the measurement data, the value of the cost function decreases. 
         [0040]    As another example of PFD measurement data (i.e., beat signal data) to be inverted,  FIG. 6A  depicts a test model (here, the Marmousi model) with a single source and the associated multiple receivers depicted at the top thereof.  FIG. 6B  depicts the phase data of S 1  and S 2  as receive by the receivers.  FIG. 6C  depicts the beat signal data (i.e., Φ(S 1 /S 2 )) as a result of the PFD technique described hereinabove. As is clearly visible, the beat signal data (i.e., Φ(S 1 /S 2 )) includes far fewer phase-wraps, thus reducing the possibility of cycle-skipping. 
         [0041]    In order to further the understanding of the present disclosure, a merely exemplary and non-limiting PFD FWI operation is depicted in  FIGS. 7A-7D .  FIG. 7A  depicts the true seismic velocity model for the exemplary PFD FWI operation, here depicted as the Marmousi model.  FIG. 7B  depicts the initial model used. Data is extracted from a shot such that S 1  is at 5 Hz and S 2  is at 5.5 Hz and is used with the PFD inversion as discussed above, resulting in the model shown in  FIG. 7C . This model is used as the starting model for a FWI operation, which after several operations utilizing single frequency data results in the final model result depicted in  FIG. 7D . 
         [0042]    In order to further the understanding of this disclosure,  FIG. 8  depicts a flow chart of a FWI operation consistent with embodiments of the present disclosure.  FIG. 8  is intended as an example and is not intended as limiting the scope of the disclosure in any way. As previously discussed, seismic data is recorded ( 20 ). In some embodiments, the seismic data may be pre-processed ( 30 ). Pre-processing may include, for example and without limitation, filtering of the recorded data. 
         [0043]    The two frequencies f 1  and f 2  for beat signal construction may be determined ( 40 ) based on, for example and without limitation, knowledge of initial velocity model ( 10 ) and the quality and frequency range of the pre-processed data. The pre-processed time domain seismic data may be converted to frequency domain data ( 50 ) using signal processing techniques. The frequency domain data at frequency f 1  and f 2  may be extracted as M 1  and M 2  ( 60 ). The initial velocity model may be input to a seismic wave propagation engine to generate the simulated frequency domain data from signals S 1  with the frequency f 1  and S 2  with the frequency f 2  ( 70 ). The cost function and the data residual of AFD method or PFD method may then be calculated using the measurement data M 1  and M 2  and the simulated data for signals S 1  and S 2  ( 80 ) as discussed herein above. The gradient of the AFD or PFD cost function and the velocity update direction are calculated ( 100 ) to obtain an updated velocity model ( 120 ). The updated velocity model may be input into the seismic wave propagation engine to obtain a new set of simulated data S 1  and S 2  to be used for cost function and data residual evaluation ( 140 ). The updated velocity model and the data fitting are judged ( 160 ). If the data fitting is satisfactory and the update velocity model is accepted, the beat tone FWI process finishes ( 180 ). Otherwise, using the calculated gradient of the cost function, the velocity model may be updated ( 170 ) and used as the input to operations  100 ,  120 ,  130 ,  140 , and  160 . These operations may repeat until the data fitting is satisfactory, and the reconstructed velocity model is satisfactory. 
         [0044]    In some embodiments, additional operations may be included in the above described FWI operation. For example and without limitation, in some embodiments, the initial velocity model ( 10 ) may be obtained through well logs, velocity tomography procedures, or any other velocity analysis techniques. In some embodiments, in operation  40 , besides initial velocity model and pre-processed seismic data, other available information (e.g., seismic migration image, geology information) may be used to help determine f 1  and f 2 . 
         [0045]    In some embodiments, in operation  70 , the seismic wave propagation engine may be based on a time domain method or a frequency domain method. In some embodiments, it may be a finite-difference method, finite-element method, integral-equation method, or other numerical methods suitable for wave propagation simulations. 
         [0046]    In some embodiments, in operations  80  and  100 , extra regularization terms may be added to further stabilize the FWI processes. 
         [0047]    In some embodiments, in operation  120 , the updated velocity model may be further modified based on a priori information such as, for example and without limitation, water bottom profile or salt body geometry. 
         [0048]    The foregoing outlines features of several embodiments so that a person of ordinary skill in the art may better understand the aspects of the present disclosure. Such features may be replaced by any one of numerous equivalent alternatives, only some of which are disclosed herein. One of ordinary skill in the art should appreciate that they may readily use the present disclosure as a basis for designing or modifying other processes and structures for carrying out the same purposes and/or achieving the same advantages of the embodiments introduced herein. One of ordinary skill in the art should also realize that such equivalent constructions do not depart from the spirit and scope of the present disclosure and that they may make various changes, substitutions, and alterations herein without departing from the spirit and scope of the present disclosure.