Method for continuous estimation by non-linear signal matching of seismic traces
This is a continuation of application Ser. No. 831,382 filed Feb. 20, 1986, now abandoned. 
    
    
     BACKGROUND OF THE INVENTION 
     1. Field of the Invention 
     The invention relates to the digital signal processing of seismic data to perform certain pattern recognition techniques and, more particularly, but not by way of limitation, it relates to a process for extraction of the continuous RMS and interval velocities by examining the moveout information from common depth point traces at different offsets as measured by non-linear signal matching of selected individual traces. 
     2. Description of the Prior Art 
     The prior art methods for velocity analysis or estimation of continuous velocity determine measures of coherency, such as semblance as a function of time and velocity, as between selected seismic trace data of known offset relationship, and the velocities are then manually picked from an appropriate printout. A general discussion of the prior procedures appears in the publication &#34;Interval Velocities From Seismic Reflection Time Measurements&#34; - P. Hubral and T. Kray, Society of Exploration Geophysics Press, 1980. The seismic velocities can be estimated by exploiting two pieces of information, the change of reflectivity with offset, and the moveout information with offset. Most of the practical methods attempt to exploit the moveout information with offset due to the difficulties in implementation of the first mentioned procedure, i.e. attempting to establish change of reflectivity with offset. Prior methods compute contour plots of the semblance coherency measures and selected velocity relationships may then be manually picked by visual observation of the contour plots. Quantum jumps are usually made in both time and velocity variables in the computation of such contour plots. Such prior methods require an impractically large amount of computation in order to provide only poor resolution of continuous velocity estimation, and they require a considerable amount of manual involvement in order to estimate velocities for each basement point throughout the data compilation. 
     SUMMARY OF THE INVENTION 
     The present invention relates to a method for determining RMS and interval velocity data utilizing non-linear signal matching of common depth point trace data. The method functions to establish selected beginning, ending and intermediate times on the CDP trace group in preparedness for subsequent non-linear signal matching procedure. The routine then functions to estimate derivatives of inverse of the RMS velocity squared for subsequent computation of RMS velocities by post processing of the non-linear signal matching estimates. Given this data, the procedure then establishes an estimate for interval velocities and proceeds to generation of interpretive display, e.g. an output of the seismic section under investigation with multi-color overlay indicating the RMS or interval velocity data. 
     Therefore, it is an object of the invention to develop a technique for the high resolution estimation of elastic constants for an earth cross-section from multi-component, multi-sensor seismic data. 
     It is also an object of the invention to develop a computationally efficient process for estimation of continuous RMS and interval velocities from seismic trace data. 
     It is still further an object of the present invention to provide a method for determination of seismic trace data velocity relationships that are more accurate and better suited for use in related seismic processing operations such as normal moveout corrections, data migration, and the like. 
     Finally, it is an object of this invention to provide an automatic data processing technique for use with seismic section data that provides more accurate continuous velocity estimation. 
    
    
     Other objects and advantages of the invention will be evident from the following detailed description when read in conjunction with the accompanying drawings which illustrate the invention. 
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 is a graphic illustration in idealized form showing a trace match curve indicative of path optimality with local constraints; 
     FIG. 2 is a general flow diagram illustrating the basic procedures of the present invention; 
     FIG. 3 is an idealized showing of plural seismic traces of known offset and time with time window designations; 
     FIG. 4 is a graphic showing of moveout curves for increasing time and offset; 
     FIG. 5 is a graphic illustration of a match curve of CDP traces at different offsets; 
     FIG. 6 is a graphic depiction illustrating the values of the variable u at the ends and check points of a match curve; 
     FIG. 7 is a flow diagram showing the computational flow steps for the continuous estimation of velocities; 
     FIG. 8 is an idealized illustration showing normal moveout estimation for the case of dipping reflectors; 
     FIG. 9 is a typical data printout for instantaneous slope estimates of a match curve; 
     FIG. 10 is a typical data printout for RMS velocity estimates as obtained by non-linear signal matching of CDP traces; 
     FIG. 11 is a typical data printout of interval velocity estimates as obtained by non-linear signal matching of CDP traces; 
     FIG. 12 is a photographic reproduction of a data printout showing a processed seismic section with RMS velocities indicated in color overlay; and 
     FIG. 13 is a photographic reproduction of trace data for the same seismic section with interval velocities illustrated in color overlay. 
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     In accordance with the stated objectives, the present invention must function to match common depth point traces of selected offsets as derived from a line of prestack data. Since the event spacings of traces at different offset positions will vary non-linearly with time, non-linear signal matching is performed to derive moveout information from which to further estimate the RMS and interval velocities. The velocities within the intervals defined by specific interfaces are called interval velocities. These intervals usually pertain to homogeneous stratum. The RMS velocity is defined as the square root of the average of travel time weighted squared interval velocities. 
     The method first establishes a criterion for signal matching. Consider common depth point traces p j  (t) and q j  (t), l≦j≦J, at offsets x 1  and x 2  of J consecutive basement points, where x 2  is greater than x 1 . The trace q j  (t) is a non-linearly stretched version of the trace p j  (t). The stretch function or match curve μ j  (t) that matches the trace p j  (t) with the trace q j  (t) contains the relative moveout information as between the traces at the respective offsets x 1  and x 2  of basement point j. Let p j  (t) be modeled in terms of q j  (t) as 
     
         P.sub.j (t)=q(μ.sub.j (t))+e.sub.j (t),jεJ      (1) 
    
     where e j  (t) is an additive noise term. Letting the term be gaussian with mean α j  and variance σ j   2  as 
     
         e.sub.j (t)˜N(α.sub.j,σ.sub.j.sup.2)     (2) 
    
     and if it is also independent from sample to sample and from trace to trace, the log likelihood function of the data can be written as ##EQU1## where n j  is the number of points on the match curve μ&#39; j  (t). Differentiating L with respect to α j  and σ j  and equating the resulting expressions to zero, we get ##EQU2## Substitution of equations (4) into equation (3) yields an expression for log likelihood function in terms of μ j  (t). That is, ##EQU3## The estimate of μ j  (t) maximizing equation (5) gives maximum likelihood estimate of the match curve. Since log is a monotonic function of its argument, assuming the match curve is constant over a small number of basement points, the u(t) that minimizes ε 2 , where ##EQU4## yields a maximum likelihood estimate for the match curve μ(t). If n is the number of points on the match curve, and in particular, if we consider only one basement point, the criterion for the maximum likelihood estimate of μ(t) is given by the minimization of ε2, where ##EQU5## 
     FIG. 1 illustrates in idealized form the manner in which computationally efficient estimation of μ(t) or match curve 10 is effected from a consideration of non-linearly related traces 12 and 14, i.e. traces q(t) and p(t), respectively. Thus, if we place trace Q(t) 12 on the y-axis and trace p(t) 14 on the x-axis, and we know the samples of beinning point B 16 and the samples of ending point E 18 that match, then it can be postulated that any match curve 10 that goes through points 16 and 18 will define some stretch of one trace relative to the other. The matching of traces is equivalent to estimating a match curve μ(t) that is optimal with respect to the criterion of equation (7). 
     Estimation of the match curve 10 can also be viewed as finding an optimal path from the beginning point 16 to the ending point 18 by using the principle of path optimality. This may be stated as follows. If the optimal path from point B 16 to point C 18 goes through some intermediate point a j , i.e. one of points 20 somewhere within bounds 22, then the optimal path from point B to point C, includes as a portion of it the optimal path from B to a j , or from point 16 to some point 20. This principle can be used to propagate the path sequentially from the beginning point B 16 to the ending point E 18. If we let ε 2  (B, C) be the minimum value of the criterion for the path from B to C, then by using the principle of path optimality it can be computed in terms of the paths from point 16 to 20 as follows. ##EQU6## 
     The size and the shape of the local search region 22 containing the a j  points 20 depend upon the apriori information that one has about the range of local slopes of the particular match curve. 
     FIG. 2 illustrates the basic steps that are performed in the present data processing method. Initially, flow stage 24 examines the input CDP trace data in order to establish times for beginning points, ending points and intermediate check points of a possible match curve. Flow stage 26 then performs non-linear signal matching of the collection of CDP trace data and this includes the function of estimating derivatives for the inverse of the square of the RMS velocity. The RMS velocities are then computed in flow stage 28 by post processing of the match estimates, and interval velocities are estimated by the procedure of flow stage 30 for subsequent generation of output display at stage 32. 
     Referring again to flow stage 24, selection of the end and check points on the CDP traces is done utilizing a semblance criterion to select the corresponding timings on the subject CDP traces for use as the end and check points in non-linear signal matching processing stage 26. Minimum and maximum velocities and velocity increment are selected by the process operator, and suitable time window is selected for semblance calculations on a near offset trace. 
     FIG. 3 illustrates in an idealized form the moveout of a time window for minimum and maximum velocities on a CDP gather. Thus, CDP traces 1 . . . n are shown in relative offset and, for a given event indication, the time window is defined as between the faster velocity moveout lines 34a, 34b and the slower velocity moveout lines 36a, 36b. Thus, the semblance or factor SMB in a time window may be defined as ##EQU7## where s i  (j) is the i th  sample of the time window on the j th  CDP trace. Let x i  be the offset distance of the i th  CDP trace with x 1  being that of the near offset trace. For any particular velocity V the corresponding time on the traces at offsets x 1  and x i  are related as ##EQU8## 
     A semblance factor is calculated for each velocity in the chosen velocity interval. The velocity corresponding to the maximum semblance value is then taken as an estimate of the velocity at the central point of the time window, and equation (10) is then used to find the corresponding times on the traces for subsequent use in non-linear signal matching. The check points are chosen near the dominant events in order to provide maximum resolution. 
     The process then moves to perform non-linear signal matching as in flow stage 26 (FIG. 2). The timings on the traces in a CDP gather at offsets 0, x 1  and x 2  are illustrated in FIG. 4. The curves represent variation with offset of the trace times at the beginning point T B  as indicated by curve 38 through intermediate times to curve 40 at the ending point T E . If these are the matched times, then the time at any offset can be used to reference the velocity. For example, it may be written 
     
         v(T.sub.20)=v(T.sub.2x.sbsb.1)=v(T.sub.2x.sbsb.2)          (11) 
    
     where v(T 20 ) is the RMS velocity at the zero offset time T 20 . 
     FIG. 5 illustrates the general shape of a match curve 42 that was obtained by putting the trace at offset x 1  along the x-axis and the trace at offset x 2  along the y-axis, where x 2  is greater than x 1 . At shallow depths, the velocity will be lower and hence there will be a greater moveout with the increase of the offset as indicated by the beginning point 44 whose ordinate is much longer than it&#39;s abscissa. Since, in general, the velocity will increase with depth, the moveout will be smaller with increase in time such that the ordinates of the points 46, 48 and 50 along match curve 42 will approach their abscissa. The slope of the match curve 42, starting with a low value, will approach unity (line 51) with the increase in time. 
     The match curve 42 possesses information about the relative travel times of the traces and hence information about the velocities. If the time T x  at offset x and the time T o  at zero offset are matched times, the relationship between T x , T o , the offset distance x and the RMS velocity v(T x ) is given as 
     
         T.sub.x.sup.2 =To.sup.2 +x.sup.2 u(T.sub.x)                (12) 
    
     where 
     
         u(T.sub.x)=(1/v.sup.2 (T.sub.x))                           (13) 
    
     Using equation (12), the times T 1x .sbsb.1, T 2x .sbsb.2 on the trace at offset x 1  can be written as 
     
         T.sub.1x.sbsb.1.sup.2 =T.sub.10.sup.2 +x.sub.1.sup.2 u(T.sub.1x.sbsb.2)(14) 
    
     and 
     
         T.sub.2x.sbsb.1.sup.2 =T.sub.20.sup.2 +x.sub.1.sup.2 u(T.sub.2x.sbsb.2)(15) 
    
     Subtraction of equation (14) from equation (15) yields 
     
         (T.sub.2x.sbsb.1.sup.2 -T.sub.1x.sbsb.1.sup.2)=(T.sub.20.sup.2 -T.sub.10.sup.2)+x.sub.1.sup.2 u(T.sub.2x.sbsb.2)-u(T.sub.1x.sbsb.2)(16) 
    
     Similarly, for the times on the trace at offset x 2   
     
         (T.sub.2x.sbsb.2.sup.2 -T.sub.1x.sbsb.2.sup.2)=(T.sub.20.sup.2 -T.sub.10.sup.2)+x.sub.2.sup.2 u(T.sub.2x.sbsb.2)-u(T.sub.1x.sbsb.2)(17) 
    
     and subtracting equation (17) from equation (16), there derives 
     
         (T.sub.2x.sbsb.1.sup.2 -T.sub.1x.sbsb.1.sup.2)=(T.sub.2x.sbsb.2.sup.2 -T.sub.1x.sbsb.2.sup.2)-(x.sub.2.sup.2 -x.sub.1.sup.2)(u(T.sub.2x.sbsb.2))-u(T.sub.1x.sbsb.2))   (18) 
    
     In equation (18), setting 
     
         T.sub.2x.sbsb.2 →T.sub.2x.sbsb.2 →T.sub.x.sbsb.2 and T.sub.2x.sbsb.2 →T.sub.2x.sbsb.2 →T.sub.x.sbsb.2 
    
     there is obtained ##EQU9## 
     Equation (19) relates the slope of the reciprocal of the square of the RMS velocity to the slope of the match curve, and to the matched times, and to the offset distances of the traces. If the RMS velocity increases with time, the quantity ##EQU10## will be less than zero. At any time (t), the variable u can estimated as ##EQU11## Equation (13) can then be used to obtain estimates of the RMS velocities. The zero offset time can be estimated from the matched times of the traces at offset x 1  and offset x 2  as follows. From the moveout equation there is 
     
         T.sub.x.sbsb.1.sup.2 =T.sub.o.sup.2 +x.sub.1.sup.2 /v.sup.2 (T.sub.x.sbsb.1)(23) 
    
     
         T.sub.x.sbsb.2.sup.2 =T.sub.o.sup.2 +x.sub.2.sup.2 /v.sup.2 (T.sub.x.sbsb.2)(24) 
    
     From equations (23 and 24), it may be derived that ##EQU12## 
     Flow stage 28 (FIG. 2) then functions to effect post processing of the non-linear signal match estimates to derive continuous velocity extraction. The purpose of post processing is to minimize the estimation errors of signal matching and to constrain the estimated velocity values at the end and at the respective check points to the given values. The velocities are estimated in such a way that the abrupt discontinuities in the estimated values are minimized subject to the constraint that the estimated velocities agree with the values given at the ending and all check points. 
     Let there be l-1 check points. The beginning point, the ending point and the intermediate check points are graphically illustrated in FIG. 6 for a typical representation of the curve u(t). Designate t o  as the beginning time, t 1  the ending time, Δt the sampling interval, and u&#39;(t) the non-linear signal matching estimate of the derivative of u(t). 
     Further, if u m .sbsb.2 is the value of u(t) at time (t o  +m 2  Δt); and n i  is set to the number of estimates of u&#39;(t) in the time interval t i-1  to t i  (including t i=1  but excluding t i ); and if u&#39;(t) and u(t) are the estimates of u&#39;(t) and u(t) that are obtained by post processing; then the quantities u&#39;(t) and u(t) may be estimated using the relations ##EQU13## 
     The coefficients a j  are chosen so that the estimates u are constrained to the values given at the end and check points and further constrained so that they are smooth. If ##EQU14## then constraining of the u values at the end and the check points will yield the following set of equations 
     
         HA=v                                                       (29) 
    
     where for -m 1  ≦j≦m 2  and 2≦i≦l, the elements of H and v are given by ##EQU15## 
     The smoothness of the estimated u values can be measured in terms of ε 2 , where ##EQU16## If 
     
         ε.sub.1.sup.2 =ε.sup.2 /(Δt).sup.2   (32) 
    
     then equation (31) can be put into the matrix form as 
     
         ε.sub.1.sup.2 =A.sup.T UA                          (33) 
    
     where the (i, j) th  element of the matrix U is given by ##EQU17## The coefficient vector A that minimizes ε 1   2  of equation (33) subject to the constraints of equation (29) can be obtained as 
     
         A-U.sup.-1 H.sup.T (HU.sup.-1 H.sup.T).sup.-1 v            (35) 
    
     The u values are estimated as ##EQU18## 
     Estimation of the interval velocities is carried out in flow stage 30 (FIG. 2) as the RMS velocity is related to the interval velocities of the respective subterranean earth layers and the respective layer travel times as follows ##EQU19## where V i  and t i  are the interval velocity and one way travel times of earth layer i, and V RMS  is the RMS velocity up to the layer n. From equation (37), the DIX equation, the interval velocity of layer n can be written as ##EQU20## 
     The computations of flow stages 26, 28 and 30 of the flow chart of FIG. 2 are further illustrated in the flow chart of FIG. 7. Thus, the match curve estimate is made in flow stage 52 and the match curve slopes are determined in flow stage 54. The flow stage 28 (FIG. 2) consists of a multiple of individual routines as u&#39;(t) and ηl i  are determined in flow stages 56 and 58. The matrix H is then formed in stage 60, as at equations (29 and 30), and the matrix U is derived in stage 62 as per equations (33 and 34). The co-efficient vector A (equation 35) is determined in flow stage 64, and the RMS velocities are then finally estimated in stage 66. Flow stage 30 then functions to estimate the interval velocities in accordance with the functions of equations (37 through 39). 
     The velocity estimates are considerably affected by the effects of dip in the reflected horizons. FIG. 8 illustrates the relationship and parameters to be considered for normal moveout estimation for dipping strata. Energy from source S at earth surface 70 is directed along ray path 72 to be reflected at dipping bed 74 for return along path 76 to receiver point R on surface 70. Thus, from the relationships of FIG. 8, it can be written 
     
         (vt).sup.2 =4h.sup.2 +x.sup.2 -4hx cos(π/2+θ) 
    
     
         (40) 
    
     Rearranging equation (40), there is obtained ##EQU21## and using equation (41) while proceeding similar to equation (19), it may be stated as ##EQU22## If the changes in the velocities are small compared to the changes in the square of the velocities, the last term of equation (42) can be neglected. Then equation (42) reduces to equation (19). 
     The zero offset time with dip can be estimated as follows. Letting T x .sbsb.1, T x .sbsb.2, and T o  be the matched times; then for the case of dipping beds, using equation (41), it can be written that ##EQU23## From equations (43) and (44), there obtains ##EQU24## Thus, knowing T x .sbsb.1 and T x .sbsb.2, the calculation of T o  involves the angle of dip; however, if we use one trace at offset zero, and the other trace at offset x 1 , the calculations are essentially independent of dip to a first order approximation. 
     FIGS. 9 through 13 illustrate output results of the present method as utilized in a specific seismic application. A group of common depth point data as derived from a marine survey was the subject matter. The seismic survey data consisted of 81 CDP gathers, of forty-eight fold redundancy, and the data was processed to provide estimation of RMS and interval velocities. For the 81 CDP gathers, signal match traces were selected centering on offsets 3 and 24 of each DCP gather, this providing a significant and ample offset spacing. Thus, trace data at offset 3 and the trace data at two offsets on each side thereof were processed on the x-axis, and five offset traces centering on offset 24 were processed on the y-axis in order to establish match curve information for the respective CDP group using non-linear signal matching. The matched times at the beginning of the traces, end of traces, and for each of two intermediate check points, were estimated initially using semblance criterion. 
     FIG. 9 illustrates computer printout information for the instantaneous slope estimates of a typical match curve as determined at flow stage 54 (FIG. 7). The data values are read sequentially horizontally across the rows and vertically down the columns row by row, and it can be generally discerned that the slopes of the match curve data start with a lower value and increase with increasing time approaching unity at the trace ending point. 
     FIG. 10 shows the printout data for the RMS velocity estimates for the same data. The RMS velocities are determined in flow stage 66 which essentially functions around the evaluation of equation (36). Finally, the printout data of FIG. 11 illustrates the interval velocity estimates as determined in flow stage 30 (FIG. 7), such interval velocity data being estimated utilizing equations (38 and 39). 
     The estimated RMS and interval velocities can then be used to generate displays of color overlays on seismic sections. Such displays then allow the interpreter to delineate the potential hydrocarbon traps through deduction. That is, any lithologic unit has a characteristic rock velocity which will remain constant until either the rock composition changes or until its porosity varies. A change in rock composition or porosity will be revealed by a change in transit time, i.e. the interval velocity for the stratum. The individual lithologies or porosity changes can be outlined by contouring the boundaries at the top and base of each member that has a velocity that contrasts with that of an adjacent member. 
     FIG. 12 is a seismic compressional wave section 80a wherein the RMS velocity data is displayed by means of color overlay. The section 80a shows the eighty-one CDP gathers in equal spacing across the survey line from left to right, and this survey passed through a well position 82 which proved to be a gaseous but non-commercial well. The color scale 84 progresses from 5500 feet per second through 7000 feet per second in color gradation as follows: blue-86, light blue-88, green-90, yellow-92, red-94 and indigo-96. Thus, the color overlay for RMS velocities of FIG. 12 shows largely blue-86 in the upper regions above a generally light blue-88 line. The lower regions thereunder range from green-90 to yellow-92. Since the gas zone will usually have a low interval velocity, it can be noted beneath well position 82 that the blue-86 portion extends nearly to the bottom of the seismic section 80, i.e. down to about 2.4 seconds travel time. Keep in mind that seismic event data is also present for visual interpretation with the color representation. 
     FIG. 13 illustrates the same seismic section 80b as it includes a color overlay indicating interval velocities. In this case, the color scale 100 from blue to indigo extends from 5500 feet per second to 9500 feet per second. The color overlay from about 1.5 seconds through 2.6 seconds travel time is again largely blue-86 with increased interval velocity zones showing up as light blue-88 intermixed with greens-90, yellow-92 and reds-94. It also may be noted beneath well position 82 a region of low interval velocity as indicated by blue-86 with a bright spot light blue-88 indication overlying concise events in the vicinity, i.e. at about 1.8 seconds. This is a classic gas zone indicator as high porosity is evidenced. Lower strata shows scattering of higher velocity structure by presence of yellows, reds and indigo. 
     The foregoing discloses a novel method for processing of seismic data whereby velocity analysis is effected relatively faster over a section of seismic data. The present invention enables not only the better and more complete derivation of respective RMS and interval velocities, but also the enhanced interpretive tool consisting of the color overlayed seismic section wherein the velocity relationships are denoted. The method utilizes non-linear signal matching to establish a common match curve relationship that is a basic function of all traces of the section, and it brings a new tool to the geophysicist for evaluating seismic data in order to better identify subsurface structures. 
     Changes may be made in the combination and arrangement of steps as heretofore set forth in the specification and shown in the drawings; it being understood that changes may be made in the embodiments disclosed without departing from the spirit and scope of the invention as defined in the following claims.