Abstract:
Source-receiver compression is used to help design surveys and mitigate the computational costs of data set inversion. The source-receiver compression is based on data redundancy and sensitivity. More particularly, a compressed source array is produced for minimum redundancy and maximum sensitivity to reservoir model parameters. The synthesized transmitter array has a reduced number of sources, thereby reducing the number of forward model simulations needed to carry out the inversion. Furthermore, the data collected at the receivers employed in the survey can be compressed. This has the implication of reducing the computational cost of constructing the Jacobian matrix and inverting the corresponding Hessian matrix.

Description:
BACKGROUND OF THE INVENTION 
     The present invention is generally related to surveying subterranean formations to facilitate recovery of natural resources, and more particularly to designing surveys and reducing the computational overhead associated with processing survey data. 
     Developments in numerical computation techniques have motivated research on Full-Waveform Inversion (FWI) approaches for geophysical applications. For example, R. G. Pratt, C. Shin, and G. J. Hicks, Gauss-newton and full newton methods in frequency-space seismic waveform inversion, Geophysical Journal International, 13, 341-362 (1998); C. Shin, K. Yoon, K. J. Marfurt, K. Park, D. Yang, H. Y. Lim, S. H. Chung, and S. Shin, Efficient calculation of a partial derivative wavefield using reciprocity for seismic imaging and inversion, Geophysics, 66, 1856-1863 (2001); A. Abubakar, P. van den Berg, and J. T. Fokkema, Towards non-linear inversion for characterization of timelapse phenomena through numerical modelling, Geophysical Prospecting, 51, 285-293 (2003); S. Operto, J. Virieux, P. Amestoy, J.-Y. L&#39;Excellent, L. Giraud, and H. Ben-Hadj-Ali, 3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver, A feasibility study: Geophysics, 72, SM195-SM211 (2007); D. Vigh, and E. W. Starr, 3D prestack plane-wave fullwaveform inversion, Geophysics, 73, VE135-VE144 (2008); W. Hu, A. Abubakar, and T. M. Habashy, Preconditioned non-linear conjugate gradient method for seismic full-waveform inversion, in Expanded Abstracts, U018, 71st EAGE Conference &amp; Exhibition (June 2009); W. Hu, A. Abubakar, and T. M. Habashy, Simultaneous multifrequency inversion of fullwaveform seismic data, Geophysics, 74, R1-R14 (2009); Abubakar, W. Hu, T. M. Habashy, and P. M. van den Berg, Application of the finite-difference contrast-source inversion algorithm to seismic full-waveform data, Geophysics, 74, WCC163-WCC174 (2009). FWI seismic inversion generally deals with processing a large size data set, which tends to be a processing bottleneck. A large number of sources in the survey contributes to a large computational cost in running the forward simulator a number of times corresponding to the number of sources. A large number of receivers contributes to the computational cost of constructing the sensitivity (Jacobian) matrix as well as in inverting the Hessian matrix in a gradient-type inversion approach. Recently a simultaneous-source encoded FWI approach has been proposed to reduce the number of sources used in the inversion, see H. Ben-Hadj-Ali, S. Operto, and J. Virieux, Efficient 3d frequency-domain full waveform inversion (fwi) with phase encoding, in Expanded Abstracts, P004, 71st EAGE Conference &amp; Exhibition (2009); F. J. Herrmann, Y. Erlangga, and T. T. Y. Lin, Compressive sensing applied to fullwave form inversion, in Expanded Abstracts, 5016, 71st EAGE Conference &amp; Exhibition (2009); J. Krebs, J. Anderson, D. Hinkley, R. Neelamani, S. Lee, A. Baumstein, and M. D. Lacasse, Fast full-wavefield seismic inversion using encoded sources Geophysics (2009). In this approach a large number of physical sources are converted into one simultaneous source or several simultaneous sources by summing up the individual physical sources using a phase encoding technique (see S. A. Morton, and C. C. Ober, Faster shot-record depth migrations using phase encoding, in Expanded Abstracts, 1131-1134, SEG Annual Meeting (1998)). This phase encoding technique has also been applied for the prestack migration as in L. A. Romero, D. C. Ghiglia, C. C. Ober, and S. A. Morton, Phase encoding of shot records in prestack migration, Geophysics (2000). It has been pointed out that this approach is more sensitive to noise than the standard FWI method. J. Krebs, J. Anderson, D. Hinkley, R. Neelamani, S. Lee, A. Baumstein, and M. D. Lacasse, Fast full-wavefield seismic inversion using encoded sources Geophysics (2009) attempted to reduce this noise problem by changing the encoding scheme that is used in the inversion in each inversion iteration. However, this results in the increase number of iterations needed by the simultaneous-source encoded FWI approach. Several techniques for estimating a subsurface electromagnetic model by iteratively minimizing the difference between observed and simulated data are described in T. M. Habashy and A. Abubakar, A general framework for constrained minimization for the inversion of electromagnetic measurements, Progress in Electromagnetics Research, PIER 46, pp. 265-312, (2004) and the references therein. 
     SUMMARY OF THE INVENTION 
     The inventors have recognized that part of a survey data set can be redundant or lacking in sensitivity to relevant reservoir model parameters. For example, a non-optimally designed survey can result in such redundancy or lack of sensitivity. The present invention is predicated in part on recognition that the computational cost of inverting a set of data can be mitigated by detecting and removing from the calculation data that is redundant or lacking sensitivity. 
     In accordance with an aspect of the invention, a computer program product comprises a computer usable medium having a computer readable program code embodied therein, said computer readable program code adapted to be executed to implement a method of compression for inversion of survey data associated with an array of physical transmitters and an array of physical receivers, said method comprising: calculating a compressed data set in which a number of representative transmitters, receivers, or some combination thereof is reduced by calculating a weighted sum of physical transmitters, physical receivers, or both; performing inversion; and outputting a result in tangible form. 
     In accordance with another aspect of the invention, an apparatus comprises: an array of physical transmitters; an array of physical receivers; and one or more processing devices that function alone or in combination to cause the transmitters to produce excitations that are introduced to a subterranean formation and detected by the receivers, and to process data from the receivers by calculating a compressed synthetic data set in which a number of representative transmitters, receivers, or some combination thereof is reduced by calculating a weighted sum of physical transmitters, physical receivers, or both, perform inversion, and output a result in tangible form. 
     Advantages of embodiments of the invention include reduced computational cost and computer memory storage. This advantageously reduces the time required to produce a desired result from a survey and the size of the model that can be processed. Another advantage is that the invention is model-independent in the sense that it does not rely on knowledge of a reservoir model. Other advantages may be apparent in view of the detailed description and figures. 
    
    
     
       BRIEF DESCRIPTION OF THE FIGURES 
         FIG. 1  illustrates apparatus for analyzing a subterranean formation to facilitate recovery of natural resources. 
         FIG. 2  is a flow diagram that illustrates a method for inversion using source-receiver compression. 
     
    
    
     DETAILED DESCRIPTION 
       FIG. 1  illustrates apparatus for surveying a subterranean formation to facilitate recovery of natural resources. An analyzer and control unit  104  is placed in communication with at least one array of excitation sources (transmitters)  106  and at least one array of receivers  108 . The invention is not limited to any particular type of excitation source and may be used with excitation sources including but not limited to electromagnetic (EM), seismic, sonic, gravimetric and pressure. Depending on how the formation is to be surveyed, boreholes  100 ,  102  may be present in the formation. However, the invention is not limited to any particular transmitter-receiver placement and may be used for measurements including but not limited to surface-to-surface, surface-to-borehole, borehole-to-surface, cross-well, and single well. Consequently, arrays transmitters and receivers may be disposed at the surface, in one or more boreholes or wells, or combinations thereof. The analyzer and control unit may include processing hardware and computer-readable memory for storing and operating computer program code that controls the transmitters, receives survey data via the receivers, stores the received data in tangible form, and processes the received data to produce one or more models that are stored in tangible form. Although depicted with a single block, functions of the analyzer and control unit  100  may be performed by different physical devices which may be operated at different locations, e.g., a device that gathers data at the worksite and another device that processes the data at a data processing facility. 
     Referring now to  FIGS. 1 and 2 , a set of data for an array of transmitters and receivers can be arranged in the form of a matrix  D , where the element D ij  corresponds to the measured response at the j-th receiver due to the i-th transmitter. The data matrix  D  is a rectangular matrix of dimension M×N where M is the total number of transmitters and N is the total number of receivers. The case of interest in this description is for M≧2 and N≧2. Through a singular-value decomposition scheme, the data matrix  D  can be expanded as follows: 
                       D   _     _     =       ∑     k   =   1     L     ⁢           ⁢       λ   k     ⁢       v   _     k     ⁢       u   _     k   †                 (   1   )               
where † indicates complex conjugate transpose and L=min{M,N}.  v , ū and λ are the left-singular vectors, the right singular vectors and the singular values of the data matrix  D , respectively. The following orthonormality relationships hold:
 
   v     i   †   ·  v     j =δ ij   (2)
 
 ū   i   †   ·ū   j =δ ij   (3)
 
From eqs. (1), (2) and (3):
 
   D ·ū   k =λ k     v     k   (4)
 
   D ·  v     k =λ k   ū   k   (5)
 
and from eqs. (4) or (5) and eqs. (2) and (3):
 
                           λ   k     =       ⁢         v   _     k   †     ·       D   _     _     ·       u   _     k                   =       ⁢         u   _     k   †     ·         D   _     _     †     ·       v   _     k                     (   6   )               
Construction of an initial matrix is illustrated in step  200 .
 
     Source Compression 
     In view of the above analysis, if the receiver domain matrix  R =  D   † ·  D  is full rank, then  R  will be represented by all its eigenvectors {ū i , i=1, 2, . . . , N}, since their corresponding eigenvalues are non-zero. This will be the case when all receiver measurements provide linearily independent data. In the case when some of the data are redundant,  R  will be rank-deficient and can then be represented by the eigenvectors {ū i , i=1, 2, . . . , L} that have non-zero eigenvalues, where L&lt;N. From eq. (5): 
                               u   _     k     =       ⁢       1     λ   k       ⁢           D   _     _     †     ·       v   _     k           ,           ⁢     k   =   1     ,   2   ,   …   ⁢           ,   L                 =       ⁢       1     λ   k       ⁢       ∑     i   =   1     M     ⁢           ⁢       v     k   ,   i       ⁢       r   _     i   *             ,           ⁢     k   =   1     ,   2   ,   …   ⁢           ,   L                 (   7   )               
where  r   i  is the i-th row of the data matrix  D  whose elements are the receiver measurements corresponding to the i-th transmitter. Equation (7) can be interpreted as follows: v k,i  provides the transmitter weights (in terms of the unit vectors  v   k ) needed to compress the N redundant receiver measurements into L (L&lt;N) linearly independent (non-redundant) measurements corresponding to the L eigenvectors ū k  of the receiver domain matrix  R  that have non-zero eigenvalues. If  J   i  (i=1, 2, . . . , M) represents the moment of the i-th transmitter used in the survey, a compressed set of L (L&lt;M,N) sources,  j   k , given by a weighted sum of  J   i , according to eq. (7), can be constructed as follows:
 
                         j   _     k     =       ∑     i   =   1     M     ⁢           ⁢       v     k   ,   i     *     ⁢       J   _     i           ,     
     ⁢     k   =   1     ,   2   ,   …   ⁢           ,   L           (   8   )               
Consequently, it is possible to avoid running forward modeling code M times corresponding to the M physical sources that were used in the survey. More particular, it is only necessary to run the forward modeling code L times (L&lt;M,N) corresponding to the L compressed sources that are synthesized from the physical sources according to eq. (8). Construction of a set of L sources is depicted in step  202 . Running the forward modeling code L times is depicted in step  204 .
 
     Receiver Compression 
     If the transmitter domain matrix  S =  D ·  D   †  is full rank, then  S  will be represented by all its eigenvectors {  v   i , i=1, 2, . . . , M}, since their corresponding eigenvalues are non-zero. This will be the case when all transmitters provide linearily independent data. In the case when some of the data are redundant,  S  will be rank-deficient and can then be represented by the eigenvectors {  v   i , i=1, 2, . . . , L} that have non-zero eigenvalues, where L&lt;M. From eq. (4): 
                               v   _     k     =       ⁢       1     λ   k       ⁢         D   _     _     ·       u   _     k           ,           ⁢     k   =   1     ,   2   ,   …   ⁢           ,   L                 =       ⁢       1     λ   k       ⁢       ∑     i   =   1     N     ⁢           ⁢       u     k   ,   i       ⁢       c   _     i             ,           ⁢     k   =   1     ,   2   ,   …   ⁢           ,   L                 (   9   )               
where  c   i  is the i-th column of the data matrix  D  whose elements are the measurements at the i-th receiver for all M sources. Equation (9) can be interpreted as follows: u k,i  provides the receiver weights (in terms of the unit vectors ū k ) needed to compress the M redundant transmitter firings into L (L&lt;M) linearly independent (non-redundant) firings corresponding to the L eigenvectors  v   k  of the transmitter domain matrix  S  that have non-zero eigenvalues. The vector of measured receiver response, {  c   i , i=1, 2, . . . , N}, for a particular i-th receiver incorporating all transmitter firings, can then be compressed into a set of L (L&lt;M,N) data sets given by a weighted sum of  c   i , according to eq. (9), as follows:
 
                         C   _     k     =       ∑     i   =   1     N     ⁢           ⁢       u     k   ,   i       ⁢       c   _     i           ,     
     ⁢     k   =   1     ,   2   ,   …   ⁢           ,   L           (   10   )               
Compressing the data into a set of L data points is depicted in step  206 . Consequently, that rather than using N row receiver responses in the inversion that correspond to the N physical receivers that were used in the survey, it is only necessary to use L (L&lt;M,N) compressed data sets constructed from the row data according to eq. (10). Performing the inversion is depicted in step  208 .
 
     Overall Data Compression Strategy 
     In view of the description above it will be appreciated that an overall data compression strategy includes constructing a compressed set of L (L&lt;M,N) sources,  j   k , given by a weighted sum of the M physical sources,  J   i , that were used in the survey, according to eq. (8), as follows: 
                         j   _     k     =       ∑     i   =   1     M     ⁢           ⁢       v     k   ,   i     *     ⁢       J   _     i           ,     
     ⁢     k   =   1     ,   2   ,   …   ⁢           ,   L           (   11   )               
By doing so, the forward modeling code need only be run L times instead of M times. After computing the N receiver responses corresponding to the N physical receivers that were used in the survey, the data is further compressed into a set of L (L&lt;M,N) data points using the weighted sum of eq. (10):
 
                         C   _     k     =       ∑     i   =   1     N     ⁢           ⁢       u     k   ,   i       ⁢       c   _     i           ,     
     ⁢     k   =   1     ,   2   ,   …   ⁢           ,   L           (   12   )               
In this case  c   i  (i=1, 2, . . . , N) is an L×1 vector corresponding to the measured receiver responses for a particular i-th receiver incorporating all of the L transmitter firings corresponding to the L compressed sources, {  j   i , i=1, 2, . . . , L}. By doing so, the number of data points used in the inversion is reduced from N to L for every source firing. Hence, as a result of both the transmitter and receiver compressions, the total number of data points that will be used in the inversion is reduced to L×L instead of M×N. The reduced data matrix, denoted by  d  is given by:
 
   d =  V     †   ·  D ·  U     (13)
 
where  V  is an M×L matrix whose columns are the vectors  v  and  U  is an N×L matrix whose columns are the vectors ū.
 
     The Inversion Strategy 
     Beginning at step  200 , an initial data matrix  D  can be either constructed based on the total field or based on the scattered field obtained from the difference between the total field and the field response of the background medium corresponding to the starting model. The second option is used for the purposes of this description. During the inversion process of step  208 , the minimization is performed on a cost function defined by the differences between the elements of the compressed measurement data matrix  d =  V   † ·  D ·  U  and those of the compressed simulated data  d   s  obtained employing the data compression strategy described above that projects the simulated data onto the compressed data space through the weighted sum of eq. (12) using the compressed set of sources,  j   k , defined by eq. (11). After convergence, the data misfit in the uncompressed space of the physical survey (corresponding to the M transmitters and N receivers) is computed as depicted in step  210 . If the data misfit is below the desired tolerance as determined at step  212 , the iteration is stopped as indicated by step  214 . Step  214  may also include outputting or storing a model in tangible form. If the misfit is larger than the desired tolerance as determined at step  212 , the background medium is then set to this model iterate (at which the iteration converged) as indicated by step  216  and the data misfit is then used to construct the new data matrix in step  200 . A new iteration cycle then proceeds as described above and the process is repeated until the data misfit in the uncompressed physical space is below the desired tolerance. 
     As mentioned above, step  214  may include outputting or storing a model in tangible form. The step can alternatively include a survey design where the number of source firings is limited, or where the survey duration is limited, i.e., reduced relative to a standard survey. Reducing the number of source firings reduces the time required to collect the data. The technique can also be used to do a simultaneous-source survey to obtain a maximum sensitivity in the data. 
     In the case when M&lt;&lt;N, the above scheme will lead to an over-compressed receiver data which may cause the iteration to slow-down. This can be circumvented by augmenting the data matrix,  D , with additional rows that correspond to the data matrix constructed at the previous iteration cycle(s). This expands the size of the basis upon which the receiver array measurements are projected for the purpose of compression. These additional rows correspond to the scattered field data relative to a different background medium (associated with the model inverted from the previous iterations cycle(s)). 
     Similarly, in the case when N&lt;&lt;M, the above scheme will lead to an over-compressed source array which may also cause the iteration to slow-down. This situation can be circumvented by augmenting the data matrix,  D , with additional columns that correspond to the data matrix constructed at the previous iteration cycle(s). This expands the size of the basis upon which the source array is projected for the purpose of compression. These additional columns correspond to the scattered field data relative to a different background medium (associated with the model inverted from the previous iterations cycle(s)). 
     While the invention is described through the above exemplary embodiments, it will be understood by those of ordinary skill in the art that modification to and variation of the illustrated embodiments may be made without departing from the inventive concepts herein disclosed. Moreover, while the preferred embodiments are described in connection with various illustrative structures, one skilled in the art will recognize that the system may be embodied using a variety of specific structures. Accordingly, the invention should not be viewed as limited except by the scope and spirit of the appended claims.