Robust Hydraulic Fracture Monitoring (HFM) of Multiple Time Overlapping Events Using a Generalized Discrete Radon Transform

R OBUST HYDRA UL IC FRA CT URE MONITORING (HFM) OF MUL TIPLE TIME O VERL APPING EVENTS USING A GENERALIZED DISCRETE RADON TRANSFORM Gr e gory Ely and Shuchin Aer on T ufts University , Dept. of ECE, Medford MA ABSTRA C T In this work we propose a nov el algorithm for multiple-e vent localization for Hydraulic Fracture Monitoring (HFM) through the exploitation of the sparsity of the observ ed seismic signal when rep- resented in a basis co nsisting of space time propagators. W e prov ide explicit co nstruction of these propagators using a forw ard model for wa ve propagation which depends non-linearly on the problem pa- rameters - the unknown source location and mechanism of fracture, time and extent of even t, and the locations of t he receiv ers. Under fairly gen eral assumptions and an appropriate discretization of these parameters we first b uild an ov er-complete d ictionary of generalized Radon propagators and assume that the data is well represented as a linear superpo sition of t hese propagators. Exploiting this structure we propos e sparsity penalized algorithms and workflo w for super- resolution extraction of time overlapping multiple seismic ev ents from single well data. 1. INTRODUCTION Accurate seismic hydraulic fracturing monitoring (HFM) can miti- gate many of the en vironmental impacts by providing a clear real- time image of where the fractures are occurring outside of the shale and ho w ef fi ciently they a re fo rmed within the gas deposit. Although simple in principle, real time monitoring of hydraulic fracturing is extremely difficult to perform successfully due to high noise lev- els generated by the pumping equipment, anisotropic propagation of seismic wav es through shale, and the multi-l ayered stratigraphy leading to comple x seismic ray propagation, [1, 2, 3]. In addition the complex ity of the source mechanism affects the wave polariza- tion at the 3-axis seismometers, [4 ] introducing extra parameters in the system. 2. PROBLEM SET -UP The physical set-up is sho wn in Figure 1. A typical seismic array of (say) L three-axis seismometers sample the instantaneous displace- ment across the three axes at a sampling rate between 1-8 kHz. The j th detector in the array records a trace, see Figure 2, along the three- axes that wil l be a combination of m seismic even ts observed i n noise which can be ef fectiv ely modeled as Additiv e White Gaussian Noise (A WGN), see [5]. Therefore we hav e, Y j ( t ) = m X i =1 (( S ij ( t ) + N ( t )) (1) where, Y j ( t ) =   Y j x ( t ) Y j y ( t ) Y j z ( t )   S ij ( t ) =   S ij x ( t ) S ij y ( t ) S ij z ( t )   (2) In the most gen eral case of anisotropic formation S ij ( t ) is a co mbi- θ ϕ r fracture horizontal well θ min θ max r min r max ϕ min ϕ max St r i k e a n g l e Di p An g l e r a n g e 3 - a x is se n so r b o re h o l e w i t h se i em o m e t er s Fig. 1 . This figure sho ws t he typical hydraulic fracturing setup in which a seismic array detects the emitted signal from a fracture oc- curring within the search vo lume. nation of compressional ( ρ ) wav e component and shear componen t which further con sists of vertical shea r component, s v , and horizon - tal shear component, s h , as incident on the detector along the three axes, i.e., S ij ( t ) = S ij ρ ( t ) + S ij s h ( t ) + S ij s v ( t ) (3) For a given a wav e type such as ρ , the signal at the seismometer induced by the wav e, S ij ρ ( t ) , is completely determined by the sig- nal waveform of the compressional wave at the source, s iρ ( t ) , the arriv al time f or the compressional wa ve of the i th e vent at the j th detector , τ ρ ij , the compressional polarization unit vector , P ij ρ , rep- resenting the direction o f particle mov ement at the detector . Usually the si gnal waveform s iρ ( t ) i s time compact due to finite duration of the seismic ev ent. Therefore the signal at the seismometer is the signal at the source, s iρ ( t ) , delayed by a time, τ ρ ij , and projected onto the three axes by t he polarization unit vector , P ij ρ . Therefore, mathematically we can write, S ij ρ ( t ) = δ ( t − τ ρ ij ) ∗ ( P ij ρ s ij ρ ( t )) = P ij ρ s ij ρ ( t − τ ρ ij ) (4) The time and location of the event, say { t i , l i } , and velocity profile of the st ratigraphy and the seismic source mechanism completely determine the arriv al time of the wave and the polarization vector at each of t he detector . This is expressed i n t erms of a forward model function, { P ij ρ , τ ρ ij } = f j ( t i , θ i , φ i , r i ) | {z } l i ) (5) where we have characterized the location of the even t at l i in terms of the strike, dip and range relativ e to an absolute co-ordinate sys- tem. Giv en t hat t he seismic ev ent at a one location will not affect the signal emitted from a seismic ev ent at another location, we can assume that the observed data Y j ( t ) at the j − th r ecei ver resulting from e vents that emit all three wav e t ypes can be vie wed as a lin- ear combina tion of signals observed in A WGN. In this work we will focus on using the compressional wav es for HFM as these waves typically hav e higher amplitude than the shear waves. Howe ver , our approach can be easily extended to handle the processin g of shear wa ves. Using equation (5), we can rewrite the Equation (1) for the observ ed trace due to compressional wav es, Y j ( t ) as follo ws, Y j ( t ) = m X i =1 P ij ρ s i ρ ( t − τ ρ ij ) + N ( t ) (6) Note t hat here we focus on isotropic fracture mechanism where the moment tensor is kno wn sav e for the scaling f actor . Nev ertheless our approach can be extende d to jointly estimate the time, location and moment tensor for other mechan isms such as double couple, [2]. 3. SP ARSITY PENALIZED FRAMEWORK FOR HFM T o ill ustrate the sparsity of the HFM signal, we begin by defining a propagator , Γ i ρ = { Γ i j ρ ( t ) } L j =1 , which correspon ds to noiseless data at the receiv ers as excited by a hypothetical seismic ev ent i (say) at location l i and time t i with an impulse source signal, i.e., Γ i j ρ ( t ) = δ ( t − τ ρ ij ) ∗ P ij ρ (7) W ith this notion of a propagator we construct a basis for a general- ized Radon transform as follo ws. W e discre tize the search volume in space and time with locations within the search v olume say V in- dex ed from 1 to n V and time of the ev ents in this volume indexed from 1 to n t V . Then for each of the discretized locations and time we construct a propagator using Equation (7). For sake of expo si- tion we vectorize the space time propagators by stacking the t ri-axial time traces Γ i j ρ ( t ) for each receiv er j as column v ectors obtaining a 3 × T × L vector for each propagator . W e denote this vector by Γ i ρ (:) . (NO TE - here (:) denotes the ve ctorization operation simi- lar to the MA TL AB operator (:) ). Running ov er the l ocations and time indices o ver the e vents, we collect the vectorized propag ators as columns of a matrix denoted by Φ ρ . The data is also vectorized to a vector , Y ρ (:) which is a 3 × T × L ve ctor in a simil ar manner as t he propagators. Then a generalized discrete Radon t ransform of the data is gi ven by , R ( Y ) = Φ † ρ Y ρ (:) (8) where Φ ρ = h Γ 1 ρ (:) , Γ 2 ρ (:) , . . . , Γ i ρ (:) , . . . , Γ n V × n t V ρ (:) i (9) and † denotes t he conjugate transpose. This linear operator , Φ † ρ , transforms the vectorized trace data, Y (:) to a function of source location and the time at which an eve nt occurs. It is essentially sim- ilar to generating a slant-stack (or time-mov eout) map of the data, albeit along t he non-linear trajectories. Note that because the noise is spatially incoherent the transform wi ll spread the noise energy e venly across all of the possible propagators. From t his transform we can see that the trace data becomes sparse when viewed in the in the generalized Radon domain (fi gure 2 ). Ho wev er because the dictionary is over complete, we see a number of sidelobes for just a single seismic even t, making it difficult to determine the l ocation of the ev ent in heav y noise. As stated in the problem statement we assume that the signal at the source is ti me compact. Therefore we assume that the signal s i for an even t i is supported on a set T i = [ t min , t max ] ⊂ [0 , T ] for some total time T of microseismic stimulation. In this case the source signals s i ρ ( t ) can be viewed as weighted sum of sev eral delta functions, s i ρ ( t ) = X t k ∈T i ρ η i ρ t k δ ( t − t k ) (10) where the vector η t k determine the shape of t he source wa veform. Here we would like to point our that in contrast to methods that as- sume knowledge of the signal wave form, [5], here we have not as- sumed any kno wledge of the wav e shape. Howe ver , we assume that the search volume contains t he set T i ρ for the ev ent under consider- ation. s i ρ ( t ) = X t k ∈ [0 ,T ] η i ρ t k δ ( t − t k ) : η i t k = 0 ∀ t k / ∈ T i ρ (11) Giv en our definition of a propagator , assuming no geometrical spreading, we can define the arriving observed vectorized traces Y ρ (:) du e t o an e vent as a sum of propagators with the same source location and v arying ev ent times. Therefore we can wri te, Y ρ (:) = Φ ρ X (:) + N (12) where X ( :) ∈ R n V · n t v is the coefficient vector (with elements η t k i ) corresponding to the spatio-temporal volume of possible e vents. If we were to reshape the true vector X (:) into a two dimensional ma- trix X the dictionary weights corresponding to a single ev ent would be sparse along the location dimension while compact along the tem- poral dimension, see Fi gure 2 . In this way the signal exh ibits sparsity across the ro ws of X and is thus said to be simultaneou sly sparse. 4. SP ARSITY PENALIZED ALGORITHM FOR EVENT LOCALIZA TION Moti vated by the ro w-sparse structure of the data, we use the al- gorithm presented in [6] for a high spatio-temporal resolution map- ping of the micro-seismic ev ents. The algorithm corresponds to the follo wing mathematical optimization problem also kno wn as group sparse penalization, [7] in the literature. ˆ X = ar gmin X || Y (:) − ΦX (:) || 2 + λ || X || 1 , 2 (13) where λ is a sparse tuning factor that controls the group sparseness of X versus the signal error . Where group or row spareness i s the measure of n umber o f non-z ero ro ws of X . W e can in duce this mini- mizing the t he ℓ 1 norm of the eu clidian norm o f the ro ws (same time indices) of X . The parameter λ is chosen depen ding on the noise lev el and the anticipated number of ev ents. Fo r the application at hand we use the method proposed in [8][9] for selection of a good regularization parameter . In this context our method is different tha n that used in [10 ] in that we are exploiting the row-sparsity of the generalized discrete radon transform. 4.1. Algorithm wor kflow The geometry of the linear seismic array l ends itself to dividing the task of localization into two subtask s, viz., • Localization in the sagittal plane (estimate of θ and r ); and • Determination of the strike angle φ . In order to determine the location o f the seismic ev ent in t he sagittal plane the two horizontal axes (x and y) of t he traces are combined 1500 1600 1700 1800 0 2 4 6 8 10 12 p−p = 8.5e−1 SNR 15 SNR 15 Radon 5 10 15 20 20 40 60 80 100 120 1500 1600 1700 1800 0 2 4 6 8 10 12 p−p = 1.7e+0 SNR 5 SNR 5 Radon 5 10 15 20 20 40 60 80 100 120 Fig. 2 . This figure sho ws the traces f or a synthetic ρ -wav e across fo r tw o SNRs: 15 a nd 5. Also sho wn are the (gen eralized) Rado n transforms. Note the signal sparsity in the Radon domain which is concentrated at a single position index an d among a few time indices. 1400 1500 1600 1700 0 2 4 6 8 10 12 p−p = 9.4e−1 x 1400 1500 1600 1700 0 2 4 6 8 10 12 p−p = 5.3e−1 y 1400 1500 1600 1700 0 2 4 6 8 10 12 p−p = 7.0e−1 z Fig. 3 . This figure sho ws the x,y , and z traces for all 6 synthetic even ts in the presence of noise. into a single quantity Y h ( t ) allowing the data to be observed in the Y h ( t ) an d Y z ( t ) d omain, Y ′ ( t ) , where Y j h ( t ) = p Y j x ( t ) 2 + Y j y ( t ) 2 Y ′ j ( t ) =  Y j h ( t ) Y j z ( t )  (14) Furthermore, t he transformation is applied to each dictionary ele- ments and the modified dictionary Φ hz is constructed. By con- structing the diction ary in this domain bo th the n umber of dictionary elements are reduced by the number of possible strike angles and decreases the size of each dictionary entry by the number of ti me samples. Ho weve r , because the amplitude difference across the x and y comp onents are e liminated in this procedure, all e vents hav ing the same str ike will have the same amplitude in the radon domain. The ℓ 1 , 2 minimization that was described in Equation 13, can no w be applied to the t races i n the Y ′ ( t ) space. The minimization now takes the form of 15 in which Φ hz is the compressed dictionary and X hz is the d ictionary weights and results in the estimation of signals time and sagittal location support, ˆ X hz . ˆ X hz = ar gmin X hz || Y hz (:) − Φ hz X hz (:) || 2 + λ || X hz || 1 , 2 (15) Once the minimization has been performed the r ange and depth of the location can then be determined by thresholding the dictionary coef ficients. T o locate the ev ents or a si ngle ev ent we take the ℓ 2 norm of eac h ro w across time to generate a co lumn ve ctor whose en- tries represen t the total signa l energy at a gi ven location . The largest entries of these vectors correspond to location of the seismic ev ent. Once a set of seismic ev ents hav e been located in the sagittal plane the algorithm has i dentified the ev ents’ str ike, range, and time with some uncertainty , the ev ent’ s true location i s effecti vely constrained to sev eral semi-toruses of constant range and dip as the only un- kno wn parameter is strike. Because of the possible locations of the e vent hav e been constrained to a small subset of t he entire search volume , the computational load of performing the minimization on the full three axes traces has been significantly decreased. In order to estimate the str ike of a seismic ev ent a subset of the original dic- tionary , Φ xy z , consisting of propagators in the full axes setup only along the semi-toruses. Once this small dictionary is constructed a ℓ (1 , 2) minimization is then applied to the fu ll three dimensional data to estimate the signal’ s strike support ˆ X xy z . ˆ X xy z = ar gmin X xyz || Y xy z (:) − Φ xy z X xy z (:) || 2 + λ || X xy z || 1 , 2 As with the sagittal l ocalization, we again t ake the ℓ 2 norm of each ro w of the dictionary coefficients , ˆ X xy z , across time t o generate a column vector . Because the number of ev ents has been determined in sagittal localization, the strike of each eve nt can be determined by pickin g the correspond ing number of largest peaks of the column vector . 5. PERFORMANCE EV ALU A TION ON SYNTHETIC D A T A T o illustrate the performance of the proposed method 6 seismic e vents were generated in short succession after with a temporal spacing of 5 milliseconds, see Figure 3. F urthermore, the seismic e vents follow a near linear path occurring roughly 50 meters apart. In this simulation a search volume with a depth from 750 to 1050 meters, a range f rom 750 to 1050 meters, and set of strike angles from 0 to 30 de grees. A spatial resolution of 50 meters and .1 degree s was used. Each of the seismic ev ents generated a 150Hz compressional ri cker wav elet with the same amplitude across all Radon transform location index time index 5 10 15 20 10 20 30 40 50 60 dictionary coeffcients time index 5 10 15 20 10 20 30 40 50 60 0 0.5 0 10 20 30 40 50 60 L 2 row norms Radon Transform location index time index 5 10 15 20 20 40 60 80 100 120 140 160 180 dictionary coeffcients time index 5 10 15 20 20 40 60 80 100 120 140 160 180 0 0.05 0.1 0 50 100 150 L 2 row norms Fig. 4 . Left: This figure shows the resulting dictionary coefficients of the sparse penalized mi nimization for the sagittal localization and the radon transform of the observed data. Right: This figure shows t he resulting dictionary coefficients of the sparse penalized minimization of the strike minimization and the strike radon transformed data. T he ev ents’ locations are easily identifiable after applying the minimization. The 6 distinct horizontal bands in this image correspond to the dictionary elements from the 6 dif ferent sagittal locations. e vents. Although t hese traces do not overlap along some of the detectors, at many detectors the ev ents fall on top of each other , making it difficult to determine the arriv al times of each ev ent and corrupting the relativ ely amplitude across the x and y axes (fi gure 3). As in the multiple even ts simulation the same ℓ (1 , 2) minimiza- tion and thresholding operations were applied to the data. Figure 4 shows the resulting radon transform and dictionary coefficients. After applying the minimization it easy to determine the strike and sagittal location of each of the 6 e vents. In order to characterize the performance of our algorithm in t he presence of noise, the simulation w as repeated 25 times with a single seismic ev ent occurring i n the mi ddle of the search volume with the same detector configuration bu t with varyin g S NR’ s from 1 to 15. In addition, t he sagittal resolution was increases to 5 meters. T he results of the application of the algorithm to the varying SNR cases are shown in fi gure 5. T he algorithm is able to accurately localize seismic ev ents in the sagittal to plane to within 10 meters for signal to noise ratios abov e 8. Howe ver , below a signal to noise ratio of 5 t he algorithm beco mes unstable the algorithm consistently estimates the e vent as occ urring on the edge of the search volum e 50 meters aw ay from the true e vent’ s location, driving t he appare nt track accurac y to 50 meters second. 6. CONCLUSION & FUTURE WORK W e demonstrated a nov el method based on sparsity pen alized recon- struction methods for accurate HF M. Our method assumes minimal kno wl edge of the waveform signatures and also exploits the tempo- ral information in the signal. Currently our processing works on the compressional and shear wav es separately . In future w e will extend this framewo rk to jointly process these arri vals and exploit the de- pendencies in polarization to impro ve the accuracy . 7. A CKNO WLEDGEMENTS The authors will like to ackno wledge many helpful discussions and suggestions from Dr . Sandip Bose at Schlumberger -Doll Research, Cambridge, MA. 8. REFERENCES [1] Leo E isner and Peter M. Duncan, “Uncertainties i n passi ve seismic monitoring, ” The L eading Edge 28 , vol. 28, pp. 648– 655, 2009. [2] K eiiti Aki and Paul G. Richards, Quantitative Seismology , 2nd Edition , Univ ersity Science Books, 2002. 5 10 15 3 4 5 6 7 8 9 10 saggital error SNR Accuarcy std (meters) range depth 5 10 15 2 4 6 8 10 12 Strike and time error SNR Strike std (degrees) 5 10 15 2 3 4 5 6 7 time std (degrees) strike time error Fig. 5 . Left figure sho ws the accuracy in sagittal position estimation with SNR averaged o ver 25 synthetic traces at each S NR. Right fig- ure sho ws the accuracy in strike and time estimation o f as a fun ction of SNR. [3] Peter M. Shearer , Intr oduction to Seismology , Cambridge Uni- versity Press, 200 9. [4] C. H. Chapman and W . S. L eaney , “ A ne w moment-tensor decomposition for seismic events in anisotropic media, ” Geo- physical J ournal International , vol. 188, no. 1, pp. 343–370, 2012. [5] Qiuhua Liu, S. Bose, H.-P . V alero, R. G. Shenoy , and A. Ounad- jela, “Detecting small amplitude signal and transit times in high noise: Application to hydraulic fracture monitoring, ” in IEEE Geoscience and Remote Sensing Symposium , 2009. [6] J. A. T ropp, A. C. Gi lbert, and M. J. Strauss, “ Algorithms for simultaneous sparse approximation. part II: Con vex rel ax- ation, ” Signal Pr ocessing , special issue on Sparse appr oxi- mations in signal and image pr ocessing , vol. 86, pp. 572–588, April 2006. [7] A. Majumdar and R .K. W ard, “Fast group sparse classifica- tion, ” Electrical and Computer Engineering, Canadian J our- nal of , vo l. 34, no. 4, pp. 136 –144, fall 2009. [8] S. Aeron, S. Bose, H- P V alero, and V . Sali grama, “Broadband dispersion extraction using simultaneous sparse penalization, ” IEEE T ransactions on Signal Pr ocessing , 2011, T o appear , av ailable on IEE E Explore. [9] S. Bose, S. Aeron, and H-P V alero, “Joint multi -mode disper- sion extraction in fourier and space time domains, ” Geophysi- cal J ournal International , submitted. [10] I. V . Rodriguez, M. Sacchi, and Y . J. Gu, “Simultaneo us re- cov ery of origin time, h ypocentre location and se ismic moment tensor using sparse representation theory , ” Geophysical Jour - nal International , 2012 .