Abstract:
A rigorous full waveform inversion of seismic data is presented that does not require any source information, therefore, eliminating potential errors involved in the physical and mathematical approximation of source input functions and their coupling to unrepeatable geological sites. A set of seismic trace data is collected in a an input data set that is first Fourier transformed in its entirety into the frequency domain. A normalized wavefield is obtained for each trace of the input data set in the frequency domain. Normalization is done with respect to the frequency response of a reference trace selected from the set of seismic trace data. The normalized wavefield is source independent, complex and dimensionless. The normalized wavefield is shown to be uniquely defined as the normalized impulse response, provided that a certain condition is met for the source. This property allows construction of the inversion algorithm disclosed herein, without any source or source coupling information. The algorithm minimizes the error between data normalized wavefield and the model normalized wavefield. The methodology is applicable to any 3-D seismic problem, and damping may be easily included in the process. A proof of principle of the invention is demonstrated using a simple 2-D scalar problem.

Description:
CROSS REFERENCE TO RELATED APPLICATIONS  
       [0001]    This application claims priority to U.S. Provisional Patent Application No. 60/456,168, filed Mar. 18, 2003 entitled “Method and article of manufacture for source-independent full waveform inversion of seismic data”, which is hereby incorporated by reference in its entirety. 
     
    
     STATEMENT REGARDING FEDERAL FUNDING  
       [0002] This invention was made with U.S. Government support under Contract Number DE-AC03-76SF00098 between the U.S. Department of Energy and The Regents of the University of California for the management and operation of the Lawrence Berkeley National Laboratory. The U.S. Government has certain rights in this invention. 
     
    
     
       REFERENCE TO A COMPUTER PROGRAM  
         [0003]    Not applicable.  
         BACKGROUND OF THE INVENTION  
         [0004]    1. Field of the Invention  
           [0005]    The present invention pertains generally to methods used in the inversion of data, more particularly to the inversion of seismic data, and still more particularly to the full wave seismic inversion of seismic data without using a source term.  
           [0006]    2. Description of the Relevant Art  
           [0007]    Common seismic industry practice is that velocity structure is estimated by analyzing the travel time of a set of seismic signals. In cross hole and surface-to-borehole applications, typical approaches involve ray tracing in which the ray may be straight or curved depending on the degree of resolution desired, and more recently the Fresnel volume approach. The travel time approach is fundamentally a high-frequency approximation with a maximum resolution on the order of a wavelength, or a fraction (5%) of the well separation in some practical cases. The usefulness of the result obtained from the ray tomography may be limited if the objective is to better understand the petrophysical and hydrological ‘properties’ of the soil and rock: an increasingly important subject in characterizing petroleum and geothermal reservoirs, with potential environmental applications.  
           [0008]    A better alternative to the travel time approach appears to be full waveform inversion. A number of recent studies on this subject suggest that full waveform inversion may provide improved resolution of velocity and density structures. Amplitude and phase information contained in the waveform are both sensitive to intrinsic energy loss of the propagating wave through dispersion, and to the petrophysical properties of the material through which the wave propagates. These sensitivities are utilized by full waveform analysis to form an attractive tool for investigating hydrological and petrophysical properties of a geophysical medium.  
           [0009]    Full waveform inversion has a major impediment to easy implementation. In all field applications, the effective source waveform and the coupling of the medium between the source and the receiver are not very well understood. The effective source waveform and its coupling to the medium under study are generally taken together as a lumped “source term”. The source/medium coupling problem can be resolved to some extent by using a good velocity approximation in geophysical applications; however, the measured signal cannot generally be accurately calibrated, rendering full waveform inversion a technical challenge.  
         BRIEF SUMMARY OF THE INVENTION  
         [0010]    This invention provides for an iterative method of using normalized waveform inversion to obtain a model parameter describing one or more physical properties of a medium, the method comprising the steps of: a) inputting one or more each of NS source and NG measurement spatial locations; b) measuring time domain data at each of the NG measurement locations resulting from an input waveform at one of the NS source locations propagating through a medium, i) for each of the NS source locations, ii) thereby forming a time domain measurement data set D ji   d (t), (j=1˜NG; i=1˜NS), where d denotes a data tensor; c) Fourier transforming the time domain measurement data set D ji   d (t) to create a measurement spectral data set D ji   d (ω) having frequency and amplitude information for each of the NG measurement locations; d) normalizing the measurement spectral data set D ji   d (ω) to create a normalized data wavefield T ji   d (ω); e) modeling the medium using an iterated model parameter m describing one or more physical properties of the medium, the NS source and the NG measurement spatial locations used as respective model input and model response spatial locations contained within the model of the medium, and initializing the iterated model parameter m with corresponding one or more known bulk properties of the medium being modeled, said modeling step comprising: i) creating a measurement model by: (1) applying a delta function source collocated with the i th  NS source, (2) modeling the response at the NG measurement locations, using a time domain modeling method, to create a synthetic medium response at the j th  receiver due to the i th  source, P ji   m (t), m denoting the model response, (3) repeating the applying and modeling steps at each of the NS source locations and NG measurement locations until the measurement model is full, and ii) Fourier transforming the model response P ji   m (t), to obtain a frequency domain synthetic response P ji   m (ω); iii) forming a normalized modeled wavefield using the frequency domain synthetic response T ji   m (ω))=P ji   m (ω)[P 1i   m (ω)] −1 ; and f) minimizing a weighted error, between the normalized data wavefield T ji   d (ω) and the normalized modeled wavefield T ji   m (ω) of the response of the medium, to a level below an error bound, i) said weighted error met by using the iterated model parameter m , known as the a minimized model parameter m .  
           [0011]    In the method described above the modeled normalized wavefield T ji   m (ω) is independent of the input waveform applied to any of the NS source locations, and may be used for a variety of applications, so long as they are capable of modeling the physical medium under test. Such modeling may utilize finite elements of various types, boundary elements, mixtures of both, and/or discrete finite differences. The results of the models at the iterative step concluding the minimization results in the minimizing model parameter m, which may be directly stored in a computer-readable medium, or immediately processed into a display of the minimizing model parameter m as a graphical representation of one or more of the physical properties of the medium measured, and optionally stored as a data set or data structure in a computer-readable medium. Although the thrust of this invention is highlighted herein to the field of seismic analysis, it is not restricted to such field.  
           [0012]    The methods described above may be further stored as a computer program in at least one computer-readable medium, to allow for execution on a computer.  
           [0013]    The method of normalized waveform inversion uses minimization of weighted error φ(m) that is calculated by a square of the weighing matrix W d  multiplied by the difference between the normalized data wavefield T ji   d (ω) and the normalized modeled wavefield T ji   m (ω), or φ(m)=∥W d (T ji   m (ω)−T ji   d (ω))∥ 2  summed over all frequencies, NS sources, and NG receivers. The weighing matrix W d  is an inverse of the standard deviation of the normalized data wavefield T ji   d (ω), or other appropriate weighting method indicative of the convergence of the minimization.  
           [0014]    The detailed normalization step described above is further comprised of the steps of: a) choosing one of the measurement spectral data set D ji   d (ω) at one of the measurement locations, such as D 1i   d (ω), i) said chosen measurement spectral data set D 1i   d (ω) known as a reference measurement; ii) dividing the measurement spectral data set D ji   d (ω) at a particular frequency in the measurement spectral data set D ji   d (ω), iii) by the reference measurement D 1i   d (ω) at the same frequency, b) to create a normalized data wavefield T ji   d (ω).  
           [0015]    The error bound of the minimizing step is met by first assigning an initial weighted error E 0 , then dividing each subsequent iterated weighted error E k  by the initial error E 0  to create a relative error  
           E   k       E   0       .                         
 
           [0016]    When the relative error  
         E   k       E   0                           
 
           [0017]    is below a pre-established error bound, preferably below 10 −5 , more preferably below 10 −6 , still more preferably below 10 −7 , and possibly below 10 −8 , the iterations of the minimization step conclude. Typically an error bound on the order of 10 −7  works particularly well. The model parameter m that achieved this error bound is then the solution to the minimization problem, and the wavefield inversion is complete.  
           [0018]    In another embodiment of this invention, modeling takes place not in the time domain, but rather in the frequency domain as a modal analysis. In this case, the delta function input, which has an infinite number of Fourier frequencies, is truncated at the highest frequency present in the model data. In this implementation, the frequency domain synthetic response P ji   m (ω) is directly obtained in the frequency domain. As with the previous methods, the frequency-based model, minimized model parameter m solution, and plots of the physical properties of the model parameter m solution may be stored on a computer-readable medium for further computations on a computer.  
           [0019]    In yet another embodiment, an article of manufacture comprises a computer-readable medium having the method of full waveform inversion described above placed upon it for subsequent execution on a computer.  
       
    
    
     BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWINGS  
       [0020]    The invention will be more fully understood by reference to the following drawings, which are for illustrative purposes only:  
         [0021]    [0021]FIG. 1A is a 2-D velocity model using a fault model in a background of 3000 m/s constant velocity  
         [0022]    [0022]FIG. 1B is the inversion result using impulse response of the model of FIG. 1A.  
         [0023]    [0023]FIG. 1C is the inversion result using normalized wavefield obtained for the model of FIG. 1A.  
         [0024]    [0024]FIG. 2 is a root mean square convergence showing the root mean square (rms) error and the associated Lagrange multipliers as a function of the number of iterations during the full waveform inversion. 
     
    
     DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT  
       [0025]    Defined Terms  
         [0026]    Computer means any device capable of performing the steps developed in this invention to result in an optimal waterflood injection, including but not limited to: a microprocessor, a digital state machine, a field programmable gate array (FGPA), a digital signal processor, a collocated integrated memory system with microprocessor and analog or digital output device, a distributed memory system with microprocessor and analog or digital output device connected with digital or analog signal protocols.  
         [0027]    Computer readable media means any source of organized information that may be processed by a computer to perform the steps developed in this invention to result in an optimal waterflood injection, including but not limited to: a magnetically readable storage system; optically readable storage media such as punch cards or printed matter readable by direct methods or methods of optical character recognition; other optical storage media such as a compact disc (CD), a digital versatile disc (DVD), a rewritable CD and/or DVD; electrically readable media such as programmable read only memories (PROMs), electrically erasable programmable read only memories (EEPROMs), field programmable gate arrays (FGPAs), flash random access memory (flash RAM); and remotely transmitted information by electromagnetic or optical methods.  
         [0028]    Physical property means a derived or observable, measurable characteristic of sample, including, but not limited to speed of sound in a specified direction, density, permeability, porosity, resistivity, permittivity, Young&#39;s modulus, tensor stress-strain properties, speed of sound propagation, shear modulus, amount of water present, and amount of oil present.  
         [0029]    Bulk property means a physical property of a medium spatially averaged over a volume of interest.  
         [0030]    Time domain measurement data set means a collection of data measured in the time domain such as a seismogram, possibly including spatial information of one or more source and/or receiver locations.  
         [0031]    Wavefield means modeled or measured data, or combination of such data, that obey differential equations describing wave motion in time domain or in frequency domain.  
         [0032]    Normalized data wavefield means measured data at any receiver point normalized by measured data at a reference point in the frequency domain.  
         [0033]    Normalized modeled wavefield means modeled (calculated) data at any receiver point normalized by modeled data at a reference point in the frequency domain.  
         [0034]    Introduction  
         [0035]    Source-independent full waveform inversion of seismic data can lead to high-resolution imaging of subsurface geophysical parameters that would otherwise have not been possible using other conventional methods. Detailed information of distribution of geophysical parameters are crucial in understanding the spatial distribution of gas and fluid, and their movement when repeated measurements are made over time. Potential applications of the proposed method include, but are not limited to exploration of oil, methane hydrate, and geothermal resources, characterization and monitoring of various petrological reservoirs, CO 2  sequestration, and environmental remediation processes.  
         [0036]    Normalized Wavefield Analysis  
         [0037]    Without loss of generality, let us assume a field survey involving NS transmitter positions and NG receiver positions of arbitrary configurations, but with known spatial locations. Additionally, the spatial coordinates of the boundary of the medium are also known. In acquiring data, the source and receiver deployment configuration is an important part of survey design for successfully achieving the desired survey objective. The full waveform inversion scheme described herein is general, since any arbitrary configuration is acceptable: surface or single borehole reflection, surface-to-borehole or borehole-to-surface conventional processing of time-lapse Vertical Seismic Profiling (VSP) data, cross-hole, or combinations of these.  
         [0038]    To develop the full 3-D normalized wavefield inversion method, a full tensor measurement is required with three linearly independent component measurements at each receiver spatial position for each one of three transmitting events at each source spatial position. The data tensor may be in the form of pressure, displacement, velocity or acceleration, and may be described in general as  
           D   ji   d ( t )= R   j ( t ){circumflex over (x)}P ji   d ( t ){circumflex over (x)} S   i ( t ),  j =1 ˜NG; i= 1 ˜NS )  (1)  
         [0039]    The superscript d indicates data and each constituent in this equation is a (3×3) tensor. Equation (1) states that the data D ji   d ( t ) recorded at the j th  receiver position is a multiple convolution of the actually transmitted source S i (t) that includes source system function and the radiation pattern caused by source-medium coupling, the impulse response P ji   d (t) of the medium at the j th  receiver due to the i th  source, and the receiver system function R j (t) including the medium-receiver coupling. In the following method of analysis, we will remove the R j (t) by assuming that the receiver (geophone) calibration is known and the effect of medium-receiver coupling to data is relatively minor compared to that of the source waveform. The source function S i (t) is generally not well understood, and is a subject of considerable ongoing research just by itself, in order to improve the quality of data sampling. The impulse response P ji   d (t) is the solution to the governing differential equation of the medium under test, with an impulse source in time at the i th  position. If we Fourier transform Equation (1), FT{(D,S,P)(t)}→(D,S,P)(ω), with the R j (t) term neglected, we get  
           D   ji   d (ω)= P   ji   d (ω) S   i (ω)  (2)  
         [0040]    The transformed frequency domain data contains all of the information the time series originally contained. In the time-series data, longitudinal and transverse waves of primary arrivals are followed by events of converted modes, multiple reflections, and all other scattering waves caused by heterogeneities in the sampling medium. The data at the j th  receiver position may then be described with Equation (2) as  
                   D     j                 i     d          (   ω   )       =         {       d   aji   d     ,     d   bji   d     ,     d   cji   d     ,     }          (   ω   )       =       [           d     1      aji     d           d     1      bji     d           d     1      cji     d               d     2      aji     d           d     2      bji     d           d     2      cji     d               d     3      aji     d           d     3      bji     d           d     3      cji     d           ]          (   ω   )           ,           (   3   )                               
 
         [0041]    Multiple subscripts ‘lkji’ are used to describe data elements. The last two subscripts ‘ji’ only indicate receiver (j=1, 2, . . . , NG) and source (i=1, 2, . . . , NS) positions, and the first two subscripts ‘lk’ are used to describe Cartesian field components l=(1, 2, 3). Elements of the tensor source function S i (ω) in equation (2) are unknown. The tensor source function is composed of three source vectors [s ai ,s bi , s ci ](ω) and each source vector may in turn be composed of three unknown Cartesian components; i.e., s ai =(s 1ai ,s 2ai , s 3ai ) T . So, the source tensor at the i-th position may be fully described by  
                 S   i          (   ω   )       =         {       s   ai     ,     s   bi     ,     s   ci       }          (   ω   )       =       [           s     1      ai             s     1      bi             s     1      ci                 s     2      ai             s     2      bi             s     2      ci                 s     3      ai             s     3      bi             s     3      ci             ]            (   ω   )     .                 (   4   )                               
 
         [0042]    The tensor impulse response of the medium that relates the diagonal impulse source at the i th  transmitter position to the measurements at the j th  receiver position may be written as  
                 P     j                 i     d          (   ω   )       =         {       p     1      ji     d     ,     p     2      ji     d     ,     p     3      ji     d     ,     }          (   ω   )       =       [           p     11      ji     d           p     12      ji     d           p     13      ji     d               p     21      ji     d           p     22      ji     d           p     23      ji     d               p     31      ji     d           p     32      ji     d           p     33      ji     d           ]            (   ω   )     .                 (   5   )                               
 
         [0043]    To define the normalized wavefield (alternatively known as the transfer function of the medium under test), we first select the reference point (which may be any of the measurement locations), say j=1 for convenience. The tensor normalized data wavefield (or data transfer function) is now defined by the data at j=1˜NG, normalized by that of the reference point:  
                       T   ji   d          (   w   )       =           D     j                 i     d          (   w   )            [       D     1                 i     d          (   w   )       ]         -   1                   =         P     j                 i     d          (   w   )                      S   i          (   w   )            [       S   i          (   w   )       ]         -   1            [       P     1      i     d          (   w   )       ]         -   1                       =             P     j                 i     d          (   w   )            [       P     1      i     d          (   w   )       ]         -   1       =       [           t     11      ji                  d             t     12      ji                  d             t     13      ji                  d                 t     21      ji                  d             t     22      ji                  d             t     23      ji                  d                 t     31      ji                  d             t     32      ji                  d             t     33      ji                  d             ]          (   w   )           ,                 (   6   )                               
 
         [0044]    where the convention of the subscripts for the transfer function is the same as that used for impulse response. In Equation (6) the source term S i (ω) divides out to unity, so the transfer function is defined as the normalized impulse response of the medium; hence the uniqueness of the transfer function. The necessary condition for the source term to divide out is that the determinant of the source matrix is non-zero. In other words, the source events need to be linearly independent. This condition can be met, in principle, even if source events consist of explosions, so long as their Cartesian constituents are linearly independent.  
         [0045]    In this procedure we assume that for each source, i=1˜NS, NG measurements are made simultaneously. If for some reason all NG data cannot be taken simultaneously, a simultaneous data set may be simulated if neighboring subsets share at least one overlapping receiver position. Since we do not know the medium, the impulse response itself is not known either at this point.  
         [0046]    Full Waveform Inversion Using Normalized Wavefield Analysis  
         [0047]    To obtain the numerical solution for the media impulse response for scalar and elastic wave equations one needs to spatially discretize the constitutive parameters, and apply finite difference, finite element, or integral equation (boundary element) techniques to solve the governing differential equation of the discretized system. Using the finite difference or finite element method on a mass/spring/damper model of the medium under test, the assembled system of equations, including the damping, takes a general form:  
           M{umlaut over (p)} ( t )+ C{dot over (p)} ( t )+ Kp ( t )= g δ( t )  (7)  
         [0048]    where the field vector p (t) is the discretized wavefield, M is the mass matrix, K is the stiffness matrix, and C is the damping matrix. If there are a total of N unknowns in the discretization, then all M, C, and K matrices are N×N square, the field vector p is N×1, and the load vector g is also N×1 whose entries are all zero except for ‘1’ values at the source locations. Boundary conditions are implicitly included by the known spatial coordinates, and any level of continuity of the derivatives of the boundary using well-known finite element, boundary element, or finite difference methods, depending on the model solution method utilized.  
         [0049]    The reduced system of equations may be solved in the time domain, typically by using coupled first-order differential equations on a staggered grid, or in the frequency domain after Fourier transforming Equation (7), FT{p(t)}→p(ω), into  
         −ω 2   Mp (ω)+ iωCp (ω)+ Kp (ω)=g  (8)  
         [0050]    Next, we will show that the normalized wavefield defined by Equation (6) is all that is required for full waveform inversion. In the inversion the objective error functional typically consists of a measure of the error between the measured data and the resulting model fit. In this case, the normalized wavefield can be used to evaluate the error term. For a given model one can generate synthetic data using Equation (8), and then obtain the model normalized wavefield similar to the data normalized wavefield described by Equation (6).  
         [0051]    Formally, the synthetic impulse response at the j th  receiver due to the i th  source may be obtained and designated as P ji   m (ω),  
                 P     j                 i     m          (   ω   )       =       [           p     11      ji     m           p     12      ji     m           p     13      ji     m               p     21      ji     m           p     22      ji     m           p     23      ji     m               p     31      ji     m           p     32      ji     m           p     33      ji     m           ]            (   ω   )     .               (   9   )                               
 
         [0052]    where the superscript m indicates model. Accordingly, the model normalized wavefield (or transfer function) is obtained from the numerical solution for the given velocity model  
                 T     j                 i     m          (   ω   )       =               P   ji   m          (   ω   )            [     P     1                 i     m     ]         -   1            (   ω   )       =       [           t     11      ji     m           t     12      ji                  m             t     13      ji                  m                 t     21      ji                  m             t     22      ji                  m             t     23      ji                  m                 t     31      ji                  m             t     32      ji                  m             t     33      ji                  m             ]            (   ω   )     .                 (   10   )                               
 
         [0053]    The inversion procedure begins with a measure of the error between the normalized wavefield (or transfer function) model T m  and the data T d , with the subscripts to the normalized wavefield dropped  
         φ( m )=∥ W   d ( T   m   −T   d )  (11)  
         [0054]    The misfit error between data and model transfer functions at the reference point is always zero (t lk1i   m −t lk1i   d =0). In setting up the error functional of Equation (11), both real and imaginary parts are separated, so the actual number of data parts used for the inversion is NEQ=2×NFREQ×NT×(NG−1)×9, and the computation is done in real, or extended real arithmetic. The variable NFREQ is the total number of frequencies. The matrix W d  is an NEQ×NEQ weighing matrix that is the inverse of the standard deviation of the normalized data wavefield T d .  
         [0055]    For inversion we consider Gauss-Newton method by first expanding the objective functional, Equation (11), into a Taylor series  
         φ( m+δm )=φ( m )+γ m   T   δm+ 0.5 δm   T   H   m   δm+ 0{(δ m )}  (12)  
         [0056]    Here, δm is a perturbation to the minimizing model parameter m, γ m  is an M×1 column matrix consisting of elements  
       {         ∂   φ       ∂     m   p         ,     p   =     1   ∼   M         }                         
 
         [0057]    with M being the total number of parameters to be determined, and is compactly written as γ m =2J T W d   T W d (T m −T d ) and H m  is an M×M square matrix (a weighted and constrained Hessian) consisting of elements  
       {           ∂   2        φ         ∂     m   q            ∂     m   p           ,   p   ,     q   =     1   ∼   M         }                         
 
         [0058]    written as  
         H   m     =       2                   J   T          W   d   T          W   d        J     +       O        (       ∂   J       ∂   m       )       .                             
 
         [0059]    The last term above,  
         O        (       ∂   J       ∂   m       )       ,                         
 
         [0060]    is literally interpreted as the change in the partial derivative of data (transfer function in this case) due to a change in the velocity. This term is small if either the residuals are small, or the forward differential equation is quasi linear. Furthermore, it is usually difficult to compute, so is generally dropped. For each frequency and source the sensitivity function (Jacobian) J is a {2×(NG−1)×9}×M rectangular matrix: {2×(NG−1)×9} because real and imaginary parts have been separated for each of 9 tensor elements. For example, for the i th  source and a fixed frequency, the entries to the Jacobian corresponding to the j th  receiver and the p th  parameter may be obtained as  
                   (           J     p   ,   l   ,   k   ,     (       2   *   j     -   1     )     ,   i                 J     p   ,   l   ,   k   ,     (     2   *   j     )     ,   i             )     =       (         real           imaginary         )            ∂               ∂     m   p              t     l                 k                 j                 i     m         ;   l     ,     k   =   1     ,   2   ,   3   ,     
        with           (   13   )                   ∂     t     l                 k                 j                 i     m         ∂     m   q         =           ∂               ∂     m   q                  t   _       l                 k                 j                 i     m            p     1                 i     m              =       1          p     1                 i     m                 (         ∂       t   _       l                 k                 j                 i     m         ∂     m   q         -           t   _       l                 k                 j                 i     m            p     1                 i     m                     ∂          p     1                 i     m                           ∂     m   q             )                 (   14   )                               
 
         [0061]    Here, |p 1i   m | is the determinant of the tensor impulse response at the reference receiver point and {{overscore (t)} lk1i   m , l, k=1, 2, 3} are the result of multiplication of p ji   m (ω) and the adjoint of p 1i   m (ω) (see Appendix A). The partial derivatives of the normalized wavefield can be explicitly written as a function of the partial derivatives of the impulse responses (see Appendix B). Partial derivatives of the impulse responses may be obtained using the forward modeling results using equation (8). The full waveform inversion of seismic data based on equations (11) through (14) will not require the knowledge of the actual source waveform, which is an advantageous feature of this method.  
         [0062]    The error functional that will be minimized consists of the error term, Equation (12), and a constraint that will have a smoothing effect on the variation of the model in the updating process. Specifically, it may be written as  
         Φ( m+δm )=φ( m+δm )+λ∥ W   m   δm∥   2   (15)  
         [0063]    where λ is the Lagrange multiplier that controls relative importance of data errors and the behavior of the parameter variation, and W m , which is an M×M controlling matrix. When the W m  matrix is diagonal it has an effect of keeping the parameter from changing from the current one. On the other hand, if the matrix represents a gradient operator its effect is to spatially smooth out the changes. Minimization of the functional of Equation (15) with respect to perturbation in model parameter results in a system of normal equations  
         ( J   T   W   d   T   W   d   J+λW   m   T   W   m )δ m=−J   T   W   d   T   W   d ( T   m   −T   d )  (16)  
         [0064]    from which the model parameter at the (k+1)-th iteration is updated to m k+1 =m k +δm.  
         [0065]    The iteration process stops when the change in the model parameter δm is below a preset tolerance, typically given in terms of root mean square (rms).  
         [0066]    Construction of the Source Function  
         [0067]    Once the inversion is successful, the functional description of the actual source waveform can be established by deconvolution. The deconvolution can be achieved by dividing the data with the impulse response of the synthetic model that has been obtained through the inversion process. From Equation (2), the deconvolution tensor may be estimated to be  
                 C     j                 i            (   ω   )       =               D   ji   d          (   ω   )            [     P   ji   m     ]         -   1            (   ω   )       =       [           c     1      aji             c     1      bji             c     1      cji                 c     2      aji             c     2      bji             c     2      cji                 c     3      aji             c     3      bji             c     3      cji             ]            (   ω   )     .                 (   17   )                               
 
         [0068]    Notice that the elements of the impulse response of the model P ji   m (ω) are independent of the source events (a, b, c), but the elements of the deconvolution tensor C ji (ω) are related to the source events. If the data and the inversion process were exact, the deconvolution function C ji (ω) is actually the frequency spectrum of the source, and it must therefore be identical for all receiver position. Since neither the data nor the inversion process are exact, the source function may be obtained by taking a statistical mean of the deconvolution function at different receiver positions. For example the first element of the source function in Equation (4) may be estimated using the root mean square approach  
                 s     1                 a                 i            (   ω   )       =           1     N                 G              ∑     j   =   1       N                 G                         {       c     1      aji            (   ω   )       }     2           .             (   18   )                               
 
         [0069]    Since the function C ji (ω) is itself complex, the averaging process must be done in amplitude and phase, or real and imaginary parts, separately. Inverse Fourier transforming, FT −1 {S i (ω)}→S i (t), one obtains the effective source function in time:  
                 S   i          (   t   )       =         {       s   ai     ,     s   bi     ,     s   ci       }          (   t   )       =       [           s     1      ai             s     1      bi             s     1      ci                 s     2      ai             s     2      bi             s     2      ci                 s     3      ai             s     3      bi             s     3      ci             ]            (   t   )     .                 (   19   )                               
 
         [0070]    Note that each one of three events (a,b,c), for example s ai (t)={s 1ai ,s 2ai ,s 3ai } T (t), will produce a complete description of the effective source at the transmitter location i as a function of time. The function S i (t) describes the source at the i th  transmitter position that would most likely have produced data D ji   d (t), j=1˜NG, shown in Equation (1). The reconstructed source function is the effective source at the time of the event, and will also contain the spatial and temporal radiation patterns.  
         [0071]    Full Waveform Inversion of 2-D Acoustic Velocity  
         [0072]    The inversion method described here has been tested using a simple 2-D acoustic model. Let us consider the impulse response governed by a 2-D acoustic wave equation in the frequency domain  
                     ∇   2          p        (     x   ,     x   s     ,   ω     )         +         ω   2         v   2          (   x   )              p        (     x   ,     x   s     ,   ω     )         +     δ        (     x   -     x   s       )         =   0     ,           (   20   )                               
 
         [0073]    where the impulse response p is the scalar pressure wavefield, v the velocity, and (x,x s ) are the field and source positions in 2-D. The source is an impulse point source expressed as a 2-D spatial delta function δ(x−x s ) located at x s , and is also a delta function at t=0 in the time-domain formulation. A finite-element modeling (FEM) scheme is used to generate the synthetic impulse response. The model parameter is the acoustic velocity in each of the rectangular elements used for the FEM solution. Following the procedure described in the previous section, the scalar synthetic normalized wavefield is obtained from the numerical solution for the given velocity model  
                   t     j                 i                  m            (   ω   )       =         p     j                 i     m          (   ω   )           p     1                 i     m          (   ω   )           ,     (       j   =     1   ∼     N                 G         ;     i   =     1   ∼     N                 S           )     ,           (   21   )                               
 
         [0074]    The inversion procedure starts with an objective functional, Equation (15), reduced to handle scalar problem. The number of equations is NEQ=2×NFREQ×NT×(NG−1), with all computations done in real arithmetic. The related sensitivity functions are  
                 (           J     p   ,     (       2   *   j     -   1     )     ,   i                 J     p   ,     (     2   *   j     )     ,   i             )     =       (           real                        imaginary         )            ∂               ∂     m   p              t                  j                 i                    m           ,   with                     ∂     t                  j                 i                    m           ∂     m   p         =           ∂               ∂     m   p                p     j                 i     m       p     1                 i     m         =       1     p     1                 i     m            (         ∂     p     j                 i     m         ∂     m   p         -         p     j                 i     m       p     1                 i     m              ∂     p     1                 i     m         ∂     m   p             )           ;                 i   =     1   ∼     N                 S         ;     j   =     2   ∼   NG       ;     p   =     1   ∼     M   .                                     
 
         [0075]    [0075]FIGS. 1A shows a model  100  used for the test. The model is a broken dipping fault in a background of material  105  having a 3000 m/s constant velocity. The fault consists of a 6 m thick low velocity (2500 m/s) layer  130  overlain by another 6 m thick high velocity (3500 m/s) layer  140  starting an x position of 30-60 meters, and discontinuously dipping respectively toward the x=−55 meters position with low velocity section  135  and high velocity overlay  145 . A cross-hole configuration is used for the exercise, with the borehole at x=−45 m for the transmitter (Tx) borehole  110  and the other at x=45 m for the receiver (Rx) borehole  120 . A total of 21 line sources starting at elevation z=0 are used with an equal vertical separation of 9 m, and the same number of receivers and separations for the receivers. Each of the source pressure wavefields computed at 21 receiver positions have been normalized by the first pressure wavefield, resulting in 21 transfer functions. The number of frequencies used is 10; starting from 10 Hz to 100 Hz, linearly separated by 10 Hz. The shortest wavelength in the background is 30 m, so the maximum resolution is expected to be on the order of 7-8 m if we take into account the wavelength as a measure of resolution.  
         [0076]    Using the numerical solution of the model as synthetic data, the inversion was initiated with an initial model of 2850 m/s uniform whole space. A grid of uniform cell size, 3 m by 3 m, was used throughout. The inversion domain was 120 m by 180 m, containing a total of 2400 velocity parameters.  
         [0077]    The size of the matrix from Equation (15) is relatively modest for the test model, so it was solved by using QR decomposition with successive Householder transformations. The Lagrange multiplier λ is automatically selected in the inversion process. It starts with executing a given number, say nl, of inversions using nl different multipliers that are spaced appropriately. The same Jacobian is used at this step. As a result nl updated parameter sets are produced, followed by nl forward model calculations resulting in nl error terms. Among these, we choose the model and parameter λ giving the lowest error term result.  
         [0078]    To demonstrate the validity of the inversion method disclosed here, we first carried out conventional inversion using impulse response, with the result is shown in FIG. 1B, which uses the same numbering scheme of FIG. 1A. Additional features of a slightly slower velocity  150  and a slightly higher velocity  160  are pointed out as artifacts of the limited data set reconstruction. Using standard conventions, we assume that the source waveform is known, so that the data can be reduced to an impulse response, with an inversion carried out by using the impulse response.  
         [0079]    It should be pointed out that this numerical simulation has an ideal impulse response, with ideal coupling of the impulse to the medium. This situation is far from common in the field, however, does provide a comparison with the transfer function methods of this invention, which does not depend on source-to-medium coupling.  
         [0080]    We additionally obtained the velocity structure by using the transfer function approach, with the result shown in FIG. 1C, which also uses the same numbering scheme of FIG. 1A. Additional features of FIG. 1B showing a slightly slower velocity  150  is no longer readily apparent in the transfer function approach of FIG. 1C. However, the slightly higher velocity  160  of FIG. 1B is still somewhat apparent, but reduced in FIG. 1C. Again, such artifacts are due to the limited data set reconstruction. In this exercise we used nl=3 in each iteration to select parameter update and Lagrange multiplier. After 5 iterations, two of the results appear almost identical. Slight differences may have been due to the fact that the transfer function approach has one less data sample than the impulse response approach because one data sample was used to normalize the others. The behaviors of the rms errors are also similar as discussed below.  
         [0081]    [0081]FIG. 2 compares in chart  200  the rms errors of both the impulse response  210  to the transfer function method  220 . Note that the Lagrange multiplier has also changed as with iteration. Coordinate points  211 ,  212 ,  213 ,  214 , and  215  are the first five iterations of the impulse response  210 . Similarly, coordinate points  221 ,  222 ,  223 ,  224  and  225  are the first five iterations of the transfer method  220  described herein. It is seen that the transfer method generally yields a better result, and is quite similar in response to the ideal impulse response  210  method. The quality of inversion could be improved by using higher frequency data and by using a higher density deployment of transmitters and receivers, but the method has been successfully demonstrated.  
         [0082]    Conclusions  
         [0083]    An innovative, rigorous full waveform inversion method has been disclosed, and the validity of the method successfully demonstrated using a simple 2-D synthetic model. An important advantage of the invented method is that full waveform inversion of seismic data may be accomplished without the source information by using the properties of the normalized wavefield describing the system (typically a geological sample) under test. Under the proper combination of sources and receivers, the normalized wavefield is shown to be uniquely determined in terms of the normalized impulse response. The methodology is readily extended to include general 3-D problems. As an important byproduct, it has also been shown that the source function can be reconstructed once the full waveform inversion is completed. The reconstructed source function describes the effective source, not the source system output, at the time of event, including the spatial and temporal radiation patterns.  
         [0084]    All publications, patents, and patent applications mentioned in this specification are herein incorporated by reference to the same extent as if each individual publication or patent application were each specifically and individually indicated to be incorporated by reference.  
         [0085]    The description given here, and best modes of operation of the invention, are not intended to limit the scope of the invention. Many modifications, alternative constructions, and equivalents may be employed without departing from the scope and spirit of the invention.  
         [0086]    Appendix A—The Normalized Wavefield  
         [0087]    The relationship between the normalized wavefield and the impulse response is derived. We first rewrite the normalized wavefield defined in equation (6), or equation (10). After dropping the superscript d, the source position index i, and the angular frequency ω, we write  
                 T   j     =           P   j          [     P   1     ]         -   1       =     [           t     11      j             t     12      j             t     13      j                 t     21      j             t     22      j             t     23      j                 t     31      j             t     32      j             t     33      j             ]         ,           (     A        -        1     )                               
 
         [0088]    along with the impulse response, equation (5), as  
               P   j     =       [           p     11      j             p     12      j             p     13      j                 p     21      j             p     22      j             p     23      j                 p     31      j             p     32      j             p     33      j             ]     .             (     A        -        2     )                               
 
         [0089]    The matrix operation in equation (A-1) may be expanded to  
                     T   j     =           P   j          [     P   1     ]         -   1       =              [           t     11      j             t     12      j             t     13      j                 t     21      j             t     22      j             t     23      j                 t     31      j             t     32      j             t     33      j             ]     =         [           p     11      j             p     12      j             p     13      j                 p     21      j             p     22      j             p     23      j                 p     31      j             p     32      j             p     33      j             ]          [           p   111           p   121           p   131               p   211           p   221           p   231               p   311           p   321           p   331           ]         -   1                       =              1          P   1                 [           p     11      j             p     12      j             p     13      j                 p     21      j             p     22      j             p     23      j                 p     31      j             p     32      j             p     33      j             ]                            [               p   221          p   331       -       p   321          p   231                   p   321          p   131       -       p   121          p   331                   p   121          p   231       -       p   221          p   131                       p   311          p   231       -       p   211          p   331                   p   111          p   331       -       p   311          p   131                   p   211          p   131       -       p   111          p   231                       p   211          p   321       -       p   311          p   221                   p   311          p   121       -       p   111          p   321                   p   111          p   221       -       p   211          p   121               ]                   =              1          P   1                 [             t   _       11      j               t   _       12      j               t   _       13      j                   t   _       21      j               t   _       22      j               t   _       23      j                   t   _       31      j               t   _       32      j               t   _       33      j             ]         ,                 (     A        -        3     )                               
 
         [0090]    where the determinant of the reference impulse response is given by  
         | P|=p   111   p   221   p   331   +p   121   p   231   p   311   +p   131   p   321   p   211   −p   311   p   221   p   131   −p   321   p   231   p   111   −p   331   p   121   p   211   (A-4)  
         [0091]    and elements of the reduced tensor normalized wavefield are  
           {overscore (t)}   11j   =p   11j (p 221   p   331   −p   321   p   231 )+ p   12j ( p   311   p   231   −p   211   p   331 )+p 13j ( p   211   p   321   −p   311   p   221 ) 
           {overscore (t)}   12j   =p   11j (p 321   p   131   −p   121   p   331 )+ p   12j ( p   111   p   231   −p   311   p   131 )+p 13j ( p   311   p   121   −p   111   p   321 ) 
           {overscore (t)}   13j   =p   11j (p 121   p   231   −p   221   p   131 )+ p   12j ( p   211   p   131   −p   111   p   231 )+p 13j ( p   111   p   221   −p   211   p   121 )  
           {overscore (t)}   21j   =p   21j (p 221   p   331   −p   321   p   231 )+ p   22j ( p   311   p   231   −p   211   p   331 )+p 23j ( p   211   p   321   −p   311   p   221 ) 
           {overscore (t)}   22j   =p   21j (p 321   p   131   −p   121   p   331 )+ p   22j ( p   111   p   331   −p   311   p   131 )+p 23j ( p   311   p   121   −p   111   p   321 ) 
           {overscore (t)}   23j   =p   21j (p 121   p   231   −p   221   p   131 )+ p   22j ( p   211   p   131   −p   111   p   231 )+p 23j ( p   111   p   221   −p   111   p   121 )  
           {overscore (t)}   31j   =p   31j (p 221   p   331   −p   321   p   231 )+ p   32j ( p   311   p   231   −p   211   p   331 )+p 33j ( p   211   p   321   −p   311   p   221 ) 
           {overscore (t)}   31j   =p   31j (p 321   p   131   −p   121   p   331 )+ p   32j ( p   111   p   331   −p   311   p   131 )+p 33j ( p   311   p   121   −p   111   p   321 ) 
           {overscore (t)}   31j   =p   31j (p 121   p   231   −p   221   p   131 )+ p   32j ( p   211   p   131   −p   111   p   231 )+p 33j ( p   111   p   221   −p   211   p   121 )  (A-5)  
         [0092]    Appendix B—The Partial Derivatives of the Normalized Wavefield  
         [0093]    In Appendix A, each element of the normalized wavefield was shown to be a rational function of combination of elementary impulse responses. Let us rewrite equation (14) for the partial derivatives of the elementary normalized wavefield with respect to the q th  parameter m q . After dropping the superscript m and the source position index i, we get  
                   ∂     t   lkj         ∂     m   q         =         ∂     ∂     m   q                  t   _     lkj            p   1              =       1          p   1                 (         ∂       t   _     lkj         ∂     m   q         -           t   _     lkj            p   1                   ∂          p   1              ∂     m   q             )           ,           (     B        -        1     )                               
 
         [0094]    where the determinant is given by equation (A-4). To illustrate the process of obtaining partial derivatives of normalized wavefields, we will choose an element, say {overscore (t)} 31j , which is given in equation (A-5). Partial derivatives of |P 1 | and {overscore (t)} 31j  with respect to the parameter m q  are:  
                         ∂          p   1              ∂     m   q         =                  ∂     p   111         ∂     m   q              p   221          p   331       +       p   111            ∂     p   221         ∂     m   q              p   331       +       p   111          p   221            ∂     p   331         ∂     m   q           +                                ∂     p   121         ∂     m   q              p   231          p   311       +       p   121            ∂     p   231         ∂     m   q              p   311       +       p   121          p   231            ∂     p   311         ∂     m   q           +                                ∂     p   131         ∂     m   q              p   321          p   211       +       p   131            ∂     p   321         ∂     m   q              p   211       +       p   131          p   321            ∂     p   211         ∂     m   q           -                                ∂     p   311         ∂     m   q              p   221          p   131       -       p   311            ∂     p   221         ∂     m   q              p   131       -       p   311          p   221            ∂     p   131         ∂     m   q           -                                ∂     p   321         ∂     m   q              p   231          p   111       -       p   321            ∂     p   231         ∂     m   q              p   111       -       p   321          p   231            ∂     p   111         ∂     m   q           -                                ∂     p   331         ∂     m   q              p   121          p   211       -       p   331            ∂     p   121         ∂     m   q              p   211       -       p   331          p   121            ∂     p   211         ∂     m   q                          
        and           (     B        -        2     )                         ∂       t   _       31      j           ∂     m   q         =                  ∂     p     31      j           ∂     m   q              (         p   221          p   331       -       p   321          p   231         )       +                              p     31      j            (           ∂     p   221         ∂     m   q              p   331       +       p   221            ∂     p   331         ∂     m   q           -         ∂     p   321         ∂     m   q              p   231       -       p   321            ∂     p   231         ∂     m   q             )       +                                ∂     p     32      j           ∂     m   q              (         p   311          p   231       -       p   211          p   331         )       +                              p     32      j            (           ∂     p   311         ∂     m   q              p   231       +       p   311            ∂     p   231         ∂     m   q           -         ∂     p   211         ∂     m   q              p   331       -       p   211            ∂     p   331         ∂     m   q             )       +                                ∂     p     33      j           ∂     m   q              (         p   211          p   321       -       p   311          p   221         )       +                              p     33      j            (           ∂     p   211         ∂     m   q              p   321       +       p   211            ∂     p   321         ∂     m   q           -         ∂     p   311         ∂     m   q              p   221       -       p   311            ∂     p   221         ∂     m   q             )       .                   (     B        -        3     )                               
 
         [0095]    Substituting (B-2) and (B-3) into (B-1), we get the partial derivative of t 31j  with respect to the parameter m q . So, partial derivatives of the normalized wavefield may be evaluated in terms of a combination of partial derivatives of the impulse responses:  
                     ∂     t   lkj         ∂     m   q         =     f        (       ∂     p   lkj         ∂     m   q         )         ;   l     ,     k   =   1     ,   2   ,     3   ;     j   =   1       ,   2   ,   3   ,   …              ,   NG   ,           (     B        -        4     )                               
 
         [0096]    where NG is the number of receiver measurement positions.