Expedient processing and waveform inversion of seismic data

Abstract

A method for expediently processing and inverting elastic wave data to reduce the amount of data and to determine a physical properties model of the material medium and a source properties model. The data are processed to generate waveforms containing the phase difference between compressional- and shear-wave arrivals using auto-correlation, cross-correlation, or deconvolution of said data sensed at each of an arrangement of sensors, whereby said lengthy elastic wave data records are reduced substantially in time. Said waveform data are thereafter inverted using waveform inversion by modifying the source term in the equation of motion, wherein the source term is mathematically expressed as a product of time-independent source properties volume defined at every location in space within said material medium and a space-independent source-time function, whereby no prior knowledge of the number of sources, spatial distribution of source location, source amplitude, or source focal mechanism is needed.

Claims

1. A method for compressing the amount of elastic wave data, for eliminating the source-time functions, and for simulating composite impulse responses, comprising the steps of: (a) recording elastic wave data in at least one spatial location that includes at least one directional component of particle displacement, velocity, or acceleration and said elastic data originating from at least one source selected from the group consisting of man-made and natural sources, and (b) using at least one operation selected from the group consisting of auto-correlation, cross-correlation, and deconvolution, on said elastic data, thereby generating waveforms containing the phase difference between compressional-wave mode and slow-shear-wave mode, or any two different elastic wave modes, selected from the group consisting of compressional-wave mode, fast-shear-wave mode, and slow-shear-wave mode, and (c) discarding the acausal part, and (d) truncating the causal part at a predetermined time value ascertained by the maximum expected time interval between the said elastic wave modes originating from said sources, whereby said elastic wave data are compressed, and whereby either the phase or both the amplitude and phase of the source-time functions of sources generating said elastic waves are eliminated, and whereby the said computed waveform at each sensor represents the composite impulse response between the sources and the said sensor.

2. The method in claim 1, wherein said elastic wave data are pre-processed to separate said elastic wave modes.

3. The method in claim 1, wherein the amplitude about zero time-lag of said computed waveforms is suppressed.

4. The method in claim 1, wherein said elastic wave data are pre-processed to generate envelope functions in time, whereby the sign of the source radiation pattern is eliminated.

5. A system for compressing the amount of elastic wave data, for eliminating the source-time functions, and for simulating composite impulse responses, comprising: (a) a data recording system comprising of one or more computers, data recording hardware, and one or more elastic wave sensors that include at least one directional component of particle displacement, velocity, or acceleration, that record elastic wave data in at least one spatial location and said elastic data originating from at least one source selected from the group consisting of man-made and natural sources, and (b) a processor that uses at least one operation selected from the group consisting of auto-correlation, cross-correlation, and deconvolution, on said elastic data, thereby generating waveforms containing the phase difference between compressional-wave mode and slow-shear-wave mode, or any two different elastic wave modes, selected from the group consisting of compressional-wave mode, fast-shear-wave mode, and slow-shear-wave mode, and (c) wherein the processor further discards the acausal part, and (d) wherein the processor further truncates the causal part at a predetermined time value ascertained by the maximum expected time interval between the said elastic wave modes originating from said sources, whereby the system compresses said elastic wave data, and whereby the system eliminates either the phase or both the amplitude and phase of the source-time functions of sources generating said elastic waves, and whereby the system simulates the composite impulse response between the sources and each sensor.

6. The system in claim 5, wherein the processor pre-processes said elastic wave data to separate said elastic wave modes.

7. The system in claim 5, wherein the processor further suppresses the amplitude about zero time-lag of said computed waveforms.

8. The system in claim 5, wherein the processor pre-processes said elastic wave data to generate envelope functions in time, thereby eliminating the sign of the source radiation pattern.

9. A method for efficient imaging of acoustic or processed elastic wave data using the inhomogeneous acoustic wave equation $\begin{array}{cc}m\left(x\right)\frac{{\partial }^{2}u\left(x,f\right)}{\partial t^2}-{\nabla }^{2}u\left(x,t\right)=f\left(x,t\right),& \end{array}$ where x represents space, t is time, m(x) represents the material medium properties, μ(x,t) is the wavefield that may be a scalar or a vector field, ƒ(x,t) represents the sourcing function in space and time, said method comprising the steps of: (a) modifying said sourcing function in said inhomogeneous acoustic wave equation to the mathematical product of a time-independent source attribute defined at every location in space within said material medium and a space-independent source-time function given by $m\left(x\right)\frac{{\partial }^{2}u\left(x,t\right)}{\partial t^2}-{\nabla }^{2}u\left(x,t\right)=A\left(x\right)S\left(t\right),$ where S(t) is the said space-independent source-time function, A(x) represents the said time-independent source attribute, (b) recording acoustic wave data or using wave data obtained through processing in at least one spatial location, and (c) generating an initial source attributes model and an initial material medium properties model, said initial source attributes model optionally comprising only null values at all spatial locations within said material medium, and (d) simulating synthetic wave data using the source properties model and the material medium properties model, and (e) computing an objective function that is a measure of the discrepancy between said simulated wave data and said recorded or pre-processed wave data, and (f) computing the gradients of said objective function with respect to the source attributes model and material medium properties model given by $\begin{array}{cc}{ℊ}_A\left(x\right)={\int }_0^T{u}^{†}\left(x,T-t\right)S\left(t\right)\mathrm{dt},\mathrm{and}{ℊ}_m\left(x\right)=m^{†}\left(x\right){\int }_0^T{u}^{†}\left(x,T-t\right)\frac{{\partial }^{2}u\left(x,t\right)}{\partial t^2}\mathrm{dt},& \end{array}$ where μ.sup.t represents the adjoint wavefield, g.sub.A(x) is the gradient of said objective function with respect to said source attributes model, g.sub.m(x) is the gradient of said objective function with respect to said material medium properties model, and m.sup.t(x)=2/V(x).sup.3 if m(x)=V(x), and m.sup.t(x)=−2/V(x) if m(x)=1/V(x), and m.sup.t(x)=−1 if m(x) 1/V(x).sup.2, and (g) updating said source attributes model and said material medium properties model using
A(x).sub.updated=A(x)−αa.sub.Ag.sub.A(x), and
m(x).sub.updated=m(x)−α.sub.mg.sub.m(x), where A(x).sub.updated is the updated source attributes model, m(x).sub.updated is the updated material medium properties model, α.sub.A and α.sub.m are one-dimensional minimizers of said objective function along vectors A(x)−α.sub.Ag.sub.A(x) and m(x)−α.sub.mg.sub.m(x), respectively, optionally replaced by the Hessian or the pseudo-Hessian, and (h) optionally iterating steps (d)-(g) at least once more, wherein said updated source attributes model and said updated material medium properties model of step (g) arc used in step (d) to simulate synthetic wave data, thereby resulting in a further updated source attributes model and a further updated material medium properties model.

10. A system for efficiently imaging acoustic or processed elastic wave data using the inhomogeneous acoustic wave equation $m\left(x\right)\frac{{\partial }^{2}u\left(x,t\right)}{\partial t^2}-{\nabla }^{2}u\left(x,t\right)=f\left(x,t\right),$ where x represents space, t is time, m(x) represents the material medium properties μ(x,t) is the wavefield that may be a scalar or a vector field, ƒ(x,t) represents the sourcing function in space and time, said system comprising: (a) a processor that modifies said sourcing function in said inhomogeneous acoustic wave equation to the mathematical product of a time-independent source attribute defined at every location in space within said material medium and a space-independent source-time function given by $\begin{array}{cc}m\left(x\right)\frac{{\partial }^{2}u\left(x,t\right)}{\partial t^2}-{\nabla }^{2}u\left(x,t\right)=A\left(x\right)S\left(t\right),& \end{array}$ where S(t) is the said space-independent source-time function, A(x) represents the said time-independent source attribute, (b) a data recording system comprising of one or more computers, data recording hardware, and one or more acoustic wave sensors that record acoustic wave data in at least one spatial location, or a data storage drive containing wave data obtained through processing in at least one spatial location, and (c) wherein the processor further generates an initial source attributes model and an initial material medium properties model, said initial source attributes model optionally comprising only null values at all spatial locations within said material medium, and (d) wherein the processor further simulates synthetic wave data using the source properties model and the material medium properties model, and (e) wherein the processor further computes an objective function that is a measure of the discrepancy between said simulated wave data and said recorded or pre-processed wave data, and (f) wherein the processor further computes the gradients of said objective function with respect to the source attributes model and material medium properties model given by $\begin{array}{cc}\mathrm{ℊA}\left(x\right)={\int }_0^T{u}^{†}\left(x,T-t\right)S\left(t\right)\mathrm{dt},\mathrm{and}{ℊ}_m\left(x\right)=m\left(x\right){\int }_0^T{u}^{†}\left(x,T-t\right)\frac{{\partial }^{2}u\left(x,t\right)}{\partial t^2}dt,& \end{array}$ where μt represents the adjoint wavefield, g.sub.A(x) is the gradient of said objective function with respect to said source attributes model, g.sub.m(x) is the gradient of said objective function with respect to said material medium properties model, and m.sup.t(x)=2/V(x).sup.3 if m(x)=V(x), and m.sup.t(x)=−2/V(x) if m(x)=1/V(x), and m.sup.t(x)=−1 if m(x)=1/V(x).sup.2, and (g) wherein the processor further updates said source attributes model and said material medium properties model using
A(X).sub.updated=A(x)−α.sub.Ag.sub.A(X) and
m(x).sub.updated=m(x)−α.sub.mg.sub.m(x), where A(x).sub.updated is the updated source attributes model, m(x).sub.updated is the updated material medium properties model, α.sub.A and α.sub.m are one-dimensional minimizers of said objective function along vectors A(x)−α.sub.Ag.sub.A(x) and m(x)−α.sub.mg.sub.m(x), respectively, optionally replaced by the Hessian or the pseudo-Hessian, and (h) wherein the processor further optionally iterates steps (d)-(g) at least once more, by using said updated source attributes model and said updated material medium properties model of step (g) in step (d), thereby further updating said source attributes model and further updating said material medium properties model.

11. A method for efficient imaging of acoustic or processed elastic wave data using the frequency-domain inhomogeneous acoustic wave equation
−(m(x)ω.sup.2+Δ.sup.2)U(x,ω)=ƒ(x,ω), where x represents space, ω is frequency, m(x) represents the material medium properties, U(x,ω) is the wavefield that may be a scalar or a vector field, ƒ(x, ω) represents the sourcing function in space and frequency, said method comprising the steps of: (a) modifying said sourcing function in said frequency-domain inhomogeneous acoustic wave equation to the mathematical product of a frequency-independent source attribute defined at every location in space within said material medium and a space-independent source-time function given by
−(m(x)ω.sup.2+Δ.sup.2)U(x,ω)=A(x)S(ω), where S(ω) is the said space-independent source-time function in the frequency domain, A(x) represents the said frequency-independent source attribute, and (b) recording acoustic wave data or using wave data obtained through processing in at least one spatial location, and (c) transforming said recorded acoustic or pre-processed elastic wave data into the frequency domain, and (d) generating an initial source attributes model and an initial material medium properties model, said initial source attributes model optionally comprising only null values at all spatial locations within said material medium, and (e) simulating synthetic wave data using the source properties model and the material medium properties model, and (f) computing an objective function that is a measure of the discrepancy between said simulated wave data and said recorded or pre-processed wave data, and (g) computing the gradients of said objective function with respect to the source attributes model and material medium properties model given by
g.sub.A(x)=Re{U.sup.t,*(x,ω)S(ω)}, and
g.sub.m(x)=−Re{m.sup.t(x)ω.sup.2U.sup.t,*(x,ω)U(x,ω)}, where U.sup.t represents the adjoint wavefield, * represents complex conjugation, Re represents the real part of a complex vector, g.sub.A(x) is the gradient of said objective function with respect to said source attributes model, g.sub.m(x) is the gradient of said objective function with respect to said material medium properties model, and m.sup.t(x)=2/V(x).sup.3 if m(x)=V(x), and m.sup.t(x)=−2/V(x) if m(x)=1/V(x), and m.sup.t(x)=−1 if m(x)=1/V(x).sup.2, and (h) updating said source attributes model and said material medium properties model using
A(X).sub.updated=A(x)−α.sub.Ag.sub.A(x), and
m(X).sub.updated=m(x)−α.sub.mg.sub.m(X), where A(x).sub.updated is the updated source attributes model, m(x).sub.updated is the updated material medium properties model, α.sub.A and α.sub.m are one-dimensional minimizers of said objective function along vectors A(x)−α.sub.Ag.sub.A(x) and m(x)−α.sub.mg.sub.m(x), respectively, optionally replaced by the scaled Hessian or the scaled pseudo-Hessian, and (i) optionally iterating steps (e)-(h) at least once more, wherein said updated source attributes model and said updated material medium properties model of step (h) are used in step (e) to simulate synthetic wave data, thereby resulting in a further updated source attributes model and a further updated material medium properties model.

12. A system for efficiently imaging acoustic or processed elastic wave data using the frequency-domain inhomogeneous acoustic wave equation
−(m(x)ω.sup.2+Δ.sup.2)U(x,ω)=ƒ(x,ω), where x represents space, w is frequency, m(x) represents the material medium properties, U(x,ω) is the wavefield that may be a scalar or a vector field, ƒ(x,ω) represents the sourcing function in space and frequency, said system comprising: (a) a processor that modifies said sourcing function in said frequency-domain inhomogeneous acoustic wave equation to the mathematical product of a frequency-independent source attribute defined at every location in space within said material medium and a space-independent source-time function given by
−(m(x)ω.sup.2+Δ.sup.2)U(x,ω)=A(x)S(ω), where S(ω) is the said space-independent source-time function in the frequency domain, A(x) represents the said frequency-independent source attribute, and (b) a data recording system comprising of one or more computers, data recording hardware, and one or more acoustic wave sensors that record acoustic wave data in at least one spatial location, or a data storage drive containing wave data obtained through processing in at least one spatial location, and (c) wherein the processor further transforms said recorded acoustic or pre-processed elastic wave data into the frequency domain, and (d) wherein the processor further generates an initial source attributes model and an initial material medium properties model, said initial source attributes model optionally comprising only null values at all spatial locations within said material medium, and (e) wherein the processor further simulates synthetic wave data using the source properties model and the material medium properties model, and (f) wherein the processor further computes an objective function that is a measure of the discrepancy between said simulated wave data and said recorded or pre-processed wave data, and (g) wherein the processor further computes the gradients of said objective function with respect to the source attributes model and material medium properties model given by
g.sub.A(x)=Re{U.sup.t,*(x,ω)S(ω)}, and
g.sub.m(x)=−Re{m.sup.t(x)ω.sup.2U.sup.t,*(x,ω)U(x,ω)}, where U.sup.t represents the adjoint wavefield, * represents complex conjugation, Re represents the real part of a complex vector, g.sub.A(x) is the gradient of said objective function with respect to said source attributes model, g.sub.m(x) is the gradient of said objective function with respect to said material medium properties model, and m.sup.t(x)=2/V(x).sup.3 if m(x)=V(x), and m.sup.t(x)=−2/V(x) if m(x)=1/V(x), and m.sup.t(x)=−1 if m(x)=1/V(x).sup.2, and (h) wherein the processor further updates said source attributes model and said material medium properties model using
A(x).sub.updated=A(x)−α.sub.Ag.sub.A(x), and
m(x).sub.updated=m(x)−α.sub.mg.sub.m(x), where A.sub.(x)updated is the updated source attributes model, m(x).sub.updated is the updated material medium properties model, α.sub.A and α.sub.m are one-dimensional minimizers of said objective function along vectors A(x)−α.sub.A g.sub.A(x) and m(x)−α.sub.m g.sub.m(x), respectively, optionally replaced by the scaled Hessian or the scaled pseudo-Hessian, and (i) wherein the processor further optionally iterates steps (e)-(h) at least once more, by using said updated source attributes model and said updated material medium properties model of step (h) in step (e), thereby further updating said source attributes model and further updating said material medium properties model.

Description

DRAWINGS

Figures

(1) For a more complete understanding of the invention, reference is made to the following description and accompanying drawings, in which:

(2) FIG. 1 illustrates an example embodiment of an arrangement of surface and subsurface sensors positioned to monitor incidence of P-waves and S-waves recordable within said substrate;

(3) FIG. 2 is a flow chart illustrating an example embodiment of steps of the present method through sequential inversion whereby hypocenter imaging is updateable subsequent each updated velocity model within each iteration;

(4) FIG. 3 is a flow chart illustrating an example embodiment of steps of the present method through simultaneous inversion whereby hypocenter imaging is updateable concurrent each updated velocity model within each iteration;

(5) FIG. 4 shows the compressional-wave velocity model V.sub.P used in generating the measured compressional-wave data;

(6) FIG. 5 shows the shear-wave velocity model V.sub.S used in generating the measured shear-wave data;

(7) FIG. 6 shows the velocity combination V=V.sub.PV.sub.S/(V.sub.P−V.sub.S) of said velocity models in FIG. 4 and FIG. 5;

(8) FIG. 7 shows the locations of 80 sources used in generating the elastic wave data;

(9) FIG. 8 shows a time-gated portion of the compressional-wave data recorded at the receivers on the top boundary of the model. The whole record spans over 10 minutes and is not shows here;

(10) FIG. 9 shows a time-gated portion of the shear-wave data recorded at the receivers on the top boundary of the model. The complete record spans over 10 minutes and is not shows here;

(11) FIG. 10 shows the waveforms obtained by a moving-window deconvolution of the shear- and compressional-waveforms in FIG. 9 and FIG. 8 at each receiver. The energy about zero-lag has been suppressed;

(12) FIG. 11 shows the starting velocity model for waveform inversion in the test example;

(13) FIG. 12 shows the velocity model obtained using waveform inversion; and

(14) FIG. 13 shows the source amplitude model obtained using waveform inversion.

DETAILED DESCRIPTION

(15) With reference now to the drawings, and in particular FIG. 1 through FIG. 3 thereof, example of the instant method for expediently processing seismic data for subsurface mapping with waveform inversion employing the principles and concepts of the present method for expediently processing seismic data for subsurface mapping with waveform inversion and generally designated by the reference number 100 will be described.

(16) Referring to FIG. 1 through FIG. 3 an example embodiment of the method for expediently processing seismic data for subsurface mapping with waveform inversion 100 is illustrated.

(17) The present method for expediently processing seismic data for subsurface mapping with waveform inversion 100 includes acquiring seismological waveform data 14 in situ from an arrangement of sensors, 11a and 11b, over a determinable time period. The present method is usable upon waveform data 14 acquired by single component and multicomponent sensors, as case may be. Once waveform data 14 is produced, waveform data 14 is either auto-correlated, cross-correlated, or deconvolved or decomposed into shear- and compressional-waveforms which are thereafter cross-correlated or deconvolved over said time period to generate a waveform comprising a phase difference of each compressional-wave and shear-wave incident sensed at each sensor location, and a relevant waveform data set 17 for each sensor location is thereby generated. Waveform data 14 continuously collected over an extended time period is thus reduced to relevant data sets 17 readable in seconds.

(18) Waveform inversion of each relevant data set 17 determines subsurface velocity distribution and source properties distribution across the substrate and a subsurface map is thereby generable. The inverted velocity is inversely proportional to the difference between the compressional-wave and shear-wave slowness of the substrate.

(19) Monitoring velocity distribution across the substrate during subsurface activity therein, by continuously generating relevant data sets at each sensor location, thence provides updatable mapping across the scope of said subsurface activity.

(20) FIG. 1 illustrates a diagrammatic longitudinal section of a substrate beneath surface 5 to depth 6 having an arrangement of surface sensors 11a and subsurface sensors 11b disposed appropriately to monitor seismic activity therein. Sensors 11a and 11b are connected with a data recording unit 12 wherein waveform data are storable and communicable to remote users.

(21) Compressional-wave 10 and shear-wave 9, generated at hypocenter 8, are transmitted through subsurface geologic strata 1, 2, 3, and 4, and thence recorded when incident each of sensors 11a and 11b. Incidence of compressional-wave 10 and shear-wave 9 are enumerated as waveform data 14 for cross-correlation or deconvolution 15b,c wherein a waveform is generative comprising the phase difference between said compressional-wave and shear-wave incidents whereby waveform data 14 is reduced to relevant data sets 17. Waveform data 14 is thereby compressed for expedited processing to generate subsurface high resolution mapping by waveform inversion.

(22) Referring to FIG. 2 and FIG. 3, waveform data 14 is auto-correlated, cross-correlated, or deconvolved using a sliding-window operation 15b,c and energy is suppressed at zero-lag 16 whereby sequential data are processed to delimit relevant data sets 17 from raw waveform data 14. Initial hypocenter image 24 is thus generated by reverse time migration of the relevant data using the initial velocity model 23. Independently, source-time function is estimated 18 from raw waveform data 14 and auto-correlated (or deconvolved) 20 to produce a modified source-time function 22. If shear-wave and compressional-wave decomposition 15a is performed 15c using a polarization filter, then the envelope data may be generated from the decomposed shear- and compressional-waveforms (instead of prior to the shear-wave and compressional-wave decomposition 15a).

(23) Relevant data sets 17 are then processed in conjunction with the modified source-time function 22 to compute the gradient of objective function 25. Convergence 26 of synthetic data with relevant data sets 17 updates hypocenter imaging 29 and velocity distribution 27 assessed the substrate, whereby updatable high resolution mapping of the substrate is enabled. At least one iteration is performed.

(24) When convergence 26 is achieved, the process is ended. When convergence 26 is not achieved, the velocity model is updated 27 and the hypocenter image is updated 29 to compute the gradient of objective function 25 until convergence is met.

(25) As shown in FIG. 2, updating of the velocity model 27 and the hypocenter image 29 may be performed sequentially, or, as shown in FIG. 3, concurrently, as desired. Update of the velocity model 27 results in an updated velocity model 28 and update of the hypocenter image 29 results in an updated hypocenter image 30.

(26) Thus a high resolution map is generable and updatable from continuously sourced waveform data 14 acquired at an arrangement of sensors 11a, 11b disposed appropriate to record seismic activity within a substrate, and seismic data affected during activity within said substrate is thereby usable to generate relevant data sets to which waveform inversion is specifically applied to rapidly and efficiently determine velocity distributions and hypocenter image within said substrate.

(27) Many examples of equations and mathematical expressions have been provided in this disclosure. But those with skill in the art will appreciate that variations of these expressions and equations, and related expressions and equations that can be derived from the example equations and expressions provided herein may also be successfully used to perform the methods, techniques, and work-flows related to the embodiments disclosed herein.

(28) While illustrative embodiments of the invention have been illustrated and described, it will be appreciated that various changes can be made therein without departing from the scope of the invention. For example, reduction of recorded data might not be necessary but waveform inversion, as set forth in the present invention, may be used for obtaining physical properties of a material medium and a source properties distribution. On other instances, only data reduction, as set forth in the present invention, may be used to expedite other processing and inversion.

EXAMPLE

(29) FIG. 4 through FIG. 12 present a synthetic example of reducing the volume of microseismic data and inverting the resulting waveforms for obtaining said velocity combination V and source amplitude distribution using the method described in this invention.

(30) FIG. 4 is the true compressional-wave velocity model V.sub.P that is used to generate the compressional-wave data. FIG. 5 is the true shear-wave velocity model V.sub.S that is used to generate the shear-wave data. The spatial extent of the model is 3.5 km in the x-direction and 2.37 km in the z-direction.

(31) FIG. 6 shows the true compressional- and shear-wave velocity model combination V=V.sub.PV.sub.S/(V.sub.P−V.sub.S) of the velocity models in FIG. 4 and FIG. 5 that the reduced waveforms are sensitive to.

(32) FIG. 7 shows the spatial locations of 80 sources used in generating the compressional- and shear wave data. The amplitudes of all the sources are the same and equal to unity. Receivers are located on all four boundaries of the model at an interval of 20 m.

(33) Time-gated portions of the generated compressional- and shear-wave data are shown in FIG. 8 and FIG. 9, respectively. The complete data record spans about 10 minutes in this test example and is not shown here.

(34) FIG. 10 shows the reduced waveform data obtained from the compressional- and shear-wave data in FIG. 8 and FIG. 9 using deconvolution. A time gate size of 2 sec was used to generate the waveforms in FIG. 10. Also, the energy about zero-lag has been suppressed. Note that the original recorded data spanning about 10 minutes has been reduced to only 2 seconds which will substantially reduce the computational time involved in waveform inversion.

(35) FIG. 11 shows a homogeneous starting velocity model used in waveform inversion. Because of weak velocity heterogeneities, a homogeneous starting velocity model was sufficient for this test example. Since waveform inversion is a highly nonlinear problem, a good starting velocity model may be necessary is cases with complicated velocity models.

(36) Note that no prior knowledge of source amplitude distribution was used in starting the waveform inversion.

(37) FIG. 12 shows the velocity model obtained using waveform inversion of the waveform data using the method described in this invention in FIG. 10. The source amplitude distribution obtained during the same inversion is shown in FIG. 13. The Limited-memory Broyden-Fletcher-Goldfarb-Shanno optimization algorithm [19] was used in solving the inversion problem. Comparisons of FIG. 12 with FIG. 6 and FIG. 13 with FIG. 7 demonstrate that the waveform inversion method described in this invention was successful at obtaining a reasonable model.

(38) The foregoing application is directed to particular embodiments of the present invention for the purpose of illustrating it. It will be apparent, however, to one skilled in the art, that many modifications and variations to the embodiments described herein are possible. All such modifications and variations are intended to be within the scope of the present invention, as defined in the appended claims. Persons skilled in the art will readily recognize that in preferred embodiments of the invention, at least some of the steps in the present inventive method are performed on a computer, i.e. the invention is computer implemented. In such cases, the resulting updated physical properties model and the updated source properties model may either be downloaded, displayed, or saved to computer storage.

REFERENCES

(39) [1] C. A. Tosaya, “Acoustical properties of clay-bearing rocks,” Ph.D. dissertation, Stanford University, 1982. [2] J. T. Rutledge and W. S. Phillips, “Hydraulic stimulation of natural fractures as revealed by induced microearthquakes, carthage cotton valley gas field, east texas,” Geophysics, vol. 68, no. 2, pp. 441-452, 2003. [Online]. Available: http://dx.doi.org/10.1190/1.1567214 [3] K. Chambers, S. Brandsberg-Dahl, J. Kendall, and J. Rueda, “Testing the the ability of surface arrays to locate microseismicity,” SEG Technical Program Expanded Abstracts, pp. 1436-1440, 2008. [Online]. Available: http://library.seg.org/doi/abs/10.1190/1.3059380 [4] M. Hall and J. E. Kilpatrick, “Surface microseismic monitoring of slick-water and nitrogen fracture stimulations, arkoma basin, oklahoma,” SEG Technical Program Expanded Abstracts, pp. 1562-1565, 2009. [Online]. Available: http://library.seg.org/doi/abs/10.1190/1.3255147 [5] B. Artman, I. Podladtchikov, and B. Witten, “Source location using time-reverse imaging,” Geophysical Prospecting, vol. 58, no. 5, pp. 861-873, 2010. [Online]. Available: http://dx.doi.org/10.1111/j.1365-2478.2010.00911.x [6] S. C. Maxwell, J. Rutledge, R. Jones, and M. Fehler, “Petroleum reservoir characterization using downhole microseismic monitoring,” GEO-PHYSICS, vol. 75, no. 5, pp. 75A129-75A137, 2010. [Online]. Available: http://dx.doi.org/10.1190/1.3477966 [7] J. P. Verdon, J.-M. Kendall, and S. C. Maxwell, “A comparison of passive seismic monitoring of fracture stimulation from water and co2 injection,” GEOPHYSICS, vol. 75, no. 3, pp. MA1-MA7, 2010. [Online]. Available: http://dx.doi.org/10.1190/1.3377789 [8] J. Virieux and S. Operto, “An overview of full-waveform inversion in exploration geophysics,” Geophysics, vol. 74, no. 6, pp. WCC1-WCC26, 2009. [Online]. Available: http://link.aip.org/link/?GPY/74/WCC1/1 [9] P. Lailly, “The seismic inverse problem as a sequence of before stack migrations,” in Abstracts. Society of Industrial and Applied Mathematics, 1984, pp. 206-220. [10] A. Tarantola, “Inversion of seismic reflection data in the acoustic approximation,” Geophysics, vol. 49, no. 8, pp. 1259-1266, 1984. [Online]. Available: http://link.aip.org/link/?GPY/49/1259/1 [11] R. Kamei and D. Lumley, “Passive seismic imaging and velocity inversion using full wavefield methods,” SEG Technical Pro-gram Expanded Abstracts, pp. 2273-2277, 2014. [Online]. Available: http://library.seg.org/doi/abs/10.1190/segam2014-0948.1 [12] J. Tromp, C. Tape, and Q. Liu, “Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels,” Geophysical Journal International, vol. 160, no. 1, pp. 195-216, 2005. [Online]. Available: http://dx.doi.org/10.1111/j.1365-246X.2004.02453.x [13] Q. Liu and J. Tromp, “Finite-frequency kernels based on adjoint methods,” Bulletin of the Seismological Society of America, vol. 96, no. 6, pp. 2383-2397, December 2006. [14] Y. Kim, Q. Liu, and J. Tromp, “Adjoint centroid-moment tensor inversions,” Geophysical Journal International, vol. 186, no. 1, pp. 264-278, 2011. [Online]. Available: http://dx.doi.org/10.1111/j.1365-246X.2011.05027.x [15] C. Tape, Q. Liu, A. Maggi, and J. Tromp, “Adjoint tomography of the southern california crust,” Science, vol. 325, no. 5943, pp. 988-992, 2009. [Online]. Available: http://www.sciencemag.org/content/325/5943/988.abstract [16] A. Askan, “Full waveform inversion for seismic velocity and anelastic losses in heterogeneous structures,” Ph.D. dissertation, Carnegie Mellon University, 2006. [17] R.-E. Plessix, “A review of the adjoint-state method for computing the gradient of a functional with geophysical applications,” Geophysical Journal International, vol. 167, no. 2, pp. 495-503, 2006. [Online]. Available: http://dx.doi.org/10.1111/j.1365-246X.2006.02978.x [18] O. Jarillo Michel and I. Tsvankin, “Gradient calculation for waveform inversion of microseismoic data in vti media,” Journal of Seismic Exploration, vol. 23, no. 3, pp. 201-207, 2014. [19] J. Nocedal and S. J. Wright, Numerical Optimization, 2nd ed. New York: Springer, 2006.