Abstract:
Improved finite-difference staggered grid wave propagation systems and methods. One illustrative computer-based wave field simulation method includes: applying at least one signal to a grid of model cells forming a model space, each model cell having stress values associated with stress nodes and velocity values associated with velocity nodes staggered from the stress nodes; and propagating the at least one signal as a wave into the model space by alternately updating the stress values and the velocity values to obtain a time-dependent wave field associated with the at least one signal. The stress value updating includes, for each model cell: determining spatial derivatives of the velocity values for the model cell; interpolating the spatial derivatives to multiple stress nodes within the model cell; and, for each stress node within the model cell, combining the spatial derivatives associated with that stress node to update at least one stress value associated with that stress node.

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS 
       [0001]    Not applicable. 
       STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH OR DEVELOPMENT 
       [0002]    Not applicable 
       BACKGROUND 
       [0003]    1. Field of the Invention 
         [0004]    This invention relates generally to the field of geophysical exploration for hydrocarbons. More specifically, the invention relates to a method of faster and more efficient grid processing in seismic modeling. 
         [0005]    2. Background of the Invention 
         [0006]    Scientists and engineers often employ geophysical surveys for exploration and engineering projects. Geophysical surveys can provide information about underground structures, including formation boundaries, rock types, and the presence or absence of fluid reservoirs. Such information greatly aids searches for water, geothermal reservoirs, and mineral deposits such as hydrocarbons and ores. Oil companies in particular often invest in extensive seismic and electromagnetic surveys to select sites for exploratory oil wells. 
         [0007]    Geophysical surveys can be performed on land or in water. As indicated in the example survey of  FIG. 1 , an energy source  102  near the region of interest  104  generates waves  106  that propagate into the region of interest and reflect from internal features such as bed boundaries. Some of the reflected wave  108  energy eventually reaches an array of receivers  110  on the surface  112 . A recording system  114  captures the received signals for storage and processing. The process is repeated with many different source positions and optionally with different receiver positions. Although various methods exist for converting the received wave signals into an image of the subsurface structure, the most popular such techniques employ finite difference wave field modeling, a process that propagates waves forward or backward in time using discrete time steps and fast approximations of wave function derivatives. With small enough time steps and accurate enough approximations, the finite difference wave field modeling process provides an efficient, high fidelity solution. Nevertheless, the finite difference modeling process imposes a substantial computational demand. 
         [0008]    Many efforts have been made to satisfy the computational demands of the finite difference wave field modeling process. Many combined hardware-software solutions have been implemented in the attempt to maximize runtime efficiency. In addition to exploiting different parallelization strategies offered by different hardware systems, the software is often customized to employ different wave-propagation kernels using, e.g., different discretizations, wave equation approximations, and propagation strategies. Attempts to achieve faster processing (i.e., shorter run-times) with acceptable fidelity while honoring the memory, bus bandwidth, and other system limitations, have led to innumerable kernel variations which are disclosed in literature. Yet the computational burdens remain undesirably high. 
       BRIEF SUMMARY 
       [0009]    Accordingly, there are disclosed herein finite-difference staggered grid wave propagation systems and methods having a reduced computational burden achieved by re-ordering of operations within the propagation kernel. One illustrative computer-implemented wave field simulation method includes: applying at least one signal to a grid comprising a plurality of model cells, the grid forming a model space, each model cell having a plurality of stress values associated with a plurality of stress nodes and a plurality of velocity values associated with a plurality of velocity nodes staggered from the stress nodes; and propagating the at least one signal as a wave into the model space by alternately updating the stress values and the velocity values to obtain a time-dependent wave field associated with the at least one signal. The stress value updating includes, for each model cell: determining spatial derivatives of the velocity values for the model cell; interpolating the spatial derivatives to multiple stress nodes within the model cell; and, for each stress node within the model cell, combining the spatial derivatives associated with that stress node to update at least one stress value associated with that stress node. 
         [0010]    The computational savings need not be restricted to wave propagation. One illustrative computer-based method provides a combined convolution-interpolation operation for applying values in a first grid to nodes in a second grid staggered from the first grid, the method including, for each cell in the second grid: obtaining a weighted sum of values in a window around a corresponding cell in the first grid using coefficients of a convolution filter; and interpolating the weighted sum to each of multiple nodes in the cell of the second grid. 
         [0011]    A disclosed non-transitory information storage medium further provides software for a seismic imaging method that includes: applying at least one seismic source or receiver signal to a grid comprising a plurality of model cells, the grid forming a model space, each model cell having a plurality of stress values associated with a plurality of stress nodes and a plurality of values associated with a plurality of velocity nodes staggered from the stress nodes; and propagating the at least one seismic source or receiver signal as a wave into the model space by alternately updating the stress values and the velocity values to obtain a time-dependent wave field associated with the at least one seismic source or receiver signal; repeating said applying and propagating operations to obtain multiple time-dependent wave fields; analyzing a counterpart source and receiving one or more wave fields to obtain matching information; stacking the matching information to image a subsurface region represented by the model space; and displaying the image. The stress value updating includes, for each model cell: determining spatial derivatives of the velocity values for the model cell; interpolating the spatial derivatives to multiple stress nodes within the model cell; and for each stress node within the model cell, combining the spatial derivatives associated with that stress node to update at least one stress value associated with that stress node. 
     
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         [0012]    In the accompanying drawing sheets: 
           [0013]      FIG. 1  is a schematic depiction of a seismic survey. 
           [0014]      FIG. 2  is a block diagram of a seismic survey system having a finite difference wave field simulator. 
           [0015]      FIG. 3  is a graph of illustrative digitized receive signal traces. 
           [0016]      FIG. 4  is a view of a geo-modeler employing an illustrative modeling system. 
           [0017]      FIG. 5  is a block diagram of an illustrative hardware platform suitable for geophysical modeling. 
           [0018]      FIG. 6A  is a data space representing a subsurface region of interest. 
           [0019]      FIG. 6B  is a cell of the data space having staggered grid values. 
           [0020]      FIG. 6C  is a cell of the data space including interpolated velocity values. 
           [0021]      FIG. 7  is a flow diagram of an illustrative subsurface imaging method. 
           [0022]      FIG. 8  is a flow diagram of an illustrative migration kernel of a first type. 
           [0023]      FIG. 9  is a flow diagram of an illustrative migration kernel of a second type. 
       
    
    
       [0024]    It should be understood, however, that the specific embodiments given in the drawings and detailed description thereto do not limit the disclosure. On the contrary, they provide the foundation for one of ordinary skill to discern the alternative forms, equivalents, and modifications that are encompassed together with one or more of the given embodiments in the scope of the appended claims. 
       DETAILED DESCRIPTION 
       [0025]    The seismic survey of  FIG. 1  may be performed using the system of  FIG. 2 . An array of seismic receivers  110  convert seismic waves to electrical signals that are amplified and digitized. (Illustrative signal waveforms are shown in  FIG. 3 .) A recording system  114  collects the digitized signals via a bus  202  or other communications pathway and stores the digitized signals on a non-transitory information storage medium for later processing. Typically, each digitized signal is associated with parameters  204  such as a receiver location and a shot location and such other information as the system designer deems worthwhile. Recording system  114  may perform some initial processing to filter and/or compress the data, and in at least some cases, to perform quality control. 
         [0026]    The recording system  114  provides the seismic survey data via the Internet or some other communications mechanism  206  to a data processing center  208  having sufficient computational resources for the imaging process. The data processing center includes one or more computers that may use finite difference wave field modeling to perform migration and thereby convert the recorded seismic signals into a three dimensional map or “image” of the subsurface structure which can be displayed on a monitor  210  and stored in persistent storage for later use. 
         [0027]    In the data processing center  208 , processing of the seismic survey data is coordinated by a geo-modeler such as that shown in  FIG. 4 . The geo-modeler may employ a user interface  402  of a workstation  404  to initiate processing of the seismic data and to view and analyze the resulting seismic image. The workstation  404  is part of the hardware platform of a subsurface modeling system such as that shown in  FIG. 5 . (Other hardware platforms are also available and suitable, including supercomputers, massively parallel processor networks, and heterogeneous architectures that include GPUs or other dedicated co-processors.) 
         [0028]    The illustrative hardware platform couples the workstation  404  to one or more multi-processor computers  506  via a local area network (LAN)  504 . The one or more multi-processor computers  506  are in turn coupled via a storage area network (SAN)  508  to one or more shared storage units  510 . Using the personal workstation  404 , the geo-modeler is able to load seismic survey data into the system, to configure and monitor the processing of the seismic survey data and to retrieve the results from the system, optionally in the form of volumetric images. 
         [0029]    Personal workstation  404  may take the form of a desktop computer with a display that shows graphical representations of the input and result data, and with a keyboard that enables the user to move files and execute processing software. LAN  504  provides high-speed communication between multi-processor computers  506  and with personal workstation  404 . The LAN  504  may take the form of an Ethernet network. 
         [0030]    Multi-processor computer(s)  506  provide parallel processing capability to enable suitably prompt processing of the input data to derive the results data. Each computer  506  includes multiple processors  512 , distributed memory  514 , an internal bus  516 , a SAN interface  518 , and a LAN interface  520 . Each processor  512  operates on allocated tasks to solve a portion of the overall problem and contribute to at least a portion of the overall results. Associated with each processor  512  is a distributed memory module  514  that stores application software and a working data set for the processor&#39;s use. Internal bus  516  provides inter-processor communication and communication to the SAN or LAN networks via the corresponding interfaces  518 ,  520 . Communication between processors in different computers  506  can be provided by LAN  504 . 
         [0031]    SAN  508  provides high-speed access to shared storage devices  510 . The SAN  508  may take the form of, e.g., a Fibrechannel or Infiniband network. Shared storage units  510  may be large, stand-alone information storage units that employ magnetic disk media for nonvolatile data storage. To improve data access speed and reliability, the shared storage units  510  may be configured as a redundant disk array (“RAID”). 
         [0032]    Whether the hardware platform takes the illustrated form or another suitable form, an issue common to all such platforms is the comparatively low data retrieval rate versus the rate at which the processors can perform computations. That is, most computers designed for processing large volumes of data have a memory bottleneck that limits the rate at which the processing can occur. Thus, although the computational burden associated with many parallelizable processes can be reduced by increasing the available memory, such solutions may fail to complete the processing any faster. It is often a better approach to reduce the volume of memory accesses required, even if such a strategy increases the computational burden carried by the processors. 
         [0033]    To enable the hardware platform to simulate the wave field and construct the image, the software models the region of interest as a two-dimensional or three-dimensional space divided into a grid of cells  602 .  FIG. 6A  shows an example of a 3D space. Each cell has a representative set of values that commonly are associated with a single point within the cell. However, as disclosed in literature (see, e.g., R. W. Graves, “Simulating Seismic Wave Propagation in 3D Elastic Media Using Staggered-Grid Finite Differences”, Bull. Seismological Soc. Am., v86(4) pp 1091-1106, August 1996), the use of a staggered grid yields certain efficiency advantages, at least in isotropic and orthotropic media. Before delving into the details of the staggered grid, an overview of the imaging process is helpful. 
         [0034]      FIG. 7  is a flowchart of an illustrative seismic imaging method that may be performed by the computers in the data processing center. In block  702 , the computers obtain the survey data, including the digitized signals and associated source and receiver position information. An initial model of the subsurface may also be obtained in this block, e.g., to specify an estimated density and elasticity as a function of position. Often a uniform density model is employed as the starting point and gradually refined during the imaging process. 
         [0035]    In block  704 , the computers simulate the evolution of the wave field (the spatial distribution of two interrelated wave parameters such as force, pressure, or stress; and displacement, velocity, or strain) generated by the source by migrating the source signal (or a suitable substitute) into the model space from the source position. That is, the source generates stress and particle velocity fields that propagate outward from the source position in accordance with the wave equation. A finite difference model simulates this propagation to determine the source wave field, i.e., the stress and particle velocity at each point in the model space, as a function of time. The propagation kernel is described further below. 
         [0036]    In block  706 , a similar operation is applied to the receive signals. Recognizing that the recorded signals represent waves that have reached the receiver locations, the corresponding particle velocities or stresses are propagated backward in time into the model space in accordance with the wave equation. The same finite difference propagation kernel as before can be used, albeit with a reversed time index, to determine the receive wave field as a function of time. (With the time reversal, receive waveforms get treated like source waveforms.) 
         [0037]    Because reflectors in the earth converted the source wave field into the receive wave field during the survey process, these two wave fields will match at the reflection points. Accordingly, the source and receive wave fields are analyzed as a function of time and position in block  708  to identify those regions that match and hence indicate a reflection point. This process is repeated for each seismic shot and the matching information is added up (“stacked”) to image the whole region of interest. 
         [0038]    In block  710 , the image is analyzed to identify inaccuracies in the estimated formation parameters (density and/or elasticity). Since such inaccuracies produce certain patterns in the image data, they enable the initial estimates to be refined. Blocks  704 - 710  can be iteratively repeated to gradually refine and improve the model parameters (and the resulting images). Once the model parameters converge, the computers can display the image in block  712 . Actions can then be taken based on the displayed image, including steering decisions, landing decisions, and completion decisions. 
         [0039]    A particular point of note is that the migration process may need to be repeated many times before convergence is reached, and each iteration of the migration process requires the solution of the propagation equation at each point in the model space for each time step for each shot. Hence the finite difference propagation kernel may be repeated hundreds of billions of times. Though the process includes a significant degree of parallelism which can support implementation on parallel processors, the computational burden remains high. Even small gains in the kernel&#39;s computation efficiency can yield substantial savings in execution time. Benefits of such efficiency gains in the migration process are not limited to seismic imaging. For example, migration can also be employed to simulate wave fields for a wide range of physics (acoustic waves, electromagnetic waves, fluid dynamics) and application contexts. 
         [0040]    As mentioned previously, significant efficiency gains are achieved through the use of a staggered grid. An illustrative staggered-grid is illustrated in  FIG. 6B , where a single cell  602  from the 3D space of  FIG. 6A  has normal (compressional) stress values σ xx , σ yy , σ zz , associated with the corner nodes, shear stress values σ xy , σ xz , σ yz , associated with respective face nodes, and velocity values v x , v y , v z , associated with the edge nodes. The edge nodes are offset (“staggered”) from the corner nodes by a half-cell spacing along one axis, and the face nodes are offset from the corner nodes by half-cell spacings along two axes. In isotropic or even orthotropic media, the finite difference operators that interrelate velocity and stress values are naturally centered on the node values they are used to modify, so no additional interpolation is needed. 
         [0041]    However, as explained in J. W. Rector, G. M. Hoversten, K. T. Nihei, “Seismic Modeling Engines Consortium”, Lawrence Berkley National Laboratory, Dec. 3, 2003, additional interdependencies exist in fully generalized anisotropic media. As a consequence, the propagation kernel must be modified to do one of two things. Either the kernel can employ a full stencil (such as that used in a rotated staggered grid model), which requires four times as much memory, or the kernel can employ interpolation to provide the centering as needed. Specifically, as shown by the underlined values in  FIG. 6C , the velocity values may be interpolated to the other edge nodes and to the body center node of the cell  602  to support updating of the stress values. Though such interpolation imposes an additional computational burden, it is preferable to the full stencil approach because it avoids exacerbating the memory bandwidth bottleneck and it avoids certain mesh drift instabilities. 
         [0042]    Using tensor notation, the momentum conservation and constitutive equations for the particle velocities v i  and stresses σ ij  are: 
         [0000]    
       
         
           
             
               
                 
                   
                     ρ 
                      
                     
                       
                         ∂ 
                         
                           v 
                           i 
                         
                       
                       
                         ∂ 
                         t 
                       
                     
                   
                   = 
                   
                     
                       
                         ∂ 
                         
                           σ 
                           ij 
                         
                       
                       
                         ∂ 
                         
                           x 
                           j 
                         
                       
                     
                     + 
                     
                       
                         f 
                         i 
                       
                        
                       
                           
                       
                        
                       and 
                     
                   
                 
               
               
                 
                   ( 
                   
                     1 
                      
                     a 
                   
                   ) 
                 
               
             
             
               
                 
                   
                     
                       
                         ∂ 
                         
                           σ 
                           ij 
                         
                       
                       
                         ∂ 
                         t 
                       
                     
                     = 
                     
                       
                         
                           c 
                           ijkl 
                         
                         2 
                       
                        
                       
                         ( 
                         
                           
                             
                               ∂ 
                               
                                 v 
                                 k 
                               
                             
                             
                               ∂ 
                               
                                 x 
                                 l 
                               
                             
                           
                           + 
                           
                             
                               ∂ 
                               
                                 v 
                                 l 
                               
                             
                             
                               ∂ 
                               
                                 x 
                                 k 
                               
                             
                           
                         
                         ) 
                       
                     
                   
                   , 
                 
               
               
                 
                   ( 
                   
                     1 
                      
                     b 
                   
                   ) 
                 
               
             
           
         
       
     
         [0000]    where ρ is density, i, j, k, l are indices ranging over the spatial axes X, Y, Z, x i  is distance along the given axis, f i  is force applied along the given axis (e.g., as part of a source signal), and c ijkl  is an elastic constant relating stress to strain. Pursuant to convention, implicit summation occurs over those indices that appear only on the right side of the equation, i.e., the right side of equation (1a) has an implicit summation over j ranging from 1 to 3 (axes X to Z), and the right side of equation (1b) has an implicit summation over k and an implicit summation over l, each ranging from 1 to 3. 
         [0043]    It is equation (1b) that, when reduced to finite difference form, necessitates interpolation to center the differentiated fields on the stress value nodes in a staggered grid. Taking as an example the finite difference form of the equation for σ xy , we have: 
         [0000]      σ′ xy =σ xy   +c   xyxx   D   x ( v   x )+ c   xyyy   D   y ( v   y )+ c   xyzz   D   z ( v   z )+ c   xyxy ( D   x ( v   y )+ D   y ( v   x ))+ c   xyxz ( D   x ( v   z )+ D   z ( v   x ))+ c   xyyz ( D   y ( v   z )+ D   z ( v   y )),  (1c)
 
         [0000]    where D i ( ) is the finite difference approximation of a derivative along the given axis. For the finite differences to be suitably centered on the face node for σ xy , each of the velocity field components must be interpolated. 
         [0044]    While interpolation can be a computationally inexpensive operation, each of the three velocity field components are being interpolated to complete the information at four different nodes, requiring nine interpolations within the cell. Each finite difference operator, however, operates on values from multiple cells. The number of cells (and hence the number of interpolations needed) varies with the order of the chosen finite difference operator. A finite difference operator of order N operates on values from N+1 cells. Thus, the nine finite differences for updating each of the four stress nodes of a cell require 9 interpolations to be performed in 3N+1 cells for a total of 27N+9 interpolations. For an order 3 finite difference, there are 90 interpolations needed to update one cell. Higher orders may be used for increased fidelity. 
         [0045]    While opportunities exist for re-use of interpolated values to varying degrees, additional memory is required to store those values until they can be re-used. (At the extreme, the interpolations may be performed across the whole volume prior to This trade-off has not generally been considered acceptable. Thus there continues to be 90 interpolations performed to update each cell. 
         [0046]    With this in mind,  FIGS. 8 and 9  are flow diagrams illustrating different approaches to the propagation kernel which may be used to implement migration blocks  704 ,  706  ( FIG. 7 ). The migration blocks  704 ,  706 , parallelize the migration operation by distributing pieces of the problem (e.g., different portions of the model space) to different processors. Beginning in block  802 , each processor receives its task(s), which may represent a portion of the model space, and stage the data in a local memory for expedited paging to and from the processor cache. Four indices may be used to track progress through the propagation process, the indices representing the three spatial dimensions and a time dimension. These indices are initialized in block  804 . 
         [0047]    The first data blocks are paged into cache in block  806 . This example assumes that each block of data represents an x-z slice of the data volume corresponding to a given y-index value. The set of blocks  808 - 824  are performed as the propagation kernel iterates through a given x-z slice to update the stress values. Blocks  806  and  826 - 830  cause the propagation kernel to iterate through each of the slices in the data volume to complete the stress value update for a given time step. Blocks  831 - 834  cause the propagation kernel to iterate through subsequent time steps. 
         [0048]    The set of cells  840  represent the operations performed for a given cell to update the stress values associated with that cell. In block  808 , the processor performs bi-linear interpolation of the velocity field components to the staggered velocity nodes (see  FIG. 6C ) in each of N+1 cells along the current spatial derivative axis. In block  810 , the finite difference along that axis is determined for each of the velocity field components at each of the stress value nodes in the current cell. In block  612 , the processor checks if more spatial derivative axes exists, and if so, the next axis is chosen in block  813  and the processor returns to block  808 . The axes all intercept in the current cell, so the interpolation of values for the current cell need not be repeated—only the values for the additional N cells along the selected axis need to be interpolated. 
         [0049]    Once the full set of spatial derivatives at each stress node in the current cell has been determined, in block  816  the processor combines them with the respective elastic constants to update the stress values for the current cell. See equations (1b), (1c). 
         [0050]    In block  818 , the processor determines whether the x-index has reached the end of its range, and if not the value is incremented in block  820  and control returns to block  808 . Otherwise, in block  822  the processor determines whether the z-index has reached the end of its range, and if not, the value is incremented in block  824 , the x-index is reset, and control returns to block  808 . Otherwise, the current x-z slice has been completed. In block  826 , the updated stress values are stored. In block  828 , the processor determines whether the y-index has reached the end of its range, and if not, the value is incremented in block  830 , the x and z indices are reset, and control returns to block  806  to load the next data block. Otherwise, the stress values have been updated for the full data volume, and in block  831  the processor iterates through the data volume to update the velocity field in accordance with equation (1a). Once that has been completed, a full time step for the data volume has been completed. In block  832 , the processor determines whether the time index has reached the end of its range, and if not, the value is incremented in block  834 , the x, y, z indices are reset, and control returns to block  806 . 
         [0051]    In contrast to the propagation kernel of  FIG. 8 , the improved propagation kernel of  FIG. 9  re-orders the interpolation and derivative operations. Thus, within the set of operations performed for each cell (box  850 ), the propagation kernel begins by finding the velocity components&#39; spatial derivative (finite difference) along a given axis in block  852 . For the x-axis, block  852  yields the derivatives D x (v x ), D x (v y ), D x (v z ). Due to the absence of interpolation at this step, the centering of these derivatives will vary. Derivatives D x (v x ), D y (v y ), D z (v z ), will center on the corner node, while the three face nodes will each have two associated derivatives. For example, the derivatives D x (v y ) and D y (v x ) will center on the face node associated with σ xy . Due to the arrangement of the staggered nodes, no interpolation is required for this step. Then, in block  854 , the processor performs bilinear interpolation within the current cell only to determine the derivative values at all of the stress nodes. As 4 stress nodes exist and each derivative calculated in block  852  is already centered on one of them, each derivative needs to be interpolated to three other nodes. After all spatial axes have been accounted for, the processor need only perform a total of 9 bilinear interpolations, irrespective of the finite difference order N. Thus, this re-ordering results in a savings of 27N interpolations. 
         [0052]    We note that the finite difference calculations can be considered as a particular type of filter or other convolution operation. The benefits which have now been highlighted above are similarly achievable for other filtering or convolution operations being performed on interpolated values. That is, a significant computational savings is achieved by interchanging the order of the interpolation and convolution operations. 
         [0053]    Numerous other variations and modifications will become apparent to those skilled in the art once the above disclosure is fully appreciated. For example, the details of the hardware and software implementations are subject to a wide range of variation, alternative embodiments, and different forms of optimization or approximation. Moreover, though the foregoing disclosure was provided in the context of a seismic imaging system, wave propagation may be simulated in other contexts including acoustic fields, electromagnetic fields, and fluid flow mechanics. It is intended that the following claims be interpreted to embrace all such variations and modifications.