S T R U C T U R E OF T H E E A S T E R N COAST BELT, S O U T H W E S T E R N CANADIAN CORDILLERA, F R O M REPROCESSING SEISMIC R E F L E C T I O N LINES 88-17 AND LITHOPROBE 88-14 By Weimin Zhang B . Sc. (Applied Geophysics) Chengdu Institute of Technology A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF T H E REQUIREMENTS FOR T H E D E G R E E OF M A S T E R OF S C I E N C E in T H E F A C U L T Y OF G R A D U A T E STUDIES E A R T H AND O C E A N SCIENCES (GEOPHYSICS) We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA December, 1997 Â© Weimin Zhang, 1997 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Earth and Ocean Sciences (Geophysics) The University of British Columbia 129-2219 Main Mall Vancouver, Canada V 6 T 1Z4 Date: Dec. 'Jtt>, /?<3^ Abstract As part of the Lithoprobe Southern Cordillera research program, seismic reflection data were acquired in 1988 along a 100 km profile across the Eastern and Central Coast Belts (lines 88-14 and 88-17) in a region of complex surface geology. A very good interpretation of the crustal section based on previous processing by industry has been made. However, the data were not processed specifically to enhance the features in the upper 5 s and provide a clear tie to surface geology, some of which is newly studied since the earlier interpretation. The present study involves reprocessing of these data to provide this enhancement. Severe crookedness of the line necessitated the implementation of unconventional processing techniques, including creation of several mini-profiles which cut through the midpoint scatter and application of a first-order correction for effects of reflector crossdip ( i.e., the crossline component of reflector dip). Refraction-based statics were computed using a two-step procedure consisting of initial identification and correction of systematic errors in the picked first break traveltimes, and subsequent application of a 2-D statics algorithm. Dip moveout (DMO) was applied to constant-offset gathers to effect prestack partial migration and reduce the dependence of normal moveout velocity on reflector dip. To reduce possible smearing in the upper few seconds, stacks were repeated with truncated spread lengths. The Karhunen-Loeve transform was applied efficiently and economically on the stacked seismic data to enhance coherent signals in all selected dips. The reprocessing has better delineated the subsurface geometry of the terranes and faults and expressed some of the structures which are not visible on previous processing u sections, even though the 2-D seismic images cannot completely express the 3-D complexity inferred from the surface geology. The reprocessed sections are considered in light of the new geological mapping to provide a refined interpretation of subsurface structure and its implication for local tectonics. A new geological model is developed based on this refined interpretation and is being further tested by seismic forward modeling. in f Table of Contents Abstract " List of Tables vii List of Figures viii Acknowledgement 1 i INTRODUCTION 1 1.1 Research motivation 1 1.1.1 Southern Cordillera Transect 1 1.1.2 Southern Coast Belt in southwestern British Columbia 1 1.1.3 Purpose of this study 5 1.2 1.3 2 x Geological and geophysical background 9 1.2.1 Terrane accretion 9 1.2.2 Coast Belt Thrust System 12 Geophysical framework 15 D A T A ACQUISITION A N D INITIAL PROCESSING 18 2.1 Acquisition and data quality 18 2.2 Hardware and software 19 2.3 Initial processing 22 2.3.1 Filtering 22 2.3.2 Predictive deconvolution 26 iv 2.3.3 3 4 6 31 STATICS C O R R E C T I O N S 37 3.1 Datum correction . . . 37 3.2 Refraction statics 38 3.2.1 Refraction statics algorithm 40 3.2.2 First-break picking 42 3.2.3 Surface stacking charts and first break analysis 43 DIP M O V E O U T A N D V E L O C I T Y A N A L Y S I S 49 4.1 Background and introduction 49 4.2 Implementation of DMO 51 4.2.1 Partial stack 51 4.2.2 DMO algorithm 52 4.3 5 CMP binning Velocity analysis after DMO 57 CROSSDIP C O R R E C T I O N 62 5.1 Introduction 62 5.2 Theory 63 5.2.1 3-D CMP bin traveltime equation 63 5.2.2 Simplified 2-D case 65 5.3 Crossdip correction processing stream 67 5.4 Example of application of crossdip correction 69 5.5 Spread length truncation 69 P O S T - S T A C K P R O C E S S I N G A N D F I N A L SEISMIC S E C T I O N S 6.1 73 Application of the partial Karhunen-Loeve transform to suppress random noise 73 v 7 6.1.1 Introduction 73 6.1.2 Basic theory 74 6.1.3 Computing the filtered image by partial K - L transform 77 6.1.4 Summary and results 79 6.2 Comparison of reprocessed data with contract processed data 81 6.3 Coherency filtering and final stack presentation 82 6.4 Migration 100 I N T E R P R E T A T I O N A N D DISCUSSION 104 7.1 Introduction 104 7.2 Interpretation of reprocessed Line 88-14 109 7.3 Interpretation of reprocessed Line 88-17 112 7.4 Geological model and seismic modelling of Line 88-17 118 7.5 Summary comments 119 References 126 Appendices 135 A Dip-moveout by Fourier transform 135 vi List of Tables 2.1 Acquisition parameters for Lines 88-17 and 88-14 21 4.1 Comparison of velocities from N M O and D M O processing 59 vii List of Figures 1.1 Lithoprobe Southern Cordillera Transect 2 1.2 Lithoprobe reflection profiles in southwestern British Columbia 3 1.3 Location of seismic reflection lines 88-13, 88-14 and 88-17 on a geological map 6 1.4 Terrane map of the Canadian Cordillera 11 1.5 Possible Coast Belt paleogeographies 13 1.6 Crustal cross section of the Southern Coast Belt showing the geometry and kinematics of the C B T S 14 2.1 Acquisition geometry for Lines 88-14 and 88-17 20 2.2 Processing flow for Southern Cordillera Transect Lines 88-17 and 88-14. . 23 2.3 Amplitude spectra of 2 - 6 s for trace 1 in shot gather 389 24 2.4 Comparison of shot gather with and without f-k dip filtering 25 2.5 A flowchart for predictive deconvolution using prediction error filter. . . . 30 2.6 Autocorrelogram for Shot Gather 1139 32 2.7 The Klauder wavelet for both Lines 88-14 and 88-17 32 2.8 The C M P binning for Line 88-17 33 2.9 The C M P binning for Line 88-14 35 3.1 Datum plane and near-surface structure model 38 3.2 Refraction model 39 3.3 Computation of delay times 40 3.4 Example of one shot gather with "pitch out" 44 viii 3.5 Residual traveltime plot for VPs 570 tO 800 47 3.6 Conventional INSIGHT first break traveltimes plot 48 4.1 Travel time t from source s to receiver r 49 4.2 A common midpoint plane 53 4.3 Zero-offset section for the experiment of Figure 4.1 54 4.4 The implementation of DMO algorithm by Fourier transform 55 4.5 Processing stream of velocity analysis with DMO 58 4.6 Example of DMO processing 61 5.1 Plan view showing midpoint scattering along a portion of a crooked survey line 62 5.2 Plan view of the slalom line 64 5.3 Cross-section showing a crossdipping reflector 66 5.4 Processing stream of crossdip correction 68 5.5 Figure shows the same section after application of crossdip correction. . . 70 5.6 Stack with truncated spread length 72 6.1 Flowchart of application of the partial K-L transform 80 6.2 Part of the stacked seismic section of Line 88-14 showing a dipping reflector and background random noise 82 6.3 The reconstruction of the image in Fig. 6.2 following partial K-L filtering 83 6.4 Singular value spectrum of the image in Fig. 6.2 84 6.5 The image reconstruction of the image in Fig.6.2 using 15 images. It attenuates less noise than Fig. 6.3 85 6.6 Comparison of unmigrated seismic section along eastern part of Line 88-17. 86 6.7 Comparison of unmigrated seismic section along western part of Line 88-17. 88 ix 6.8 Upper 6 s of zero-offset stacking section of Line 88-17 92 6.9 Upper 7 s of zero-offset stacking section of Line 88-14 94 6.10 Upper 6 s of final coherency filtered zero-offset stacking section of 88-17. 96 6.11 Upper 7 s of final coherency filtered zero-offset stacking section of 88-14. 98 6.12 Migrated coherency filtered seismic section of Line 88-17 102 6.13 Migrated coherency filtered seismic section of Line 88-14 103 7.1 Lines 88-14 and 88-17 initial interpretation superimposed on migrated and coherency-filtered section produced by the contractor processing 7.2 105 Alternate interpretation of Lines 88-14 and 88-17 made by Monger and Journeay (1994) 106 7.3 Geological cross-section provided by Monger and Journeay (1994) 107 7.4 Recent geologic map of the eastern Coast Belt 108 7.5 Coherencyfiltered,structure stack for the upper 7 s of Line 88-14 110 7.6 Coherencyfiltered,migrated structural section for the upper 7 s of Line 88-14 Ill 7.7 Geologic map with location of Line 88-17 113 7.8 Coherencyfiltered,migrated, structural section for the upper 6 s of Line 7.9 88-17 116 Geological model of Line 88-17 120 7.10 Seismic traveltime profile from normal incidence modelling on the geological model in Figure 7.9 7.11 Interpretation of Line 88-18 by Monger and Journeay 121 123 7.12 Interpretation of Line 88-13 by Monger and Journeay (1994). See Figure 7.3 for geological cross-sections 125 Acknowledgement I would like to thank my supervisor Dr. R . M . Clowes for his guidance and patience over the past several years. He has been teaching me Geophysics and English at the same time. Murray Journeay of the Geological Survey of Canada provided me with new digital geological data and gave me suggestions on geological interpretation. John Amor provided me with lots of assistance in all matters computational. I sincerely thank all professors in the Department for giving me so much knowledge on geophysics and teaching me how to learn. I am grateful to all departmental members for their friendship and interesting discussion. Seismic data 88-14 and 88-17 were acquired primarily with funds from the N S E R C Collabriative Special Projects and Program grant in support of L I T H O P R O B E . Financial support for my student stipend and analysis costs derived from N S E R C Research grants to Dr. R . M . Clowes. xi Chapter 1 INTRODUCTION 1.1 1.1.1 Research motivation Southern Cordillera Transect The Canadian Cordillera is part of one of the great mountain systems on Earth. The scientific program of Lithoprobe's Southern Cordillera Transect is investigating the nature, structure, timing and dynamic evolution of the Cordillera from Vancouver Island to the Rocky Mountains (Figure 1.1). The main emphasis of the program has been to integrate seismic reflection, seismic refraction, and electromagnetic images of deep crustal and lithospheric structure with an expanded understanding of the geology and geochemistry of near-surface rocks, with the ultimate goal of understanding the processes and events that have led to the westward growth of the continent. A summary of multidisciplinary results for the transect has been published recently (Cook 1995). 1.1.2 Southern Coast Belt in southwestern British Columbia The Coast Belt is one of five morphogeological belts of the Canadian Cordillera. It runs for 1700 km between latitudes 49Â° N and 62Â° N and is 100-200 km wide. Some mountain peaks reach elevations above 4 km, although most peaks are in the range of 2-3 km (Monger and Journeay 1992). The belt consists mostly of Late Jurassic to Early Tertiary granites (80-85 % by area according to Roddick 1983); dominant structural features are Mid-Cretaceous to Early Tertiary in age. The Southern Coast Belt (SCB; Figure 1.2), 1 i 8 â€¢1 CD to O <=: o WASHINGTON' 120 r NORTH AMERICAN R O C K S OF E A S T E R N CORDILLERA CENTRAL GNEISS OF. THE C O A S T PLUTONIC C O M P L E X INSULAR S U P E R T E R R A N E LITHOPROBE SEISMIC LINE INTERMONTANE S U P E R T E R R A N E C O C O R P SEISMIC LINE 100 I 200 km Figure 1.1: Lithoprobe Southern Cordillera Transect. It shows the Cordillera in southwestern Canada and the northwestern United States with the location of the major crustal reflection and refraction profiles indicated (from Clowes 1993). Chapter 1. INTRODUCTION 3 Figure 1.2: Lithoprobe reflection profiles in southwestern British Columbia (following page). Figure is taken from Varsek et al. (1993). Reflection profiles 14 and 17 are in the Western and Eastern Coast Belt, respectively. Major tectonic units are also displayed in the figure. The five morphogelogical belts of the Cordillera are displayed in the insert. Abbreviations: Cascade metamorphic core ( C M C ) , Central Nicola horst (CNH), Mount Lytton complex ( M L C ) , Northwest Cascade system (NWCS), Okanagan plutonic complex (OPC), Ashlu Creek fault ( A C F ) , Ashcroft fault system (AF), Bralorne fault system (BF), Beaufort Range fault ( B R F ) , Central Coast Belt detachment ( C C B F ) , Cowichan Lake fault (CF), Castle Pass fault ( C P F ) , Chuwanten fault (Chwf), Crescent fault (CrF), Coldwater fault ( C W F ) , Duffy Lake fault (DLF), Fraser fault ( F F ) , Hozameen fault (HF), Harrison Creek fault (HCF), Jack Mountain fault ( J M F ) , Kwoiek Creek fault ( K F ) , Marshall Creek fault ( M C F ) , Miller Creek fault ( M C F ) , Mocassin Lake fault ( M L F ) , Mission Ridge fault ( M R F ) , Methow River fault (MtF), Okanagan Valley fault (OVF), Owl Lake fault (OLF), Pasayten fault (PF), Phair Creek fault ( P C F ) , Quilchena fault (QF), Ross Lake fault ( R L F ) , Red Shirt fault (RSF), Shuskan fault (SF), Thomas Lake fault ( T L F ) , Tyaughton Creek fault (TyF), West Coast fault (WCF), Yalakom fault ( Y F ) . 8 Plutons (Meso*., Cenoz.) ^ ' â„¢ , n $ U , a r C o m x , s ! , e T9 9 North American Metasediments FJS] CtifJ Cascadia Bart Composite T.rrana '.'.â€¢.'J Intermontane â€¢ â€¢ Composite Terrane 1 pa ^ lii&j a n A â€ž 2 omblage d Successor Basin O to o o Chapter 1. INTRODUCTION 5 the focus of this study, is defined as the region between 49Â° N and 51Â° N, bounded on the west by the Strait of Georgia and on the east by the Pasayten fault (Wheeler and McFeely 1991). The SCB is divided by Journeay and Friedman (1993) into three tectonic elements, namely, the Eastern, Central, and Western Coast belts (ECB, CCB, and WCB, respectively). As part of the Lithoprobe Southern Cordillera Transect research program, nearly 350 km of deep seismic reflection data (Figure 1.2) were acquired across the southern Coast Belt, these linking previous studies in southeastern British Columbia (Cook et al. 1988) and central British Columbia (Cook et al. 1992) to Vancouver Island (Clowes et al. 1987). Varsek et al.(1993) provide an interpretation of the reflection dataset as processed by an industry contractor (lines 12-18, Figure 1.2). 1.1.3 Purpose of this study Since the interpretation by Varsek et al. (1993), new geological mapping in the eastern Coast Belt has been combined with geochronology for that region and other parts of the Coast Belt (Journeay and Friedman 1993). While a very good interpretation of existing data was made, we believe that reprocessing of parts of the dataset will enhance image quality to the extent that, when combined with the new geological information, an improved and more detailed interpretation will be possible. Lines 88-14 and 88-17 cross the central and eastern Coast Belt in a region of complex surface geology, including multiple faults and terrane fragments (Figure 1.3). Data quality is good. Originally the data were processed robustly to provide good crustal images, but not necessarily good images from 0 to 5 s to enable a clear tie to surface geology. Thus, our principal objective in reprocessing lines 88-14 and 88-17 is to provide better definition in the upper part of the seismic section and to tie results closely with surface geology, some of which is newly studied since the earlier interpretation. Chapter 1. INTRODUCTION 6 Figure 1.3: Location of seismic reflection lines 88-13, 88-14 and 88-17 on a geological map. (A) Location of seismic reflection lines 88-13, 88-14 and 88-17 on a geological map of the Coast Belt (from Wheeler and McFeely 1991). Shot point numbers are identified for lines 88-17 and 88-14. (B) Legend for geological map in (A). Chapter 1. INTRODUCTION 8 LEGEND PLUTONIC ROCKS STRATIFIED ROCKS Pleistocene and Recent Barry Tertiary (40 - 64 Ma) IIP!!!! Thick Drift Early and Late Cretaceous Cretaceous Pasayttn Group KP X P V :5 Mid-Cretaceous granodiorifc (gd) in Coast Belt Virginian Ridge Formation Early Jurassic - Late Jurassic (145 - 187 Ma) Gambler Group KGB KGP KT Granodiorifc (gd) and quartz dioritc (qd) Winthrop Formation ffljj^f KG Granodiorifc (gd) in Western Coast Belt Brokenback Hill Formation ! Peninsula Formation Taylor Creek Group OPHIOLITE COMPLEXES Permian Bralome Complex Jackass Mountain Group Lower Jurassic to Lower Cretaceous BBt.Hal|HE Graoodiorite and monzograoite Paleozoic and Mesozoic Shulaps Complex Cayoosh Assemblage Early md Middle Jurassic Harrison Lake Formation Bmm Last Creek Formation METAMORPMC ASSEMBLAGES Early Late and Late Cretaceous j â€¢â€¢ C a d w a l Â«Jer Gronp A Chism Creek Schist Undifferentiated garoet-biotifc Hurley Formation Pioneer Formation Carboniferous to Middle Jurassic Bridge River Complex Radioutrian chert, siltstonc Light to dark grey pbyllite (B) Chapter 1. INTRODUCTION 9 To achieve these improved images, a number of special procedure were used. Surfaceconsistent refraction statics were applied and several different binning geometries which cut through the midpoint scatter were considered. Severe crookedness of the line necessitated the implementation of a first-order correction for 3-D effects of reflector crossdip. Dip moveout (DMO) was applied to effect prestack partial migration and reduce the dependence of normal moveout velocity on reflector dip. Offset-limited stacks controlled the negative contributions of extensive lateral variations along the travel paths, reduced surface noise levels and enhanced continuity of reflections. The higher quality seismic images could then be related more clearly to the complex near-surface geology. 1.2 1.2.1 Geological and geophysical background Terrane accretion The Canadian Cordillera is part of the orogen in which "terranes" were first recognized as fundamental building blocks of the continent. The concept of a "terrane" is commonly involved in describing the tectonic development of the Canadian Cordillera. A terrane is a crustal block or fragment that preserves a distinctive geologic history that is different from the surrounding areas and that is usually bounded by faults. Accreted terranes are those that become attached to a continent (in this case North America) as a result of tectonic processes. It was recognized later that many terranes could be regarded as enormous thrust sheets which, stacked one upon the other, make up much of the Cordilleran crust (Monger 1985). The various terranes of the Canadian Cordillera are shown in Figure 1.4. This present structure is considered to be the result of the accretion of two allochthonous superterranes - the Intermontane superterrane and Insular superterrane - to the leading edge of the westward-travelling North American plate prior to mid-Cretaceous time (e.g., Monger and Price 1979, Coney et al. 1980). The Intermontane superterrane, which was Chapter 1. INTRODUCTION 10 positioned adjacent to ancestral North America, consists primarily of the Stikine, Bridge River, Cache Creek, and Quesnel terranes, while the Insular superterrane comprises the Alexander and Wrangellia terranes. From a tectonic perspective, the Coast Belt is a critical component of the Canadian Cordillera because it masks relationships between the Insular and Intermontane superterranes. Numerous detailed explanations of the tectonic processes involved in creating the Coast Belt through the accretion of geologic terranes have been provided (e.g., Davis et al. 1978, Monger et al. 1982, Gehrels and Saleeby 1985, and van der Heyden 1992). The generally accepted view is that of Monger et al. (1982), who described the Coast Belt as containing the suture that resulted from the mid-Cretaceous collision between an exotic Insular superterrane and the Intermontane superterrane, which was by then attached to North America. The collision followed the subduction-related closure of an intervening ancient ocean basin, vestiges of which are presumably contained in the rocks of the oceanic Bridge River and metasedimentary Methow terranes. The suture has since been overprinted by a magmatic arc related to a post-collisional east-dipping subduction zone. However, van der Heyden (1992) suggested that the Coast Belt magmatic intrusions are due to a prolonged period of east-dipping subduction ongoing since the mid-Jurassic, at which time a single "megaterrane" consisting of Wrangellia, Alexander and Stikinia was accreted to ancestral North America along a suture which is today represented by rocks of the Cache Creek and Bridge River terranes. Most recently, Monger et al. (1994) tried to reconcile the differing view points. They discussed the dominantly plutonic Canadian Coast Belt as flanked on the northwest by a Late Jurassic-Early Cretaceous intra- or back-arc basin. In the south it contains strata of the same age deposited in an accretionary fore-arc complex, east of a coeval arc. This distribution is ascribed to pre-mid-Cretaceous sinistral movement of part of the Jurassic-Cretaceous North American plate margin arc to a position west of the accretionary fore-arc complex. Figure 1.5 Chapter 1. INTRODUCTION TERRANES IN CANADIAN CORDILLERA accreted terranes shown by (A); pericratonic terranes by (P); terranes with oceanic basement but cratonic detritus by (A and P); oceanic, marginal basin, and continental margin accretionary complexes in boldface ALEXANDER (A) BRIDGE RIVER (A) CASSIAR(P) CACHE CREEK (A) CADWALLANDER (A) CHUGACH (A and P) CHILLIWACK and HARRISON (A) CRESCENT (A and P) HO (A and P) JUAN DE FUCA OCEANIC PLATE KOOTENAY (P) MONASHEE (P) METHOW (A) NISLING (P) OZETTE (A and P) PACIFIC OCEANIC PLATE PACIFIC RIM (A and P) QUESNEL (A) SHUKSAN (A and P) SLIDE MOUNTAIN (A) STIKINE (A) TAKU (A) WRANGELLIA (A) YAKUTAT (A and P) YUKON-TANANA (A and P) 500 km Figure 1.4: Terrane map of the Canadian Cordillera (adapted from Clowes, 1990) Chapter 1. INTRODUCTION 12 shows their suggested Coast Belt paleogeographies. Their hypothesis showed how both intra-arc and inter-arc basins may have been accommodated in a single orogenic system prior to Late Cretaceous and early Tertiary tectonism. It explains the termination of Jurassic-Cretaceous arc rocks and their Wrangellian basement at the south end of the Coast Belt near 49Â° N , and apparent width of the Jurassic-Cretaceous arc at this latitude. And finally they suggested that the Insular superterrane was emplaced by sinistral movements along the margin of the Pacific Ocean basin, rather than by migration across it. 1.2.2 Coast Belt Thrust System The Coast Belt Thrust System (CBTS) is a recently defined component of the Coast Belt (Journeay and Friedman 1993). It represents the leading edge of a west-vergent contractional belt that formed along the inboard margin of the Insular superterrane in late Cretaceous time. The Western Coast Belt ( W C B ) , which consists of voluminous midJurassic to mid-Cretaceous plutons that have intruded mid-Triassic to early Cretaceous arc sequences, is of low metamorphic grade. Except in the vicinity of its eastern flank, the W C B is presumed to have acted as a rigid crustal block during two stages of vigorous late Cretaceous shortening that took place along a west-vergent contractional belt, the C B T S (Figure 1.6; Journeay and Friedman 1993). In the study area, the C B T S extends across strike for approximately 35 km, encompassing both the eastern flank of the W C B as well as the western portion of the Central Coast Belt ( C C B ) , where it is manifest in a stack of east-dipping thrust sheets which carry highly metamorphosed deep levels of the Eastern Coast Belt ( E C B ) in their hanging walls. The C C B and E C B feature metamorphosed oceanic rocks of the Bridge River and Methow terranes, as well as island arc sequences of Cadwallader terrane affinity. The C B T S is defined on its eastern limit by the Central Coast Belt detachment ( C C B D ) . Chapter 1. INTRODUCTION 13 Figure 1.5: Possible Coast Belt paleogeographies. A: Continuity of Jurassic and Cretaceous plate margin arcs and accretionary complexes; restored offset on possible sinistral fault about 800 km. B: Prior to Late Cretaceous contraction; evidence for sinistral movement suggested distribution of shear across entire accretionary complex, not along single fault shown. C: Late Cretaceous and younger; displacement restored: 450 km contraction across North Cascades (McGroeder 1991), 165 km across Bowser basin (B), and 130 km on Fraser-Straight Creek fault (Monger and Journeay 1993). Not restored are > 450 km (based on geological evidence) of orogen-parallel dextral displacement in north-central British Columbia-Yukon, or 1500-3000 km of dextral displacement across entire Canadian Cordillera (based on paleomagnetic evidence), and 100-300 km of shortening (based on geologic evidence) in easternmost Cordillera. K is Kahiltna, G is Gravina basin, N is Nutzotin. (figure is adapted from Monger et al. 1994). apterl. INTRODUCTION 14 Figure 1.6: Crustal cross section of the Southern Coast Belt showing the geometry and kinematics of the CBTS. Lines 88-14 and 88-17 were run between TLF and BF. Abbreviations: Ascent Creek fault (ACF), Breakenridge fault (BF), Central Coast Belt Detachment (CCBD), Fire Creek fault (FCF), Fitesimmons Range fault (FRF), Thomas Lake fault (TLF), Terrarosa thrust (TT). (from Journeay and Friedman 1993) Chapter 1. 1.3 INTRODUCTION 15 Geophysical framework Although earlier geophysical studies helped characterize the crustal structure of the Canadian Cordillera, more recent studies carried out as part of the Southern Cordillera Transect scientific program have established the current geophysical framework. Crustal thickness for the southern Cordillera and northwestern Washington, derived from refraction data, is estimated to vary from 30 to 48 km (Clowes et al. 1995). Crustal thickness of the Coast Belt is between 34 and 37 km ( B . C . Zelt et al. 1993). The 3D velocity structure of the lithosphere in the southern Canadian Cordillera has been mapped in considerable detail and demonstrates the strong degree of 3-D heterogeneity within the crust and upper mantle. A first-order characteristic is the continuous increase in crustal velocities westward from the Foreland Belt to the Insular Belt (Clowes et al. 1995). The Southern Cordillera Refraction Experiment (SCoRE) was executed during the summers of 1989 and 1990. SCoRE'89 focussed on the Intermontane and Coast Belts, with an incursion into the eastern Insular Belt (B.C. Zelt et al. 1992); SCoRE'90 completed studies in the Coast Belt and extended coverage eastward into the Omineca and Foreland Belts (Kanasewich et al. 1994; Zelt and White 1996). Line 10 was recorded along the Eastern Coast Belt. O'Leary et al. (1993) discuss the crustal and upper mantle P-wave velocity structure determined by 2-D traveltime inversion and amplitude forward modelling. A 2-D velocity structure comprising three intracrustal reflectors was developed for the central Coast Belt. Their discussion suggested that three types of reflecting horizons may exist: (1) structural features (e.g., fault zones), (2) transition zones at the top and bottom of a region of layered porosity within the crust, and (3) lithological contrasts (e.g., the boundary of a batholith and the surrounding terrane material). Chapter 1. INTRODUCTION 16 Moho reflections (PmP) from refraction line 10 approximately coincide with the nearvertical incidence reflection Moho on both lines 88-17 and 88-13. The reflection Moho observed on lines 88-13 and 88-17 seismic sections appears as a band of enhanced reflectivity with two-way traveltime varying from 0.2 to 0.8 s in depth extent, with its top corresponding to the top of the wide-angle Moho from the refraction survey. Line 88-14 did not delineate a Moho reflector. Cook et al. (1995) present analyses of gravity and magnetic data along the corridor of the Lithoprobe Southern Cordillera Transect, which allowed characterization of anomaly patterns according to their likely sources. Long-wavelength Bouguer gravity anomalies are attributed to isostatic effects of topography, which in most areas is compensated. A large, linear magnetic high in the western Coast Belt may be caused partly by deep-seated magnetic variations associated with active subduction and depression of the isotherms, but Clowes et al. (1997) show that the high frequency, high amplitude anomalies can be modelled by upper crustal bodies with high susceptibility. The value of these data is best appreciated when they can be combined with both surface variations delineated by mapping and subsurface variations outlined by seismic reflection techniques. In many areas, this type of approach can significantly reduce the inherent ambiguities that usually plague analyses of potential field data, and provide important limitations on hypotheses of lithospheric structure and evolution. Various other geophysical studies have been carried in the Southern Cordillera. Cassidy (1995) showed receiver function analysis from teleseismic study as a powerful tool to constrain the S-velocity structure of the crust and upper mantle across the southern Canadian Cordillera. The site-specific nature of this method makes it ideal for tracking structure. Studies to date have utilized three-component broadband data from a portable array deployed across southwestern British Columbia and broadband data from the new Chapter 1. INTRODUCTION 17 Canadian National Seismograph Network. Analysis of these data demonstrates the variations in crustal thickness and complexity across the southern Canadian Cordillera, with the Moho depth varying from 35 km in the Coast Mountains, to 33 km near Penticton, to 50 km near the Rocky Mountain deformation front. Large programs of electromagnetic studies, using principally the natural-source magnetotelluric (MT) technique, were initiated and carried out within southern and central British Columbia during the last decade. Jones and Gough (1995) compiled MT data from the largest number of sites ever recorded for the study of accretionary tectonics, and showed the lateral and depth variations of MT parameters and electrical conductivity within the southern and central Canadian Cordillera. Regionally, much of the crust of the Canadian Cordillera exhibits electrical conductivities two to four orders of magnitude greater than what one would expect from dry, competent rocks which make up the crust. They interpret the Canadian Cordillera Regional conductor as saline fluids trapped at depth below the brittle - ductile transition, and that internal conductivity variations are a consequence of fracture density variation. Hyndman and Lewis (1995) summarized the surface heatflowfor the Southern Cordillera Transect and explain that the crustal thermal regime has important implications for the interpretation of the deep seismic structure. They show high heat flow in the Coast Belt; and more specifically very high and irregular heatflowin the Garibaldi Volcanic belt. Based on their study, the very high crustal temperatures in the region of the Garibaldi Volcanic belt indicate that the mid-crustal seismic reflector there is located near the regional brittle-ductile transion. In general, the WCB is characterized by several large geophysical contrasts typical of convergent margins, whereas the ECB has geophysical anomalies that are more typical of extended continental crust. Chapter 2 DATA ACQUISITION A N D INITIAL P R O C E S S I N G 2.1 Acquisition and data quality L I T H O P R O B E reflection data are acquired by an industry seismic contractor and initial processing also is carried out under contract to a processing company. Initiation, monitoring and quality control of the contract are undertaken by L I T H O P R O B E scientific personnel. Deliverable products from the contractor include formal reports and tapes of shot gathers, common midpoint ( C M P ) gathers, structural stacks and migrated stacks. A l l of the information and data are available for reprocessing of any part of the data set. As noted previously, this thesis is concerned with two lines, 88-14 and 88-17, recorded in the Coast Belt of southern British Columbia. Continuous across-strike coverage was restricted by poor road access in the western part of the rugged Coast Mountains; on several profiles elevation changes approached or exceeded 1 km. The tendency of major valleys to be oriented along structural trends resulted in some coverage with strike orientation (Figure 1.2). The structural and metamorphic culmination of the eastern Coast Mountains is traversed in part by three profiles (14/17, 13 and 18) that are each spaced about 60 km apart; this geometry allows investigation of along-strike variations in this complex region. However, this thesis only deals with Lines 88-14 and 88-17, which provide the most northerly crossing of the axial structural belt of the Coast Mountains (Figure 1.2). Sonix Exploration Ltd. performed the data acquisition. The 240-channel data set 18 Chapter 2. DATA ACQUISITION AND INITIAL PROCESSING 19 was acquired using a 'master/slave' configuration employing two 120 channel DFS-V recording units, the 'master' recording traces 1-120, and the 'slave' recording traces 121240. The asymmetric split spread recording geometry is illustrated in Figure 2.1. Acquisition parameters are listed in Table 2.1. The large number of channels allowed recording and subsequent processing of (nominally) 60-fold CMP (common midpoint) data, with a CMP spacing of 25 m. Crooked line orientations caused variable fold, sometimes as high as 100-fold, or as low as 30-fold. As evidenced by the processed structure stacks and migrated stacks produced by Western Geophysical Ltd., the data quality was generally good but varied over the lines. The contract processing did not apply some of the more expensive and time consuming processes such as refraction statics and dip moveout procedures. To enhance the seismic images along Lines 88-14 and -17, reprocessing of the data from the basic shot gathers was necessary. Performing this task was the principal effort of the thesis study. 2.2 Hardware and software A commercial seismic processing software package, INSIGHT from Inverse Theory and Application (ITA) Ltd., installed on a Sun microsystems SPARCstation 10 provided the basis for most of the application processing in this study. Hardcopy output could be directed to a 36" CalComp plotter, 11" Tektronix color plotter or a laser printer, as best suited for the specific need. A second commercial package, OUTRIDER from MicroSeis Technology Ltd., provided the facility to build geological structure models with assigned velocity-density characteristics and compute synthetic seismogram sections for comparison with observed record sections. Chapter 2. DATA ACQUISITION AND INITIAL PROCESSING Figure 2.1: Acquisition geometry for Lines 88-14 and -17. 20 Chapter 2. DATA ACQUISITION AND INITIAL PROCESSING Source 4 vibrators Energy source Vibrator model L.R.S. 315 Sweep frequency 10-56 Hz linear Sweep length 14 s Taper 0.5 s Number of sweeps 8 Source interval 100 m Drag length 90 m Receivers Mark Products, L10A Geophone type Geophone frequency 10 Hz Geophone array 12 over 50 m Number of groups 240 Group interval 50 m Coverage 60 fold Instruments DFS-V FT-1 Model 48db Pre-amp gain Gain mode IFP Sample interval 4 ms 32 s Record length Correlated record length 18 s SEG-B Recording format Tape density 1600 B.P.I. Table 2.1: Acquisition parameters for Lines 88-17 and 88-14. 21 Chapter 2. DATA ACQUISITION 2.3 AND INITIAL PROCESSING 22 Initial processing The complete processing stream applied to the two data sets is shown in Figure 2.2. The data obtained from the contract deliverables were demultiplexed common shot gathers in standard S E G - Y format. Our initial processing of these data included: (1) set-up geometry, (2) spherical spreading correction, (3) first break picking for refraction statics analysis, (4) filtering, (5) predictive deconvolution and (6) C M P binning. For (1), the geometry was set up based on survey file and observers' notes. For (2), the correction was applied based on a single regional velocity-time function that we derived from the velocity analysis information from Western Geophysical Ltd.'s initial processing. Step (3) is discussed separately in Chapter 3 because of its importance for improving phase coherency of reflectors. In the following sections, we will describe steps (4), (5) and (6). In subsequent chapters, details of the most important processing steps outlined in Figure 2.2 (highlighted by the boxes with shadows) are discussed. 2.3.1 Filtering A seismic trace or seismogram is composed of desired signal and undesired signal or noise. Our goal is to reduce the undesired signal and noise as much as possible to enhance the reflected events so they can be interpreted for their geological significance. Since the amplitude spectra of the signal and noise are often different, we can separate them on the basis of frequency. First of all, we applied a zero-phase frequency bandpass filter to the data set. It is performed as three steps: take the Fourier transform of the input data, multiply the quantity by specified bandpass weights, apply the inverse Fourier transform and return the filtered result. In our data set, amplitudes between 10.0 Hz and 55.0 Hz are 90% of the passband value, while amplitudes at frequencies less than 5.0 Hz or larger than 80.0 Chapter 2. DATA ACQUISITION AND INITIAL _c PROCESSING 23 IL Survey file SEG-Y gathers Observers' notes Pre-processing - geometry set-up - trace edit - first break, picking Spherical spreading correction F-K filtering Deconvolution Refraction statics CMP binning Prelimnary velocity analysis DMO correction Final velocity ananlysis NMO correction Crossdip correction Residual and trim statics CMP stack | F-X decon K-L filtering Data merge Coherency filtering I Migration Figure 2.2: Processing flow for Southern Cordillera Transect Lines 88-17 and 88-14. Boxes with shadows highlight processes for which details are discussed in the text. Chapter 2. DATA ACQUISITION 00 ijo ioTo 3~o!o AND INITIAL PROCESSING 4tIo SO~XI 600 7^0 800 9<[o 24 100.0 1l6.0 120.0 Frequency (Hz) Figure 2.3: Amplitude spectra of 2 - 6 s for trace 1 in shot gather 389. The strong amplitude around 48 Hz observed in many traces is a periodic noise associated with the proximity of receivers to the recording trucks. Hz are 10% of the passband value. We use this broad filter because we want to maintain as much bandwidth as possible while still rejecting high and low frequency noise. A periodic noise pattern was observed on many traces. Figure 2.3 shows an example. It was associated with the proximity of receivers to the recording trucks. The observation suggested a selective application of the filter based on individual trace characteristics. A simple zero-phase notch filter succeeded in removing a significant portion of the noise. A significant number of traces were contaminated by direct shear wave arrivals and surface waves as a linear pattern of coherent noise which cannot be removed by bandpass filtering. These types of noise usually are isolated from reflection energy in frequencywavenumber (f-k) space. Given the velocity and frequency cutoff, the noise could be attenuated by INSIGHT module L F A F . However, the typical time dips associated with these events were sufficiently steep to cause spatial aliasing. Spatial aliasing often causes Chapter 2. DATA ACQUISITION AND INITIAL PROCESSING 25 Figure 2.4: Comparison of shot gather with and without f-k dip filtering. After filtering, coherent noise as indicated by the arrows in (A) is attenuated such that a shallow reflector (indicated by arrows) can be observed in (B). poor f-k dip filtering . Yilmaz (1987) suggests that a practical approach to this problem is to apply time shifts to the data before f-k filtering so that the unwanted noise appears at lower dips, which are extended to the spatially aliased frequency components, thus eliminating the spatial aliasing effects. The time shifts then are removed after f-k filtering. Following this procedure, we did our f-k filtering. Figure 2.3 shows one shot gather before and after f-k filtering. Chapter 2. DATA ACQUISITION AND INITIAL PROCESSING 26 This procedure may not work in all data sets, since events not spatially aliased before may be spatially aliased after linear moveout. From the experience of my processing, I noticed that guided waves, which manifest themselves as linear noise on both commonshot and CMP gathers, are largely attenuated by stacking. 2.3.2 Predictive deconvolution ^ Geologic boundaries in the subsurface generally do not give rise to sharp reflections on a seismogram. In general reflections are spread out in time and are followed by reverberations which interfere with reflections from deeper layers and may even obscure them. Deconvolution enjoys widespread use by seismic data processors. It is a process that improves the temporal resolution of seismic data by compressing the basic wavelet. It is an operation which takes the seismic trace x as input and gives an estimate of the reflectivity r as output. Many ways have been proposed to do this and most have in common the fact that they form the convolution of the data trace x with an "inverse filter" or "deconvolutionfilter"and the output of that process is the deconvolved trace or reflectivity estimate. All approaches to deconvolution assume the convolutional model of the seismic trace and the methods can be divided into 2 basic types: those which assume the wavelet is known, and those which assume the wavelet is unknown (commonly done using a process called predictive deconvolution). In this study, I performed a time-variant predictive deconvolution on common shot gathers. Actually, predictive deconvolution is one type of Wienerfiltering(Robinson and Treitel 1980). It can be explained in the following way: Given the input x(t), we want to predict its value at some future time (t+a), where a is the prediction lag. The assumption is that a seismic trace often has a predictable component (multiples) with a periodic rate of occurrence. Other features, such as genuine reflections, are unpredictable. So the error series, the difference between actual output and predictive output, is essentially the Chapter 2. DATA ACQUISITION AND INITIAL PROCESSING 27 reflection series. Why did we choose predictive deconvolution? The convolutional model for the noisefree seismogram is: x(t) = w(t) * e(t) (2.1) where x(t) is the seismogram, w(t) is the seismic wavelet, and e(t) is the impulse response. Basically, we want to recover the reflection spike. Our desired output, the impulse response, is the convolution of the inverse filter with the known seismogram, referred to as spike deconvolution: e(t) = w'(t) * x(t). (2.2) Unfortunately, such a procedure is based on the assumption of a priori knowledge of the wavelet form. In practice, it is impossible to acquire knowledge of the reverberation and attenuation components of the wavelet. Converting the seismic wavelet into a spike is like asking for a perfect solution. However, because of noise in the seismogram and assumptions made about the seismic wavelet and the recorded seismogram, spiking deconvolution is not always desirable. Mathematically, spiking deconvolution is a special case of predictive deconvolution, when the predictive lag a is set to one time unit. The predictive deconvolution is performed by calculating the autocorrelation of the input seismogram and using this information to design the prediction filter or prediction error filter. The algorithm, which is based on Wiener filter theory, is outlined here. The general form of the matrix equation for an optimum Wiener filter of length n is given by (for details, see Robinson and Treitel 1980; Yilmaz 1987): Chapter 2. DATA ACQUISITION AND INITIAL PROCESSING r 0 r r\ r r t A 1 a 28 \ \ 0 ' nâ€”1 2 (2.3) r \ ^ 5n-l y where r,- are the autocorrelation coefficients of the input series, aj are the desired filter j 0 / coefficients, and gi are the crosscorrelation coefficients of the desired output with the input series. For the special case of spiking deconvolution (the idealized solution), the desired output is (1, 0, 0, 0), so # = (1, 0 , 0 , 0 ) . To get thefiltercoefficients a* for the estimate of the seismogram x(t + a) at time (t + a), the gi are reformulated in terms of other known parameters. For the prediction lag a, gi is equal to r-j , then, rewrite +a equation (2.3) as: \ / r â€¢ 0 7*1 ao r -i \ / r a \ n r 2 (2.4) r-n-2 r j \ â€¢ 0-n-l 0 / \ ) The predictionfiltercoefficients aj can be computed from equation (2.4). The predictive output yi can be obtained by convolving the predictionfiltera; with the input yi = Xi* (2.5) ^ The prediction error series e;, e,- = x i + a - yi (2.6) The remaining unpredictable part, the difference between the desired output and the predicted output, is essentially the reflection series. A simple example is illustrative. Let i=0, 1, 2, 3 and a=2. Then y 0 = a Xo 0 Chapter 2. DATA ACQUISITION AND INITIAL PROCESSING = yi 29 a\x + a x\ 0 0 + a\X\ + 2/2 = a,2X V3 03X0 -f 02^1 0 = ax 0 + 01^2 2 + 0,QX 3 and e 0 = Â£ e = xi e2 = x â€” ax ez = Â£3 â€” aia;o â€” a o ^ i e4 = 0 â€” a2a;o â€” a,\X\ â€” agx x ee, = From these, 0 2 0 0 2 0 â€” azXo â€” a x\ â€” a\x â€” a^xz 2 2 can be obtained by convolution of a filter with coefficients (1, 0, â€” a , 0 â€” o i , â€” a , â€”03) with Xj. Yilmaz (1987) described this filter as the prediction error filter, 2 which is (1, 0, 0,... â€” a , â€” a i , â€” 0 a _i). n So there are two approaches for predictive de- convolution. The predictive filter can be designed by equation (2.4) and applied through equation (2.5) and (2.6). Alternatively, the prediction deconvolution can be performed when the prediction error filter is designed and convolved with the input series, which reduces to a four-step calculation. The flowchart for predictive deconvolution using prediction error filters is shown on Figure 2.5. As noted before, given an input wavelet of a + n samples, predictive deconvolution using a prediction filter with length n and prediction lag a converts this wavelet into another wavelet that is a samples long. The first a lags of the autocorrelation are preserved, while the next n lags are zeroed out. Now, consider the real situation of an unknown source wavelet. Based on the assumption that reflectivity is a random process, the seismogram has the characteristic of the Chapter 2. DATA ACQUISITION AND INITIAL PROCESSING 30 Input Step 1: Compute (n+a) - lags of autocorrelation, where n= length of the prediction filter and cc= the prediction lag. Step 2: Compute the prediction filter series &2,ani). (ao, a i , Step 3: Design the prediction error filter y delaying the prediction filter: D (1, 0, 0,0, -ao, - a i , - a n i ) ^lf T Step 4: * i Output Figure 2.5: A flowchart for predictive deconvolution using prediction error filter. The symbol "*" refers to convolution. Chapter 2. DATA ACQUISITION AND INITIAL PROCESSING 31 seismic wavelet since their autocorrelation and amplitude spectra are similar. Autocorrelation of input seismograms rather than that of the seismic wavelet is used to design the deconvolution operator. The first isolated burst of energy that is common to all the traces should represent the autocorrelation of the wavelet (Yilmaz 1987). So, the length of the wavelet may be estimated by the length of the first transient zone on the autocorrelation. As exemplified in Figure 2.6, the wavelet length for our data is estimated to be 150 ms, which is the value of a + n. The predictive deconvolution process invokes a minimum-phase assumption, but the recorded wavelet in Vibroseis data is mixed phase. The hope for gapped predictive deconvolution is that if a suitable predictive distance (predictive lag) can be determined, earth-filtering effects can be treated without need of a phase-correction filter (Pollet et al. 1982). The optimum prediction distance in a given case will depend upon many factors, such as the level of noise and details of earth filtering. Following the suggestion of Gibson and Lamer (1984), the prediction lag a is the effective width of the central peak of the Klauder wavelet, which is 35 ms in this study (Figure 2.7). Klauder wavelet is the autocorrelation of vibroseis sweep signal. Sweep frequency, sweep length and taper are shown in Table 2.1. 2.3.3 C M P binning In order to process the dataset using two-dimensional straight line techniques, neighboring midpoints are grouped into bins and each bin is treated as a single C M P in the subsequent processing. We constructed several slalom lines, which ideally cut through the centre-of-mass of the midpoint scatter, for each profile, using a bin width of 25 m (equal to the nominal C M P interval) and a bin length of 1000 m. The C M P bin map for fines 88-17 and 14 are presented in Figure 2.8 and Figure 2.9, respectively. The related problem of crossdip associated with the crooked fines is described in Chapter 5. Figure 2.6: Autocorrelogram for Shot Gather 1139. Inferred wavelet width is indicated by the arrow in the upper left hand corner of the figure. -0.20 -0.16 -0.12 -0.08 C/5 -0.04 U 0.00 s 0.04 0.08 0.12 0.16 0.20 -1.0 Trace Amplitude Figure 2.7: The Klauder wavelet for both Lines 88-14 and 88-17. Chapter 2. DATA ACQUISITION AND INITIAL PROCESSING 33 Figure 2.8: The CMP binning for Line 88-17. The Line 88-17 profile (thick black line) is displayed in the mainfiguretogether with scattered source-receiver midpoints (blue dots). Station number are indicated in black. The slalom line (redfine)is defined by a series of points (red dots). Insert provides an expanded view of midpoint scatter near station 964 together with CMP bins. Because each bin acts as a single CMP in the subsequent processing, the algorithm attempts to tie bin numbers (red) to station numbers (black) to within a scale factor of ten. Chapter 2. DATA ACQUISITION AND INITIAL PROCESSING 35 Figure 2.9: The CMP binning for Line 88-14. The Line 88-14 profile (thick black line) is displayed in the main figure together with scattered source-receiver midpoints (blue dots). Station number are indicated in black. The slalom line (red line) is defined by a series of points (red dots). Insert provides an expanded view of midpoint scatter near station 320 together with CMP bins. I T I A L PROCESSING Chapter 3 STATICS CORRECTIONS The presence of shallow, laterally varying low-velocity weathering anomalies can cause serious problems in the processing of seismic data. Statics correction are applied to compensate for traveltime anomalies introduced by variations in near-surface structure. They comprise datum, refraction statics, residual and trim statics corrections. 3.1 Datum correction Datum correction refers the corrections to a specific datum plane (Figure.3.1). The elevation static at the kth station, e*., is given by ek = - ^ ( 3 - 1 } where z(xk) is defined as the distance between the surface and the datum plane, and V i is the replacement velocity. At station k rep z{x ) - h(x ) - d. k k (3.2) where h(x) is the surface elevation and d is the datum plane elevation. In this study, d was assigned the average of all stations' elevations. We seek to replace the near-surface structure by a medium of uniform velocity, V i, and constant surface topography, d. rep The correction was performed by means of the INSIGHT module "DATM". 37 Chapter 3. STATICS CORRECTIONS 38 Low Velocity layer Vw Bedrock Velocity V ](x) rep Figure 3.1: Datum plane and near-surface structure model. 3.2 Refraction statics Over the last few years, refraction statics have been used with great success to correct for the effects of near-surface weathering anomalies. As indicated in Figure 3.2, these anomalies are largely due to the presence of shallow, low velocity layers whose depths and/or thicknesses vary rapidly enough to cause significant distortion of seismic energy passing through them. This distortion manifests itself as both an error in the structural times of the deeper layers and as a blurring of the stacked image. The refraction statics procedure used here is based on a two-dimensional model consisting of one or more weathering layers, with variable thickness and variable surface topography overlying the bedrock refractor. The refraction analysis models the first breaks as head waves which propagate along the interface between the weathering layer and bedrock at the bedrock velocity (Figure 3.2). A n important question in estimating shot and receiver statics is accuracy of the results Chapter 3. STATICS CORRECTIONS 39 V,(X) Figure 3.2: Refraction model. as a function of the wavelengths of statics anomalies. Following the datum corrections, an implementation of statics corrections is applied in two stages: (1) refraction-based statics corrections to remove long wavelength anomalies and (2) reflection-based residual statics corrections to remove any remaining short wavelength static shifts. Residual statics corrections are needed because field statics and datum corrections almost never totally compensate for the effects of near-surface velocity variations. This is because the near-surface velocity variations are not known and therefore exact corrections cannot be made. The method of surface-consistent statics estimation based on reflections works well in accommodating short-wavelength variations, but performs poorly in handling the long-wavelength variations (Wiggins et al. 1976). The main reason for this is that the input to reflection-based statics algorithms is arrival time difference between traces, not absolute times. Such residual statics programs may improve the stacking quality while leaving errors in structural times (Hampson and Russell 1984). Refractionbased statics methods are based on absolute first-break arrival times and, in theory, are Chapter 3. STATICS CORRECTIONS Figure 3.3: Computation of delay times. Delay time at position i+n Di puted from picked travel times between i, i+n, and i+m. 40 +n can be com- able to estimate the long-period statics components. 3.2.1 Refraction statics algorithm A number of methods have been proposed to analyze first arrival times and calculate refraction statics corrections. INSIGHT REFSTAT2 module is a multi-layered approach which solves for the surface-consistent refraction statics using the time-term method of Farrell and Euwema (1984). Here we use a model to describe the details of the algorithm. In the application for our two datasets (Lines 88-14 and 88-17), the character of the first break picks along the entire line suggested that there was little first break energy associated with deep layers. 41 Chapter 3. STATICS CORRECTIONS The static corrections Cj applied to field data can be written Gi = hj/V - hi/V, (3.3) v where hj is the thickness of the weathering, considering as a single surface layer, V\ is the weathering velocity, V j is the refraction velocity of the layer under the weathering. 2 The parameter hj can be determined using the delay time Dj at station j (Gardner 1939, Hagedoorn 1959, Barry 1967). The delay time at station j is defined as the raypath time between the station and the refractor minus the travel time measured along the normal projection of this raypath onto the refractor (Figure 3.3). Dj = ja/Vt - abjV (3.4) 2i If we define sin I = V~i/V j, then Dj â€” hj co$l/v\. 2 Thus Cj can be rewritten as C^DjUV.j-V^HV.j + V,))^ (3.5) Implementation of the algorithm involves three steps: (1) Calculate refractor velocity V (XJ) 2 (Coppens 1985) Eachfirstbreak time is considered to consist of three terms: Ta = P * Xij + t d +t d si rj (3.6) where P is the reciprocal of the average refractor velocity, given by: - [(X-^r v^ x)dx]~ i p Xij is the offset. t d si and tfj are the delay times for shots and receivers, respectively. The algorithm uses the method of least squares to solve for the delay times. A linear combination of threefirstbreak traveltimes that correspond to certain shot-receiver pairs Chapter 3. STATICS CORRECTIONS 42 in the vicinity of Xj are used to calculate the velocity at the jth station. These are then averaged and smoothed to yield a slowly varying final result V {XJ). 2 (2) Determination of the layer thickness h(xj) The delay times are separated into high and low frequency spatial components by means of a nine-point median filtering operation followed by three-point smoothing. Then REFSTAT2 inverts equation 3.8 using the low-frequency component of the delay time to obtain the layer thickness h(xj). h(x )Jl t% = <- = \ s 1 J (^) Â« 2 W } (3.8) where t>i is user-specified replacement velocity. (3) Calculation of the statics correction Cj Once Hj and V j are obtained, the low frequency (long wavelength) of statics is 2 calculated by equation 3.3. 3.2.2 First-break picking Equation 3.4 requires Tij to be known. So we need to pick the first-break time Tij in common shot gathers. As in hand picking, identification offirstarrivals on common-shot gathers can be based on several criteria. The two criteria, sudden increase in energy on a trace and coherence of this increase in energy on adjacent traces, can be completed by one criterion, similarity of shape of thefirstarrivals on contiguous traces. In fact, when noise is present, it is often difficult to pick the realfirstbreak. So it is easier and more reliable to pick thefirstpositive or negative peak, depending on the polarity of the first arrival. Provided the shape of the first arrival is constant along the profile, the delay time of the peak relative to thefirstbreak can be estimated as a constant value to be subtracted from the picked times. Chapter 3. STATICS CORRECTIONS 43 The increase in the number of geophone groups in seismic acquisition during recent years and the requirement for accurate basic static corrections for high resolution records have made it necessary to develop sufficiently accurate automatic techniques for the determination of static corrections. Many papers have studied the problem of automatic picking of first arrivals (Chun and Jacewitz 1981, Hatherly 1982, Gelchinsky and Shtivelman 1983, Coppens 1985) and each technique has some inherent limitations. The INSIGHT's software module for the automatic picking of first arrivals also has its limitations. When the dip of an interface changes, the shape of the refracted wavelet changes, and the module has difficulty in tracing the first arrival, or the accuracy of the pick is not good. Application of the program is also sensitive to noise level; high signal-to-noise ratios provide better results. So we have to check each shot record for first-break picking. Sometimes, we have a "pitch out" i.e., a new event (with faster velocity) appearing ahead of the "old event". With a pitch out feature such as that shown in Figure 3.4, INSIGHT R E F S T A T 2 might have some problems (Kris Vasudevan, personal communication, 1995). So we treat the "second" break as a continuation of the "old event". Since this "anomalous" feature was not a common occurrence, our approach was to apply R E F STAT2 to the near offsets only , i.e., the near 40 stations. But when shot moves on, the "pinch out" will appear in the near offset. Any residual static effects can be recovered by application of residual static corrections. 3.2.3 Surface stacking charts and first break analysis Surface stacking charts (SSC), which are simply grids with shotpoint and geophone locations as the two axes, are used to display and analyze data-base information related to the acquisition and processing of reflection data. SSC can be used to analyze first-break picks and identify erroneous information at an early stage in processing (Milkereit 1989). SSC display can highlight the asymmetric and subsurface structure for the geography. Chapter 3. STATICS CORRECTIONS Shot 9450 0 2 0 44 File 7178 4 0 6 0 8 0 1 0 0 1 2 0 1 4 0 1 6 0 1 8 0 2 0 0 2 2 0 Figure 3.4: Example of one shot gather with "pitch out". One of the common parameters displayed in SSC is residual traveltime, defined here as: T[f' = T^-^--t -t d 'red d = e^ + 6e tj (3.9) where 7Â£ are the error-contaminated first break picks, V d is a " reducing velocity" which Te is a constant value that is representative of the regional bedrock refractor velocity, tf and tj are the delay time associated with positions i and j, e*J* are the systematic errors, and Seij are small random errors that embody any residual misfit due to the imperfect model specification and first-break measurement error. The residual traveltime is generally a smoothly varying quantity. Rapid variation indicates highly variable bedrock velocity or poorly-estimated delay times. The residual time errors can be displayed to evaluate the quality of picks and hence the quality of corrections. Through the residual traveltime analyses and based on our data, five types of errors contribute to residual traveltime anomalies that can be identified. They are Chapter 3. STATICS CORRECTIONS 45 polarity reversals, tuning errors, side-lobe picking errors, shot mispositioning, and geometry errors. The first three represent picking errors; the last two are associated with systematic errors. For each type of error, corrections can be made. For polarity reversals, the INSIGHT interactive editor "VAQ2" was used to multiply a factor of -1 on the trace amplitude, which simply reverses the polarity of that trace. For tuning errors, the anomaly is clearly characterized by the inverted " L " signature (Figure 3.5). The error was caused by an automatic picker cycle-skip over the correct first break traveltime when the site had a severe local statics problem. It was corrected by careful manual repicking. Side-lobe picking errors occurred when very strong first arrivals were recorded, resulting in high amplitude Klauder wavelet side lobes. This caused incorrect picking and occurred at several locations along the line. In each case, a careful manual repicking of the first break traveltimes rectified the situation. Identification of shot mispositioning was difficult due to the near-surface heterogeneity. We only corrected the first- break traveltimes, since the effects of this error are negligibly small for reflection arrival times except for the shallowest reflectors. To perform the correction, the first breaks for the whole shot gather are adjusted by adding a constant time to those traveltimes measured at forward receivers and subtracting a constant from traveltimes measured at reverse receivers. The amount of the correction was determined visually, by inspecting and comparing the residual difference plots before and after adjustment, so the correct adjustment restores the anomalous region to the background colors. For geometry errors, and after establishing that the residual traveltime anomalies were not caused by any other errors, the original survey file was checked. In the few cases where this occurred, the error was due to a data entry error for stations' coordinates by the field observers in the recording truck. Figure 3.5 shows a residual traveltime plot for shot gathers 570 to 800 with tuning error. We can also see the small master-slave discrepancy from it. There is a "recorder Chapter 3. STATICS CORRECTIONS 46 error" between the master and slave recorders. Because the slave recorder records the signal from the first 120 channels and the master records channels 121-240 (Figure 2.1), this error appears as a bulk time shift along the constant channel locus corresponding to channel 121. But as we will see, it was not a problem in our processing. In summary, SSC at the prestack stage provides a fast and convenient means for displaying, error-checking and analyzing large volumes of data such as trace edits, first-break traveltime picks and appraises the statics corrections. In comparison, the traditional traveltime plot (Figure 3.6), which is still commonly used in industry, is difficult to read, and it is not easy to diagnose the type of errors described above. Figure 3.5: Residual traveltime plot for shot gathers from 570 to 800. The solid arrows indicate the transition between the master and slave recorders. (A) The anomaly at station number 590 due to picking error is shown by the dashed lines. (B) After the corrections, there is no anomaly observed at station 590. Chapter 3. STATICS CORRECTIONS 48 Station No. Figure 3.6: Conventional INSIGHT first break traveltimes plot, x is pickedfirstbreak. â€¢ is calculated first break from refraction model. Thisfigureshows the difficulty in diagnosingfirst-breakand systematic errors from such plots. Chapter 4 DIP M O V E O U T AND VELOCITY ANALYSIS 4.1 Background and introduction The conventional normal moveout (NMO) is the basis for determining velocities from seismic data. The N M O and common-midpoint ( C M P ) stacking process enhances reflections having a particular moveout velocity, while more or less attenuating reflections having different slopes. Or in other words, this process also acts as a dip filter applied to the C M P stack. Fortunately, this dip-filtering action can be suppressed by applying a prestack process known as prestack partial migration, or dip-moveout. As the latter term implies, this process is a dip-dependent moveout correction that enables reflections from both horizontal and dipping reflectors to be stacked with the same N M O velocity. Stated another way, N M O velocities estimated from dip-moveout corrected seismograms Figure 4.1: Applying the law of cosines to triangle s'sr, one may express the travel time t from source s to receiver r in terms of zero-offset time to, half-offset h, velocity v, and dip 0. The result is equation 4.1 in the text, the familiar dip-corrected N M O equation. 49 Chapter 4. DIP MOVEOUT AND VELOCITY ANALYSIS 50 are independent of the dips of surface reflectors. For a single dipping reflector, we must consider the NMO equation, which relates the two-way traveltime to the characteristics of the medium. ?/, ?/~x Ah cos0 t\h) = * (0) + ^ . 2 ,, â€ž. 2 (4.1) 2 where 0 is reflector dip, v is true velocity of the medium above the reflector, and h is the half offset (Figure 4.1). If the moveout term is divided into two parts, then â€¢>,,^ ->,~x 4/i t (h) = t 0) + â€” 2 2 4h sin0 2 2 â€” . 2 The first part of the moveout (4.2) is associated with the zero-dip normal moveout (NMO), while the second part is associated with dip moveout (DMO) since it is related to reflector dip. If we define NMO time (not dip-corrected) by then we must have: t, = Â« - (4.4) Let p(t,y,h) denote the seismogram recorded at time t, midpoint y, half-offset h. We can rewrite equation 4.3 to define NMO as: NMO (t ,y,h) Pn = p(y/t + 4h /v ,y,h) 2 n 2 (4.5) 2 n and use equation 4.4 to define DMO as: DMO (t , y, h) = (y/t - Ah sin0 /v ,y, h) 2 Po 0 2 2 (4.6) 2 Pn If we think of NMO as the transformation from recording time to NMO time t , then n DMO is the transformation from NMO time t to zero offset time to. The above equations n imply that we can apply the NMO correction using the medium velocity, and then apply Chapter 4. DIP MOVEOUT AND VELOCITY ANALYSIS 51 the D M O correction. Unlike the N M O term, which is implemented in the C M P domain, the D M O term needs to be implemented in a domain in which dips can be recognized, such as the common-offset domain. From equation 4.2, we can readily assess the properties of the D M O term 4fc2 ^" fl2 â€¢ First, it has no effect on zero-offset data (h=0), regardless of dip. Second, steeper dipping structures require larger corrections. This implies that the shallower the event, the more significant is this term. Also, since lower velocities generally are found at shallower depth, the dip moveout term is more significant for shallower events where velocity is slower. Stacking the N M O - and DMO-corrected gathers yields a section that more closely represents the zero-offset section than the stacked section without the D M O correction. Based on studies by Yilmaz (1987), N M O + D M O + stack + time migration after stack is approximately equivalent to full time migration before stack. Also, since D M O correction is a migration-like process, it causes the energy to move from a C M P gather to neighboring gathers in the updip direction. 4.2 4.2.1 Implementation of DMO Partial stack The D M O algorithm operates separately on each common offset gather. A complication arises because of the irregular offset distribution associated with a crooked line geometry. We need to group the traces according to some measure of nominal offset. The procedure we follow is analogous to the approach we adopted in performing the C M P binning. "Common offset binning" which sums traces within each offset range effectively decreases the number of common offset planes and increases the fold within each common offset gather. Since this process involves a summation or "stack" in the C M P domain, the common offset binning procedure is referred to as the "partial stack". Chapter 4. DIP MOVEOUT AND VELOCITY ANALYSIS 52 The partial stack algorithm we used is a revised version of the INSIGHT " P S T K " module (Weizhong Wang, personal communication, 1994). Two improvements are achieved (Figure 4.2). First, we can infill an empty bin with a trace formed by taking the mean of the two traces that occupy the surrounding bins. This is helpful to avoid the energy smearing property, which is effected by seeding a dead trace when no trace occupies a particular bin. Second, since D M O correction does not depend on the polarity of the offset, only on the absolute value, we can use the absolute value of offset. This increases the number of traces within each common-offset bin and reduces the number of commonoffset gathers, thus minimizing the run time significantly. Along the noisy portions of the line, the increase of the number of traces in each bin helped improve the signal-to-noise ratios. 4.2.2 DMO algorithm The algorithm that was applied to effect prestack partial migration of the data and reduce the dependence of the N M O velocity on reflection dip operates in the frequencywavenumber domain. It is based on the Hale method (Hale 1984), which is computationally intensive. From Figure 4.1, we can see t v = 2(y - )sin9 0 yi (4.7) where yi denotes the intersection (not shown) of the reflector and the surface. For the subsurface model of Figure 4.1 and from equation 4.7, the slope of a single reflector is Ai 2sin0 . , -zâ€” = l -Â°J Ay v as illustrated by the zero offset section in Figure 4.3. Then the D M O transformation of 0 4 Chapter 4. DIP MOVEOUT AND VELOCITY ANALYSIS zoom box h (c) 53 -haif-oteet t Figure 4.2: (a) a common midpoint plane; (b) expanded view of (a) with definition of the offset bin; (c) result of applying partial stack to (b). Traces C, D and E are stacked to form trace G . Traces A , B and F are stacked to form trace I. Trace H is the mean of trace G and I. Chapter 4. DIP MOVEOUT AND VELOCITY ANALYSIS 54 Po('o.y. 0 / Figure 4.3: Zero-offset section for the experiment of Figure 4:1. The slope of the reflection is At /Ay = 2sind/v, where 6 is the dip of the reflector and v is the velocity. 0 equation 4.6 may be expressed as DMO (t , y, h) = (y/t 2 Po 0 Pn 0 - (Ato/Ayfh*, y, h) (4.9) Note that D M O via equation 4.9 does not require knowledge of the velocity v or dip 6. Because these parameters are typically unknown, equation 4.9 might seem to be preferred to equation 4.6. As a practical method of applying D M O , however, equation 4.9 is not very practical for two reasons. First, the slope Ato/Ay required in equation 4.9 to be used must be measured on a zero-offset section, and zero-offset sections are typically not recorded. A second and more significant problem with equation 4.9 is that more than one slope may exist at a given (to,y). The conventional process of stacking after dip-corrected N M O (using either equation 4.6 or equation 4.9 ) must enhance one slope at the expense of attenuating the other. Hale's method (Hale, 1984) solves the problem of how to treat conflicting slopes. Equation 4.9 implies that a different D M O correction is needed for each different zerooffset slope Ato/Ay. D M O , as the name implies, is a dip-dependent process. The two- dimensional (2-D) Fourier transform (over t and y) provides a particularly useful domain 0 in which to apply a dip-dependent process, because all events having a particular slope 55 Chapter 4. DIP MOVEOUT AND VELOCITY ANALYSIS Pn(tn, y, h) 1-DFFT Pn(tn, k, h) integrate over tn for all k andwo Po(wo, k, h) 2-D Inverse FFT Po(to,y,h) Figure 4.4: The implementation of D M O algorithm by Fourier transform in the (io,y)-domain transform to a single radial Une in the (u>,k) -domain given by: k Ar-o 2sind LL> Ay - = â€” = . . (4.10) v where ui is the transform variable associated with the two-way zero-offset time ioFigure 4.4 illustrates the most straightforward implementation of Hale's method. Because h appears only as a constant parameter in the D M O flow, the data should be sorted into constant-offset sections prior to their use. D M O is then most naturally applied to each constant-offset section separately. Chapter 4. DIP MOVEOUT AND VELOCITY ANALYSIS 56 For the common offset case both input and output are in midpoint, y, spatial coordinates and the general input/output relationships reduces to INPUT : P {t ,y) n OUTPUT : P (t ,y), n 0 0 (4.11) The coordinate relationships, t (t ), between common offset data and zero offset data is 0 n the dip-corrected NMO equation (Hale 1984) given by h?k 2 to = M l + 7 T - | ] 1/2 = *Â«A (4.12) The general DMO formula reduces to Po(io, y> ) = ^j h dw e- iuoto 0 J dke P (u , k, h). (4.13) iky 0 0 which is the result given by Hale (1984; see Appendix A). Figure 4.4 provides a method for applying DMO correction, not just for one slope, but for each slope At /Ay = k/ui (equation 4.10) in the zero-offset section. It proves what 0 we said before, that in the limit of zero offset (h/t â€” > 0), and zero slope (&/u> â€”â€¢ 0), n 0 DMO does nothing (equation 4.12). DMO is most significant for early reflections at large offsets (h/t n oo) and steep slopes (/e/u> â€”* oo). 0 DMO processing is time-consuming. For the sake of computational speed, we used the log-stretch DMO algorithm. It applies a logarithmic stretch to the time coordinate, proceeds with DMO correction in the frequency-wavenumber domain, and then removes the logarithmic stretch from the DMO-corrected data. The procedure is based on Hale's DMO formula as derived by Liner and Bleisten (1988). It gives the correct impulse response geometry. The algorithm produces results similar to those obtained by Hale's DMO module, but at a substantially reduced running rate (to about one third). On a Sun SPARCstation 10, 250 CDPs took about 3 hours to process. Chapter 4. DIP MOVEOUT AND VELOCITY ANALYSIS 57 We should note that the computed po(t ,y,h) should be interpreted as zero-offset 0 sections only in the sense that their reflection times are the same as those in the true zerooffset section p(t , y,h = 0). This qualification bears repeating because replacing At /Ay 0 0 with k/u> is valid only with respect to the arrival times of reflections. A more rigorous (and considerably more complicated) analysis of the transformation from nonzero to zero offset must be based on wave theory, whereas the results in this study are from a geometrical analysis of ray-paths. 4.3 Velocity analysis after DMO The C M P stack can be improved by using the velocities repicked after D M O correction (Yilmaz 1987; Hale 1984). To test this idea, consider the following procedure. First, apply inverse N M O correction to the gathers with the velocity that was used in the first N M O correction step. Then, assuming we pick the correct velocity function following D M O correction, use it on the second N M O correction. From such experiments, we arrive at the flowchart in Figure 4.5 for velocity analysis with D M O processing. The most important reason to obtain a reliable velocity function is to generate the best quality stack. There are two common methods for velocity analysis. One involves generating constant velocity scans of C M P gathers; stacking velocities are estimated from the data stacked with a range of constant velocities on the basis of stacked event amplitude and continuity. Another method is the velocity spectrum which is based on the crosscorrelation of the traces in a C M P gather, and not on lateral continuity of the stacked events. Some additional aspects of velocity analysis as used in this study should be pointed out: (1) From equation (4.1), an optimum stack would be obtained if the N M O correction Chapter 4. DIP MOVEOUT AND VELOCITY ANALYSIS (1) Data (2) Velocity Analysis (3) N M O Correction Using Flat-event Velocities (4) Data Scaling (5) D M O Correction (6) Remove Data Scaling T (7) Inverse N M O Correction Using Velocities from (3) T (8) Velocity Analysis Figure 4.5: Processing stream of velocity analysis with DMO. 58 Chapter 4. DIP MOVEOUT AND VELOCITY Time (s) Velocity from N M O (m/s) ANALYSIS 59 Velocity from D M O (m/s) C D P 2420 2.00 2.90 5.30 5938. 6146. 6458. 5796. 6013. 6392. C D P 4055 1.48 2.10 3.91 6042. 6146. 6250. 1.37 2.52 3.46 5750. 6083. 6250. 5904. 6013. 6175. C D P 8180 5688. 6013. 6121. Table 4.1: Comparison of velocities from N M O and D M O processing. were performed using a "stacking velocity" V , which is So V should be higher than the medium velocity V . After D M O , we should obtain a lower medium velocity by velocity analysis than before D M O correction. Table 4.1 shows a comparison of picked velocities with and without D M O correction. (2) Of the many coherency measures proposed over the last twenty years for velocity analysis, the "semblance" criterion (Neidell and Taner 1971) is the one most commonly used by the exploration industry. However, when data are of lower quality or structural complexity and velocity inhomogeneities wrap the hyperbolas, then constant velocity stacks are often the best way to go. Since the data used in this study are of the latter type, the constant velocity analysis procedure was used instead of semblance velocity analysis. (4) The "energy smearing" action of D M O , which moves energy from a single midpoint to neighboring ones on a common offset gather, is a common problem. We anticipate such "smearing" due to the migration-like action of D M O . Also, the D M O process, which maps Chapter 4. DIP MOVEOUT AND VELOCITY ANALYSIS 60 data from x-t space to the f-k domain (by a 2-D Fourier transform), may be susceptible to aliasing problems. Fortunately, we do not need to exploit the anti-aHasing feature of DMO in our present study on Lines 88-14 and 88-17. Based on a spatial sampling interval of 25 m, the threshold aliasing frequency is 50 Hz for Lines 88-14 and 88-17. Given the characteristics of the vibroseis sweep (frequency range 10-56 Hz), it is unlikely that DMO-induced aliasing effects are significant for our study. The DMO process could degrade the quality of stack along crooked parts of thefines.As shown in Figure 4.5, we scale the data before DMO processing and remove the scaling after DMO. With DMO processing, a better stack was obtained and velocities were dip-corrected, resulting in a better migration after stack. The imaging quality is similar to full migration before stack. Significant improvement is seen when shot gathers are stacked using DMO processed velocities (Figure 4.6). The example shows a structure stack at station number between 8810 and 10070 of Line 88-14, where (A) used velocities derived from regular NMO - corrected analysis, (B) used velocities derived from NMO+DMO - corrected analysis. The mid-crustal dipping reflector has better continuity and more distinctive amplitude in (B). Chapter 4. DIP MOVEOUT AND VELOCITY ANALYSIS Figure 4.6: Structure stack at station number between 8810 and 10070 of Line 88-14. (A) Using velocities derived from regular NMO - corrected analysis. (B) Using velocities derived from NMO+DMO - corrected analysis. Chapter 5 CROSSDIP 5.1 CORRECTION Introduction Analysis procedures described in previous chapters assumed that the survey line is straight, receiver groups and shotpoints are uniformly spaced, and group spacing and shotpoint spacing are properly related. For crooked survey lines (Figure 5.1), these assumptions are not realized; there is usually a scatter of C M P points about the line of acquisition. For dipping subsurface structures, time errors that impair the stacking process are associated with this scatter. Corrections applied to data traces to compensate for these time anomalies are called "crossdip corrections", for which crossdip refers to component of dip perpendicular to the local direction of the survey line (Figure 5.2). If the slalom line were strictly a dip line, the dynamic N M O corrections would be approximately hyperbolic and, consequently, reflection events would stack properly. If however, subsurface interfaces have components of crossdip, the time of a reflection from Figure 5.1: Plan view showing midpoint scattering along a portion of a crooked survey line. Small dots mean midpoints, and black line stands for output profile. Y is crossline offset of the midpoint from the effective profile. One C M P bin is indicated by the dashed pair of lines. 62 Chapter 5. CROSSDIP CORRECTION 63 a given interface will vary according to shotpoint-receiver displacement in the direction perpendicular to the profile line. Traces whose midpoints are relatively downdip will show longer reflection time than traces with midpoints updip. Such variations in reflection times cannot be corrected by residual statics. Consequently, when crossdip is not included explicitly as a component in the residual statics solution, reflection quality on stacked sections must suffer. 5.2 5.2.1 Theory 3-D C M P bin traveltime equation After CMP binning, or regrouping of the traces, the midpoints of traces associated with an output trace gather remain scattered. How does this midpoint scattering affect the quality of the CMP stacks when three-dimensional dipping structures are involved? First, we consider a three-dimensional model consisting of a single dipping reflector with constant medium velocity V. Perz (1993) gave the full description of the reflection equation for a crooked line CMP gather. With reference to Figure 5.2, we define 7 = VT + tan -! + tan C (5.1) 8 â€” tanlcosp + tanCsinp (5-2) 2 2 where I and C are the angles of inline dip and crossdip, respectively, and p is the sourcereceiver azimuth. If we define the "crossdip slowness" p as: ItanC P= . s 77~Â» I-) 5 3 then the reflection traveltime t with an arbitrary midpoint is given by t' = (tc+ Y) 2 2 P + ~ ( l - ( S - Y ) . (5.4) Figure 5.2: (a) Plan view. The slalom line is coincident with the x-axis and the CDP bin centre is located at the origin, (b) Three-dimensional view of the earth model. V is layer velocity. I is the angle of the inline dip. C is crossdip. Chapter 5. CROSSDIP CORRECTION 65 where B is the source-receiver offset, t is the zero-offset two-way traveltime measured at c the bin centre, and Y is the distance between CMP centre and individual CMP. For the general 3-D situation (equation 5.4), we note that all dependence on the source-receiver azimuth is confined to the second term in the equation, the "moveout" term. Considering the limiting case of zero crossdip, for a source-receiver azimuth along profile, we can rewrite equation 5.4 as: t* = tl + f ^ c o , 2 / ) (5.5) which is the standard dip-moveout equation. Considering the limiting case of zero inline dip, we can rewrite equation 5.4 as: t" = (t c + Yr + ^ P (5.6) where p simplifies to Pcrosa â€” y I) so after NMO correction, t = t + pY 0 5.2.2 c (5.8) Simplified 2-D case The amount of traveltime difference attributable only to crossdip (no inline dip) can be quantified easily from Figure 5.3. It shows a cross-section taken perpendicular to a crooked survey line at CDP k. Assume that NMO has been applied so that each unstacked trace simulates a coincident shotpoint-receiver experiment located at the shotpointreceiver midpoint corresponding to that trace. If we assume that the single corresponding reflector is planar, the difference in reflection time 8tij between a trace with its midpoint on the effective profilefineand one Chapter 5. CROSSDIP CORRECTION 66 2 _AT _ 2 sinY _ p Figure 5.3: Cross-section showing a crossdipping reflector. The effective seismic profile, perpendicular to the plane of the figure, is indicated by the dot. displaced a distance Yij (defined as transverse offset) from it is given by f>Uj = â€”Tp^-Yij where D is constant for a given k, j = k k = E> Yij (5.9) k crossdip angle at C D P k, V =velocity at C D P k k, i=shot number, j=receiver number of the trace under consideration. A straightforward geometrical argument shows that St{j is linearly related to Y^. For the "crossdip parameter" p, t = t + pY, (5.10) =^P. (5.11) e where ? Equation 5.10 suggests that the crossdip correction might be performed by stacking NMO/DMO-corrected C M P gathers along linear trajectories in (Y-t) space. Such a procedure is readily implemented using "SLNR" and "DCOR" in the existing INSIGHT software. Because of its ease of implementation, a crossdip correction based on equation 5.10 is attractive from a practical viewpoint. For the case of a multilayer model in which all interfaces are planar, Larner et al. (1979) showed that equation 5.10 is still valid. Consider normal-incidence rays from Chapter 5. CROSSDIP CORRECTION 67 two neighboring points on the earth's surface. So long as layer interfaces are planar, Snell's law dictates that the raypaths will be parallel in every layer. Also by Snell's law, the difference in reflection times will equal twice the time that it takes a wavefront propagating along the rays to travel from one surface point to the other. Because the rays are always parallel, the reflection time difference is linearly proportional to the distance between the two surface points. So for a multilayer model, 7^ is the emergence angle of the reflected rays and V*. is the velocity of the uppermost layer. Comparing the 3-D and 2-D cases (from equation 5.4, 5.6 and 5.10), we note the difference is the second term in (5.4), the "moveout" term. Therefore 3-D effects on traveltime moveout may be described in terms of increased "effective" moveout velocity. We see here that the three-dimensionality serves to further reduce the amount of traveltime moveout. We note actual P-wave velocities of the mid to lower crust are sufficiently high that moveout effects are not significant for later events. On the basis of these observations, we feel that we are justified in correcting crossdip effects, at least tofirstorder, by application of the slant stack. 5.3 Crossdip correction processing stream The crossdip correction processing stream is shown in Figure 5.4. The method used here assumes that (1) the travel time correction is independent of the azimuth of the CDP bin line, and (2) the residual moveout is tofirstorder the product of slowness and transverse offset. The basic method consists of applying NMO using the "normal" offset in the "x" direction, then do slant stacks for different slownesses using the transverse offset to indicate the distance between traces, andfinallystack using the slowness giving the highest amplitude. The transverse offset is calculated as the distance between the CMP and the CMP bin center (Figure 5.3). The direction of the line is only used to calculate Chapter 5. CROSSDIP CORRECTION YOFF 68 calculate transverse offsets r NMO NMO i correction using velocity obtained after DMO \ MHED c o p y header 19 into header 53 forward and inverse slant stack transformation for a series o f slownesses pick the best slowness D i s p l a y and p i c k in x and T and create a file analogous to velocity file DCOR apply crossdip correction to N M O e d CDP gather Figure 5.4: Processing stream of crossdip correction. the sign of the offset. For this reason, bins should not be too wide. The direction of the CDP bin line is defined only once. Thus, CDP lines should hence not change direction sharply. The "constant slowness stack" (CSS) is a straightforward method to obtain the crossdip parameter. The procedure is implemented using the SLNR module, which performs the forward slant-stack transformation by computing tau-p transforms for a series of Chapter 5. CROSSDIP CORRECTION 69 slownesses and stacking the gather for each slowness (Sacchi 1996, Ph.D. thesis); following which it performs the inverse slant-stack transformation with the application of the "RHO" filter to reconstruct t-x seismograms. The integration component of slant stack procedure enhances low frequencies. The "RHO" filter, i.e., |u;| filter (Claerbout 1985), pushes them back to their appropriate level such that the slant stacking does not cause the power spectrum of the data to change. Analogous to direct velocity picking from constant velocity stacks, the best slownesses in x and t are picked in the slant stack data and a file for running crossdip corrections is created. 5.4 Example of application of crossdip correction The crossdip correction significantly enhanced the continuity of reflections in many areas of the upper 5 s of the seismic section. Based on tests during this research, the correction is not very effective beyond 6 s. Figure 5.5 displays a portion of the stacked section in the middle of Line 88-17. Compared with (A), (B) shows improvement of the seismic image, which has better continuity and more distinctive amplitude for reflectors around 1.6 s, thereby making interpretation more reliable. Application of crossdip correction to the data set provided an overall improved image, particular for Line 88-17 along which geology is complex (Figure 1.3). 5.5 Spread length truncation Another approach used against strong heterogeneity degrading the final stack is the creation of an alternate stack using a truncated spread. The recording spread for the Southern Cordillera data is 4.0 km - 0 - 8.0 km. In order to reduce possible smearing in the upper few seconds, an alternate stack is created with truncated spread length. Figure 5.5: Figure shows the same section after application of crossdip correction. For the event at 1.6 s, the improvement in reflector continuity is quite impressive. Chapter 5. CROSSDIP CORRECTION 71 Offset-limited (3.5 km - 0 - 3.5 km) stacks controlled the negative contributions of extensive lateral variations along the travel paths, reduced surface noise levels and enhanced continuity of reflections in many areas of the upper 6 s. Figure 5.6 shows the truncated stack in western part of Line 88-17. Chapter 5. CROSSDIP CORRECTION 72 Figure 5.6: (A) Use of the full offset spread length (4.0 km - 0 - 8.0 km). (B) Use of a truncated spread length (3.5 km - 0 - 3.5 km). Note the better definition of coherent energy on the section in (B). Chapter 6 P O S T - S T A C K P R O C E S S I N G A N D F I N A L SEISMIC S E C T I O N S In this chapter, we describe the final stages of post-stack processing - random noise attenuation, data merging, coherency filtering and migration - and present the final structural and migrated stacks. 6.1 Application of the partial Karhunen-Loeve transform to suppress random noise 6.1.1 Introduction The Karhunen-Loeve (K-L) transform is a signal enhancement technique used in digital image processing (Oppenheim 1978). It is also known as the principal component transformation or the eigenvector transformation. Signal/noise ratio improvement is one of its applications. The K-L transform decomposes an image into components (principal components) that are ordered on the basis of spatial correlation. When a seismic section is considered to be an image of the subsurface (which we hope it is), the usefulness of this property is immediately apparent: spatially uncorrelated parts can be removed leaving behind a clear and coherent image. The K-L transform is an effective technique for suppressing spatially uncorrelated noise. Oppenheim (1978) and Marchisio et al. (1988) showed how the full K-L transform is used to build filters that enhance not only horizontal events but also dipping events (which can occur in the same image). Recently, a technique that combines zero-lag K-L transform and linear moveout to make the K-L 73 Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS 74 transform fast enough for typical stacked seismic data was discussed by Al-yahya (1991). This technique represents partial (or limited) filtering. The dip-limited K-L transform can be thought of as an extension of the procedure described by Jones and Levy (1987) to filter data dominated by dipping events. 6.1.2 Basic theory A seismic section which consists of M traces with N points per trace may be viewed as a data matrix X where each element Xij represents the i th X = Xn X12 ... Xiff X X2 â€¢â€¢â€¢ XN %M2 â€¢â€¢â€¢ XMN 21 2 XMl point of the j th trace. (6.1) 2 We can look at X as being made up of columns, X = C\ C 2 (6.2) N c_ ... Marchisio et al. (1988) start by building the covariance matrix C of X from the submat rices Cij, where i\Xj\ Xi Xj Xi Xj\ X{ Xj XiffXji XipfXj x CiCj â€” C, 2 2 Xu Xjjf 2 2 X{ XjN 2 2 2 ... (6.3) XiffXjjf and Xki is the 1th sample of the kth trace. The sums of the diagonals in any CiCj gives the values of the cross-covariance function between the ith trace and jth trace at all lags. Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS 75 Let C be the (M x JV) covariance matrix of the input data X computed as C â€” XX T where X is an (M x JV) matrix. Then we may write C = ITEU (6.4) T where E is an (M x M) diagonal matrix composed of the eigenvalues of C and U is an (M x M) matrix, the columns of which are composed of the eigenvectors of XX . T The Karhunen-Loeve or principal component matrix of X is given by Y, an (M x N) matrix, where Y = UX (6.5) T Because the singular value decomposition (SVD) of X is written in matrix form as X = UAV (6.6) T where V is an (JV x JV) matrix, the columns of which are composed of the eigenvectors of X X, T and A is an (M x JV) diagonal matrix with entries A i > A > A > ... > A M , where A ; 2 3 is the ith singular valve of X, we obtain Y = U UAV T T = AV T (6.7) The principal components contained in the matrix Y may be viewed as the inner product of the eigenvectors of XX T with the data. Since X may be reconstructed from the principal component matrix Y by the inverse K-L transform X = UY (6.8) then using equation (6.7), X can be written as X = UAV T (6.9) Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS 76 Equation (6.9) is identical with equation (6.4), showing that X can be reconstructed through K-L transform or through SVD. Providing we are considering a multivariate stochastic process, the SVD and K-L transformation are computationally equivalent. From a stochastic point of view, we now have one realization, X, of a two-dimensional random process. The K-L transformation and SVD ( Lanczos 1961) of X will be the same only if we assume the separability of the covariance matrix of the process into a product of the covariance matrices of the rows and columns of the image, and if these covariance estimates are computed from the one realization which is our image (Gerbrands 1981). Here, eigenimage decomposition, which depends on a deterministic decomposition using SVD, is perfectly descriptive. It should also be pointed out that in the case where the SVD and K-L decomposition are equivalent, as in the situation described above, the K-L transformation is generally referred to as the zero-lag K-L transformation since it is computed from only the zero-lag covariance matrix. Applications of the full K-L transformation to geophysical image enhancement has been discussed by Marchisio et. al (1988) who noted that since the covariance matrix contains all possible lags, for an N x N matrix, the size of the covariance matrix C is N xN . 2 2 Since the K-L transform requires eigenvector analysis, then, if we make the practical limit for the size of the matrix for this analysis to be, say 400x400, the largest seismic section we can consider forfilteringwith the full K-L transform is 20 x 20, which is extremely small. A reasonable computational threshold for a SVD on an ordinary mainframe is around 100x100 matrix, which puts a Hmit of 10x10 on the size of the original data window. So the implementation of the full K-L transformation requires the use of a supercomputer. The full K-L transform is therefore prohibitively expensive to use for filtering seismic data. Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS 77 6.1.3 Computing the filtered image by partial K - L transform As a compromise between the full K-L transform (Marchisio et. al 1988), and the zero-lag K-L transform (Jones and Levy 1987), which enhances horizontal events only, we used the partial K-L transform. In this approach, the K-L transform is performed for selected dips, and the results of the K-L filtering of all dips are added together. For each dip, linear moveout is applied to the data, i.e. d(x,t)' = d(x,t + x) P (6.10) where p is the dip in samples/trace which does not have to be an integer. After the linear moveout, the data are K-L filtered using the zero-lag crosscorrelations (computationally equivalent to the SVD method), and then inverse linear moveout is applied to restore the data to their original orientation, thus d(x, t) = d(x, t â€” px)' (611) Since the K-L filtering after linear moveout uses zero-lag cross-correlation only, it is very fast; we can filter data whose size is, say 400x200 (i.e. 200 CDPs and 1.6 s, at a sampling rate of 4 ms) in a reasonable CPU time (a few seconds). Now we turn to the problem of the actual computation of the desired image at each dip. Since we reconstruct X from a few eigenimages, the band-pass matrix or SVD-filtered data XBP is given by X = UBPABPVJP BP (6.12) where UBP, VBP and KBP represent U, V and A with first (p-1) columns and the last (r-q) columns set to zero. We compute Ug X in the case that M < N as: P Ul X P = Ul UAV T P (6.13) Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS 78 Owing to the orthonormality of the eigenvectors it follows that UlpUK = A BP (6.14) and substituting equation (6.14) into equation (6.13) we obtain Ul X =A V T P BP (6.15) Owing to the zero values along the diagonal of A^p it is clear that A V = AspVgp T BP (6.16) Using the above expression in equation (6.15), U X = AapVgp T BP (6.17) Finally, substituting equation (6.17) into equation (6.12) we obtain X p = UspUlpX B (6.18) In a similar fashion, when M > JV, we could obtain X BP = XVspVgp (6.19) The computation of SVD filtered data X p is more efficiently performed by computing B the principal eigenvectors of either XX or X X depending on whether M < N or T T M > N, respectively. When M < N, use equation (6.18). When M > N use equation (6.19). The above development may be repeated for the complex case by replacing T, the matrix transpose, with H, the complex transpose. The flowchart of the code that we prepared is shown in Figure 6.1 There are several points to be noted with respect to the computation: (1) To avoid losing the data at the edges after linear moveout, I need to zero-pad it, Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS 79 which increases the cost of the K-L transform. (2) The dip p (samples per trace) does not have to be an integer, and interpolation can be used to get better dip resolution. (3) In general, there will be several dip components in the seismic data and their dips will not be known exactly. For this reason, we should include a range of dips. (4) . When a dip is included, the signal at that dip is enhanced and some noise may also be added. This trade-off between noise and signal depends on the data. 6.1.4 Summary and results The partial K-L transform is made feasible by limiting the number of dips included in the covariance matrix of the full K-L transform. There are several advantages. First, each dip is included by a linear moveout of the original image. After this linear moveout, the image can be decomposed by the zero-lag cross-correlation (computationally equivalent to the SVD method), which is very fast. Secondly, by including several dips in the reconstruction, we need only a few image components from each dip because coherent events at those dips appear in the first few singular values. Thirdly, by judicious selection of the dips, dip filtering can be carried out simultaneously with the K-L filtering. This techniquefiltersthe data for limited dips. For each dip, linear moveout is applied to the seismic section so that events with this dip are made flat. After linear moveout, zero-lag K-Lfilteringis applied followed by inverse linear moveout. At the end, the results from all dips are added to form thefinalfiltereddata. The partial K-L transform was applied to suppress random noise in stacked sections of lines 88-14 and 17. Figure 6.2 shows a stacked section of Line 88-14 with prominent dipping reflector and background random noise. Figure 6.3 shows the reconstruction of the image in Figure 6.2 following partial K-Lfiltering.This image used thefirst5 components of the SVD process and linear moveout factor p sets from 0.1 to 0.7 units. The Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS Read Data X: M X N (M traces, N points/trace) Linear Moveout X , t + p ix Â» > x , t' M>N M<N c=x x T C = X X T C = U A C = V U XBP=UBPUBPX A V XBP=XVBPVBP \ Linear Moveout x , t ' - p i x Â» > x , t â€¢ O u t p u t f o r e a c h D a t a p i Xi: M X N (M traces, N points/trace) l<W I= w â€¢ Output Data X=Xl+...+Xi...+Xw (M traces, N points/trace) Figure 6.1: Flowchart of application of the partial K - L transform. 80 Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS 81 background random noise has been reduced while retaining the strong dipping reflector. Figure 6.3 also illustrate a potential problem, a tendency to enhance weak alignments within the dip-limited window that may not be desired signal. Thus caution must be exercised in the application of the procedure. The reason for this good reconstruction can be seen in the singular value spectrum (Figure 6.4, in which the signal and noise are well separated). By choosing only the first 5 components, most of the random noise, associated with the higher components, is removed. Indeed, the number of eigenvectors plays an important role in the reconstruction of the image. For example, Figure 6.5 shows the image reconstruction using 15 image components instead of 5. The overall background noise level is much higher than in Figure 6.3. In some instances, however, the addition of more eigenvectors is necessary for an accurate reconstruction of the signal. This K-L application works well for seismic sections in the case of horizontal events, a single dip or various dips. The INSIGHT software also has package module "KLSTAK" which performs signal-noise-ratio enhancement via the zero-lag K-L transform. Unfortunately, its application did not create good results. One reason is that its approach of zero-lag K-L transform couldn't handle seismic sections with various dipping reflectors. 6.2 Comparison of reprocessed data with contract processed data The reprocessing at UBC emphasized images of upper crustal features; contract processing emphasized images of the whole crust. Better definition of coherent energy in the upper 6 s is shown on the UBC sections. Refraction statics, cross-dip corrections accounting for 3-D effects in the vicinity of the line and use of a truncated spread length (3.5 km - 0 - 3.5 km vs. 4.0 km - 0 - 8.0 km) were the procedures which contributed most to the improved images. Specific examples of the improvements achieved are shown in Figure 6.6 and 6.7. Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS 1 21 Figure 6.2: Part of the stacked seismic section, of Line 88-14 showing a dipping reflector and background random noise. 6.3 Coherencyfilteringand final stack presentation Following application of the procedures discussed previously in this thesis, structural stack sections for the upper 6 s of Line 88-17 (Figure 6.8) and 7 s of Line 88-14 (Figure 6.9) were generated. The upper 6 s of thefinalLine 88-17 zero-offset stack section with the set of stacking velocity profile (both RMS and interval velocities) and the station elevations are shown in Figure 6.10. The upper 7 s of the final Line 88-14 zero-offset stacked section is shown in Figure 6.11. The horizontal-to-vertical aspect ratio is 1:1 assuming a constant velocity medium of 6000 m/s. The relatively poor S/N ratio presents a problem from the viewpoint of presentation 82 Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS 1 21 Figure 6.3: The reconstruction of the image in Fig. 6.2 following partial K-L filtering results and subsequent interpretation. The compressed spatial scale at which the sections are displayed often exceeds the resolution capability of the human eye, rendering detection of weak events difficult, if not impossible. Automated coherency filtering provides an objective means of isolating coherent, linear signal from a background noise field, thereby improving the display. In this study, the coherency filtering algorithm of Lane (1991) was used. For a given input point (xj,tj), the algorithm computes semblances over a lateral trace window for a range of apparent slownesses (i.e., time dips). A coherency measure equal to the semblance value raised to a user-specified exponent is then computed for each time dip. The trace amplitude at the point is smoothed by averaging the amplitudes over 83 Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS 84 the variation of singular values with no. 1.5 0.5 t- 10 20 singular value number 30 Figure 6.4: Singular value spectrum of the image in Fig. 6.2. Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS Figure 6.5: The image reconstruction of the image in Fig.6.2 using 15 images, uates less noise than Fig. 6.3. Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS 86 Figure 6.6: Comparison of unmigrated seismic section along eastern part of Line 88-17 (following page). Special aspects of processing at U B C included: (1) refraction statics (2) cross-dip correction, accounting for 3-D effects in the vicinity of the line. Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS 87 Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS Figure 6.7: Comparison of unmigrated seismic section along western part of Line 88-17 (next two pages). Special aspects of processing at UBC included: (1) refraction statics (2) cross-dip correction, accounting for 3-d effects in the vicinity of the line and (3) use of a truncated spread length (3.5 km - 0 - 3.5 km vs. 4.0 km - 0 - 8.0 km). Note the better definition of coherent energy on the UBC section. 88 Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC TWO-WAY T I M E ( S ) SECTIONS Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS 91 the trace window along the direction of maximum coherency. Finally, multiplication of this smoothed result by the maximum value of the coherency measure yields the output trace amplitude at (XJ, tj). Additional noise rejection is accomplished by the introduction of a plotting bias. A constant negative value, specified as a percentage of the average absolute amplitude computed over the entire dataset, is added to all trace amplitudes. When the coherency filtered dataset is plotted using the variable area display only (i.e., no wiggle trace), only relatively large positive values contribute to the plot. Because linear events significantly shorter than the window length are attenuated, a lateral window equal to the shortest length of coherent, linear energy observable on the final zero-offset stacked section (15 traces) was specified for the display. Semblances were computed for 15 time dips between Â±0.33 ms/m, these bounds corresponding to the maximum time dips observable on the final zero-offset section (Figure 6.8 and 6.9). Based on a visual inspection of the output of several algorithm test runs, user-selected parameters chosen were a semblance exponent of 1.0 (i.e., the coherency measure was set equal to the semblance) and a plotting bias of -200 percent. Larger magnitudes of plotting bias tended to obscure weaker events, and smaller values introduced undesirable levels of background noise. The optimal display properties associated with the coherency filtered version of lines 88-17 and 14 are shown in Figure 6.10 and Figure 6.11. Traces in thesefiguresare plotted using the "variable area only" format (i.e., without wiggle traces). Note the improved clarity of the images compared with the "standard" displays in Figure 6.8 and Figure 6.9. Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS Figure 6.8: Upper 6 s of zero-offset stacking section of Line 88-17. Horizontal-to-vertical aspect ratio is 1:1 assuming a constant velocity of 6000m/s. The station elevations are shown above the seismic section. 92 Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS Figure 6.9: Upper 7 s of zero-offset stacking section of Line 88-14. Horizontal-to-vertical aspect ratio is 1:1 assuming a constant velocity of 6000 m/s. The station elevations are shown above the seismic section. 94 Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS 95 Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS Figure 6.10: Upper 6 s of final coherencyfilteredzero-offset stacking section of 88-17. Station elevation and stacking velocity information is shown above the seismic section. Both RMS and Dix interval (Dix 1955) velocities are tabulated. Horizontal-to-vertical aspect ratio is 1:1 assuming a constant velocity of 6000m/s. Dashed line box outlines area of data in Figure 6.7. 96 Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS 97 Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS 98 Figure 6.11: Upper 7 s of final coherency filtered zero-offset stacking section of 88-14. Station elevation and stacking velocity information is shown above the seismic section. Both RMS and Dix interval (Dix 1955) velocities are tabulated. Horizontal-to-vertical aspect ratio is 1:1 assuming a constant velocity of 6000 m/s. Approximate Depth (km) CD oâ€ž cp CO o â€¢'..< i- A \, 1 .. I:' - ,>, . V â€¢i (S) aim Avn-om \> 'I h S3 Chapter 6. POST-STACK 6.4 PROCESSING AND FINAL SEISMIC SECTIONS 100 Migration The structural stacks displayed in the last section are useful, but not ideal, for the purpose of interpretation. Although partial prestack migration was achieved with the DMO procedure, and better velocity values were determined, the stack still does not present dipping sub-surface structure in their "true" subsurface position. To achieve this, the stacked data are migrated. Migration seeks to restore dipping reflection events to their true subsurface locations, by shortening and steepening them as well as by moving them in the updip direction, and to collapse diffraction hyperbolae. After considerable testing of different migration modules in the INSIGHT package (reverse time finite-difference migration, F-X migration and phase shift migration), a phase shift migration algorithm ("PSMIG-FAST") was applied to achieve repositioning of dipping reflectors approximately to their true subsurface location. Generally, there are two significant factors associated with migration performance when applied to crustal datasets. One is the velocity field; migration algorithms are extremely sensitive to errors in the velocity structures. Migration with too high velocities (or overmigration) will result in "smiles" on the display. Another factor is ambient noise, which commonly dominates deep crustal seismic sections where S/N ratio are lower and velocities are higher. This causes a smearing of amplitudes on the migrated sections. Because a single spike on the structural section migrates to a semicircle on the migrated stack, in addition to smearing effects, any sparsely distributed bursts of noise in the input section will create smiles on the migrated section. When using the poststack module "PSMIGJFAST", full (i.e., maximum) aperture width on bothfines88-17 and 14 was specified. Based on Yilmaz (1987), a narrow aperture can introduce strong smearing as spurious, nearly horizontal events. This effect Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS 101 occurs for all types of migration algorithms if the maximum dip allowed for migration is severely restricted. The amount of aperture refers to the effective width, in traces, of the diffraction hyperbola along which summation occurs (the migration process entails summing energy along a hyperbolic trajectory and placing the output result at the hyperbola apex). In order to minimize such noise problems, the migration algorithm was applied to the coherency-filtered version of the zero-offset section without application of a plotting bias. The S/N enhancement provided by the coherencyfilteringimproves the results of the migration process. For all migrations, the DMO-corrected velocity functions were used because they should be more correct than those without the DMO correction (see Chapter 4). Figures 6.12 and 6.13 display the migrated coherency-filtered seismic sections for lines 88-17 and 14, respectively. The horizontal- to-vertical aspect ratio is 1:1 assuming a constant velocity medium of 6000 m/s. We note that migrated stacks shorten and steepen reflectors in structural stacks. Compare Figure 6.10 and 6.12. The dipping reflector in the structural stack at shot point number 922 and time about 4.0 s is moved updip in the migrated stack to shot point number 764 and time about 3.5 s. Compare Figure 6.11 and 6.13. The mid-crustal ramp in the structural stack at shot point number 408 and time 6.0 s is moved updip after migration, appearing at shot point number 250 and time 5.5 s. From Figure 6.12 and 6.13, we also noted that migrated sections appeared overmigrated based on the presence of "smile" artifacts. The reason could be that migration algorithms are extremely sensitive to the velocity structure function, and/or there are spurious gaps in crustal reflector continuity. In any later study, migration velocity should be tested on segments using different percentages of the "true" interval velocity (derived from RMS velocity profile), and like constant stack velocity analysis, the best percentage will be picked based on observation, which is commonly performed in the industry. NE l?19 127a 0 CO 1 Co O to O CO to a 10 Km 88-17 Figure 6.12: Migrated coherency filtered seismic section of 88-17. Horizontal - to - vertical aspect ratio is 1-1 assuming a constant velocity of 6000 m/s. CO s s O CO ta O o o to Chapter 6. POST-STACK PROCESSING AND FINAL SEISMIC SECTIONS 103 = 88-14 10 km Figure 6.13: Migrated coherency filtered seismic section of Line 88-14. Horizontal-to-vertical aspect ratio is 1:1 assuming a constant velocity of 6000 m/s. Chapter 7 I N T E R P R E T A T I O N A N D DISCUSSION 7.1 Introduction Varsek et al. (1993) have provided a geological interpretation of the seismic sections for Lines 88-14 and 88-17 as provided by the contractor (Figure 7.1). Key features include: (1) west-dipping structures truncate and therefore post-date east-dipping structures; (2) a significant west-dipping upper crustal ramp at the center of Line 88-17 underlies the Castle Pass (CPF) and Bralorne (BF) fault systems and may be related to contractional faults found to the northeast; (3) structural continuity between the southwest end of Lines 88-17 and 88-14 is inferred from a reflective band between about 1.5 to 3.5 s that occurs on both profiles and that geometrically resembles a hanging wall anticline formed on the Thomas Lake fault (TLF). An alternate interpretation of these profiles based on the same processed sections, but taking into account new geological information and viewpoints, is given by Monger and Journeay (1994) (Figure 7.2). They suggest more complicated structures on Line 88-17 and no impression of the Thomas Lake fault on Line 88-14. In this section, we integrate the reprocessed seismic reflection sections for Lines 88-14 and 88-17 and new near-surface geological data (Figure 7.3 and Figure 7.4) to provide a revised upper crustal interpretation and new insights on the history of deformation in the southern Coast Belt. 104 Chapter 7. INTERPRETATION AND DISCUSSION 105 Figure 7.1: Lines 88-14 and 88-17 initial interpretation superimposed on migrated and coherency-filtered section produced by the contractor processing (from Varsek et al. 1993). Key features of this interpretation are addressed in the text. Abbreviations: Bralorne fault system (BF), Castle Pass fault (CPF), Marshall Creek fault (MCF), Miller Creek fault (MIF), Owl Lake fault (OLF), Thomas Lake fault (TLF), Tyaughton Creek fault (TyF), Bridge River terrane (BR), Cadwallader terrane (Cd), Gambier Assemblage (G), Late Cretaceous granitic (lKg), Jurassic to Cretaceous granitic (JKg). Figure 7.2: Alternate interpretation of Lines 88-14 and 88-17 along the Carpenter Lake-Hurley River Road made by Monger and Journeay (1994). o Chapter 7. INTERPRETATION i Â£ AND DISCUSSION 107 m E Figure 7.3: Geological cross-section provided by Monger and Journeay (1994). It shows the fault systems of part of the southeastern Coast Belt. Locations of A-B, C-D and E-F are shown in Figure 7.4. 8 Legend 8 |^ A ^ J . | MPv: Mioccnc/Plioccnc volcanics |o ^ o | Rs: Roccnc sedimentary rocks j ^ ^ J j IKI'C: Powell Creek Poimnlion ".J KSq: Silvcrqiitck Conglomerate ^ KSO: Spcnccs Bridge Group |7< r [3."T-J K G : Gambler Group Â£ ' 7 \ ^ KT: Taylor Creek Group p'T^j S3 cKJ: Jackass Mountain Croup j . ' . . | JKRM-. Relay Mountain Group HTUn. " ' arT son ^-^ c B Formolion Jl-C: Last Creek Formation j"" o "j )KC: Cnyoosli Assemblage |^\>] TT: Tynngritnn Group l'-PO [ T C U : Cadwallader Group | O i l : Bridge River Complex j^^H PliL: lir;ili><ne-ITasl Liza Complex ||H Sh: Sltulnps Complex ^ to to MCC!: Cliisiit Creek Schist j'"* Â«| EGd: lioccnc pinions j^T^j 1 K | 1 : l - n ' Cretaceous Dentlor Suite c [ T ^ / J IKS: Laic Cretaceous Scuzzy Suite o |*~ ^ * | mKgd: iniil-Crclaccous pinions j / ' ,"| cKgil: liaily Cretaceous pinions t | ^ j mlJgil: Miil/1 .ale Jurassic plulous | ^ ^ | MSL: Slullictim Scliisl j* â€¢*"â€¢ j PI": Pcrinian-Trinssic plutons r ^ F j 1TCC. Cnclic Creek Complex s â€ž J T N : Nicola Group Figure 7.4: Recent geologic map of the eastern Coast Belt (Monger and Journeay 1994). Cross-sections of A-B C-D and E-F are shown in Figure 7.3. A-B is roughly the location of Line 88-13; C-D is the east part of Line 88-17- E-F locates northeast to Line 88-17. Â§ Chapter 7. INTERPRETATION AND DISCUSSION 7.2 Interpretation of reprocessed Line 109 88-14 The reprocessing of Line 88-14 has marginally improved some aspects of the image while other features are remained similar (Zhang and Clowes 1994). The surface geology of Line 88-14 is less complicated than Line 88-17. It is underlain mainly by underformed postorogenic granitic rocks (Journeay 1990). At its western end, Line 88-14 starts in Jura-Cretaceous granites and overlies the Gambier assemblage of early Cretaceous arc volcanics which are intruded by syntectonic plutonic rocks, then crosses Miller Creek fault to end on more granites (Figure 1.3). The shortcomings associated with crustal-scale migration (Chapter 6) on Line 88-14 preclude a reliable interpretation based solely on the migrated results. So our interpretation of Line 88-14 is displayed by concurrent analysis of both the structure stack (Figure 7.5) and migrated section (Figure 7.6). A , B , C, D, E and F are reflectors identified on structure stack of Line 88-14 (Figure 7.5). Events B , C, D, E , and F are more coherent and have greater lateral extend than previous processing. In particular, event F which is interpreted by Varsek et al. (1993) to be a mid-crustal east-dipping ramp, is shown to extend across 2/3 of the section. It is inferred to be a deep counterpart to west-directed shear zones found to the west (Ashlu Creek fault). West-dipping events A and B extend into shallow crustal levels (2 s) and B can be correlated with east-directed thrusts of Birkenhead/Grizzly Pass Fault. Line 88-14 overlies the Gambier assemblage of early Cretaceous arc volcanics. The Gambier volcanics appear to have depth extend of no more than 6 km ( about 2 s). The only mapped fault within the region, the Miller Creek east-directed thrust fault, shows very little expression on either the original or reprocessed section. 110 Chapter 7. INTERPRETATION AND DISCUSSION E CD Q CD E X o \â€” CL Q. < 88-14 Figure 7.5: Coherencyfiltered,structure stack for the upper 7 s of Line 88-14. A, B, C, D E and F are the main reflectors identified on the section. Chapter 7. INTERPRETATION = = AND DISCUSSION 111 88-14 10 km Figure 7.6: Coherency filtered, migrated structural section for the upper 7 s of Line 88-14. Approximate depth was calculated based on average velocity of 6000 m/s. Chapter 7. INTERPRETATION 7.3 AND DISCUSSION 112 Interpretation of reprocessed Line 88-17 Line 88-17 starts about 6 km further east of Line 88-14 across the Pemberton Valley at the position of the Owl Creek fault, which juxtaposes Triassic-Jurassic volcanic/sedimentary rock of the Cadwallader group over the Gambier assemblage. The line continues across a series of small terrane and granitic units, all part of the Eastern Coast belt (Figure 7.7): Cadwallader (TCdp), Triassic-Jurassic Harrison terrane island arc volcanics (JH), Carboniferous-Cretaceous Bridge River terrane including the Jurassic-Cretaceous Cayoosh assemblage turbidites (JKC) and Bridge River upper Paleozoic to Mesozoic subduction complex rocks (PBC), the late Cretaceous Hurley River quartz-diorite pluton, and the Jurassic-Cretaceous Methow terrane of turbidites derived from different tectonic settings (JKCm). The Une crosses a complex series of fault systems before terminating in the Bridge River terrane (CJB). Reprocessing of reflection Line 88-17 generated a significantly improved image of the upper 6 s (to about 20 km depth). AppUcation of refraction statics, cross-dip correction and truncated spread length contributed most to the improved image (Zhang and Clowes 1995). The reprocessing has better deUneated the subsurface geometry of some of the terranes and faults, confirming their Umited depth extent, although much of the 3-D complexity inferred from surface geology is not clearly expressed in the 2-D seismic images. In Figure 7.8, we present our interpretation based on analysis of the coherency filtered, migrated reprocessed section and other presentations of the data and the most recent geological information. Figure 7.8 shows that Line 88-17 crosses a series of small terrane units and these terranes are sUced by a complex series of NW-SE striking faults. Complex geology at the surface is repUcated by the complexity of the subsurface image: "sUced up" segments of the different terrane units are juxtaposed in an irregular pattern. Individual terrane Chapter 7. INTERPRETATION AND DISCUSSION 113 Figure 7.7: Geologic map with location of Line 88-17 (next two pages). (A) Geologic map with location of Line 88-17. Digital geological data was provided by M. Journeay. Black solid line is thefieldreflection seismic line; yellow dashed one is the slalom line on which the seismic data processing is based; red numbers arefieldstation numbers. (B) Legend of this geologic map. Chapter 7. INTERPRETATION AND DISCUSSION 114 Chapter 7. INTERPRETATION AND DISCUSSION 115 LEGEND STRATIFIED ROCKS Cretaceous Pleistocene and Recent Q PLUTONIC ROCKS Thick Drift j Kgd j i mKgd j KG Granodiorite (gd) and quartz diorite (qd) Mid-Cretaceous granodiorite (gd) in Coast Belt Gambier Group Early Jurassic - Late Jurassic (145 - 187 Ma) Taylor Creek Group Granodiorite (gd) in Western Coast Belt Lower Jurassic to Lower Cretaceous gB^jSgJB Cayoosh Assemblage- Early and Middle Jurassic Harrison Lake Formation OPHIOLITE COMPLEXES Permian Triassic Bralorne Complex Cadwallader Group Hurley Formation Pioneer Formation METAMORPHIC ASSEMBLAGES Carboniferous to Middle Jurassic Early Late and Late Cretaceous Â§Â£11 Bridge River Complex OBs | Radiolarian chert, siltstone CJBm i Light to dark grey phyllile Chism Creek Schist (B) Figure 7.8: Coherencyfiltered,migrated structural section for the upper 6 s of Line 88-17. A preliminary interpretation, based on the reprocessing seismic data and new geological information made available by M. Journeay (GSC, Vancouver), is superimposed on the section. Approximate depth was converted from time scale by using average velocities for Coast Belt (Clowes et al. 1995). BR - Bridge River complex, Ca - Cayoosh assemblage, Cd - Cadwallader group, Ga/qd Gambier group intruded by quartz - diorite granites, Ha - Harrison complex. 117 Chapter 7. INTERPRETATION AND DISCUSSION slices appear to be less than about 6 km (2 s) thick, but stacking and juxtaposition of the slices develop a terrane collage that may be up to 20 km thick. The following discussion of our interpretation takes place, in a general sense, from west to east. Mapped bounding faults (Owl Creek fault and Birkenhead/Grizzly Pass fault) of the Cadwallader terrane are imaged on our reprocessed seismic section. The Owl Creek fault on the west is shallow (< 3 km) and sub-horizontal. It is truncated by the Birken- head/Grizzly Pass fault at shot station 450 and 3 km depth (about 1 s). However, the Birkenhead/Grizzly Pass fault was not clearly imaged on the previous processing, and was not interpreted by Varsek et al. (1993). Mapping indicates deeper structural levels to the south. The reprocessed seismic section confirms with Varsek et al. (1993) that west-dipping structures truncate and therefore postdate east-dipping structures. The Noel Creek fault dips moderately SW and may extend into the middle crust at 16 km depth (5 s). At surface, it locates at shot station 900, which is consistent with geological surface data (Figure 7.7). The Noel Creek fault was not imaged on the previous processing and was not interpreted by Varsek et al. (1993) either. However we consider that the Noel Creek fault is a major structure which carries terranes identified west of it (in the hanging wall) over terranes identified east of it (in the footwall). East of the Noel Creek fault, Bridge River and Cadwallader terranes (perhaps with some Cayoosh assemblage) form the upper crust. Sub-vertical fault zones (e.g., Bralorne fault) are not imaged but our interpretation indicates that they only extend a few kilometers in depth, being truncated by a major detachment. This west dipping detachment may be correlated to the contractional faults found to be the northeast. Chapter 7. INTERPRETATION 7.4 AND DISCUSSION 118 Geological model and seismic modelling of Line 88-17 The main objective of seismic interpretation is to take the migrated section and produce from it an accurate geological section. Actually, in our study, an unmigrated seismic section represents the best we can get through data processing. Anything that follows thereafter falls into the realm of interpretation. So we established a geological model and carried out normal incidence forward modelling to see how the theoretical seismic section matches the observed stacked seismic data. Such a procedure can constrain our interpretation by indicating the appropriate seismic section geometries and the relative amplitudes of the different reflectors. The main objective of our modelling is to help guide an interpretation in structurally complex areas. Outrider is a seismic modelling package which can build models, do ray tracing, create seismic sections, and also perform time-to-depth conversion by inverting seismic sections to geological sections. It has two basic types of models, depth (geological) and time (seismic). Forward modelling starts from the geological profile and generates a time section. Inverse modelling starts from the seismic profile and generates a depth section. The geological model of Line 88-17, based on the interpretation in Figure 7.8, is shown in Figure 7.9. The model is built up by three steps: (1) create the events (geological horizon surfaces and faults) and define their interfaces; (2) tie all the events together; (3) build the model by eliminating all duplicate points, checking that all events are correctly ties and forming polygonal regions. The P-wave velocities and densities were kindly provided by Murray Journeay of the Geological Survey of Canada, Vancouver. In the forward modelling part, Outrider could produce four types of structural (post stack) synthetics: diffractions (examining the accuracy of time migration), zero-offset sections (simulating unmigrated sections), image ray sections (simulating migrated sections), and wavefront sections (modelling a series of Huygens point sources to produce Chapter 7. INTERPRETATION AND DISCUSSION 119 unmigrated seismic sections). With the tool of Outrider, we did the initial forward modelling to simulate an unmigrated section by creating zero-offset synthetic data. And, we think zero-offset modelling is the fastest way to model an unmigrated section. It was generated by tracing rays upwards from the reflecting boundaries. All rays for reflection from boundaries start at normal incidence (90 degrees) to the boundaries and refract at interfaces. Figure 7.10 shows seismic traveltime profile calculated by running normal incidence modelling on Figure 7.9. Figure 7.10 shows that synthetic seismic geometries basically match the observed data (compare with Figure 6.10 and Figure 6.12). It indicates our interpretation is geometrically appropriate since the model was based on the observed data interpretation. In this thesis, we did not focus on the accuracy of the seismic amplitudes, partly due to hardware limitations for a model of this size on the computer on which calculations were made. However, I would recommend that any future analysis should combine forward modelling from geological sections and inverse modelling from seismic sections. 7.5 Summary comments Interpretation of Lines 88-14, 88-17 and 88-18 (Perz 1993) portray the Eastern Coast Belt as a crustal-scale accretionary complex comprising an imbricate stack of tectonic flakes (Figures 7.6, 7.8 and 7.11). Bralorne fault defines the fabric of this accretionary complex in the Cayoosh and Bendor ranges. Structures within this accretionary complex record a history of west-directed thrusting ( Bralorne fault (steep), Owl Creek fault (shallow)) and east-directed thrusting (Castle Pass (steep), Noel Creek fault (shallow)) (Figure 7.8). Figure 7.11 shows the interpretation of Line 88-18 by Monger and Journeay (1994). Although recorded as a continuous profile, Line 88-18 can be considered as twofinesseparated into eastern and western halves (18E and 18W, respectively) at the Chapter 7. INTERPRETATION AND DISCUSSION 120 Â£<uwi_a)vi Figure 7.9: Geological model of Line 88-17 built by Outrider . The vertical exaggeration is 2:1. P-wave velocities and densities are listed in the figure. Units are m/s and g/crro, respectively. Shear-wave velocities, which are also used for calculating the reflection coefficients, are derived by assuming Poisson's ratio is 0.3. â€¢8 150 m s e c s 250 350 450 Shot Points 650 750 550 850 950 1050 1150 500 1000 150020002500 3000 3500 4000 4500 -I 5000 5500-i 1250 500 1000 1500 I-2000 2500 3000 3500 4000 4500 5000 5500 55600 50600 45000 40000 35600 30000 25000 20000 15000 10000 5TJ6T Offset (metres) Figure 7.10: Seismic traveltime profile from normal incidence modelling. The geological model is shown in Figure 7.9. 6 S3 S3 ta O to to & s o Chapter 7. INTERPRETATION AND DISCUSSION 122 Fraser River where up to 300 km of combined dextral displacement occured on the Fraser River fault system. The upper crust is characterized by a generally west dipping fabric, and the middle crust is characterized by horizontal reflections that are quasi-continuous beneath the Fraser River fault. Within the upper crust on profile 18E, the Pasayten fault is characterized by east dipping reflections. West directed thrust motion on the Pasayten fault is the dominant structural signature on reflection data. The western margin of the southeastern Coast Mountains accretionary complex is the Coast Belt Thrust System, which is interpreted to be a part of a crustal-scale ramp across which assembled terranes of the southeastern Coast Mountains were thrust southwest-ward over supracrustal arc sequences of the Insular Superterrane in Late Cretaceous time (Journeay 1990; Journeay and Friedman 1993). Figure 7.12 shows the interpretation of Line 8813 by Monger and Journeay (1994). The most significant observation from this line is that crosscutting relationships in the upper and middle crust indicate that east-dipping structures predate wets-dipping structures. Geometric evidence for east directed motion on the west-dipping midcrustal reflectors (1) the fault system at the west side of line is apparently folded above a ramp between 3.5 and 4.5 s and (2) easting-dipping surface faults are displaced eastward relative to the east-dipping lower crustal reflectors. The ramp-flat trajectory of mid-crustal reflections seen on Lines 88-17, 88-14 and 88-13 could be interpreted to be part of east-directed thrust system (Figures 7.6, 7.8 and 7.12; Figure 7.3). These structures cut across the basal decollement of the Coast Belt Thrust system, and apparently root westward beneath Wrangellia. From overview of reflection Lines 88-13, 88-14, 88-17 and 88-18, the southeastern Coast Mountains accretionary complex is similar in both scale and overall geometry to that formed in the upper plate between the present Cascadia subduction zone and Wrangellia terrane on Vancouver Island, in which slivers of oceanic plateau and ocean floor and arc complexes (Crescent and Pacific Rim terranes) have been imbricated and Fraser Fnull Figure 7.11: Interpretation of Line 88-18 by Monger and Journeay (1994). West dipping reflectors in the footwall of the east dipping Pasayten fault exhibit cutoffs consistent with east directed motion. Restoration of up to 300 km of dextral-slip displacement on the Fraser Fault implies that it must be part of a regional thrust system. The large east dipping lower crustal ramp, referred to here as the Insular ramp, is likely the suture between the Intermontane and Insular composite terranes. Chapter 7. INTERPRETATION AND DISCUSSION underplated beneath the Cretaceous continental margin represented by Wrangellia. 7.12: Interpretation of Line 88-13 along the Duffey Lake Road by Monger and Journeay (1994), showing integration of geological and geophysical data. See Figure 7.3 for geological cross-section. to References [1] Al-yahya, K.M. 1991. Application of the partial Karhunen-Loeve transform to suppress random noise in seismic sections. Geophysical Prospecting, 39: 77-93. [2] Barry, K.M. 1967. Delay time and its application to refraction profile interpretation, in Seismic Refraction Prospecting, edited by A.W. Musgrave, Society of Exploration Geophysicists, Tulsa, OK. [3] Cassidy, J.F. 1995. Review: Receiver function studies in the southern Canadian Cordillera. Canadian Journal of Earth Sciences, 32: 1514-1519. [4] Chun, J.H., and Jacewitz, C. 1981. Fundamentals of frequency-domain migration. Geophysics, 46: 717-732. [5] Claerbout, J. F. 1985. Fundamentals of Geophysical Data Processing, Blackwell Scientific Publications, 274 pp. [6] Clowes, R.M. (editor) 1993. Lithoprobe Phase IV Proposal - Studies of the Evolution of a Aontinment. Published by the Lithoprobe Secretariat, The University of British Columbia, Vancouver, B.C., 290 pp. [7] Clowes, R.M., Baird, D.J., and Dehler, S.A. 1997. Crustal structure of the Cascadia subduction zone, southwestern British Columbia, from potential field and seismic studies. Canadian Journal of Earth Sciences, 34: 317-335. [8] Clowes, R.M., Zelt, C.A., Amor, J.R., and Ellis, R.M. 1995. Lithospheric structure in the southern Canadian Cordillera from a network of seismic refraction lines. 126 References 127 Canadian Journal of Earth Sciences, 32: 1484-1513. [9] Clowes, R.M., Brandon, M., Green, A.G., Yorath, C.J, Sutherland-Brown, A., Kanasewich, E.R., and Spencer, C.S. 1987. Lithoprobe southern Vancouver Island: Cenozoic subduction complex imaged by deep seismic reflections. Canadian Journal of Earth Sciences, 24: 31-51. [10] Coney, P.J., Jones, D.L., and Monger, J.W.H. 1980. Cordillera suspect terranes. Nature, 288: 329-333. [11] Cook, F.A. (editor) 1995. The southern Canadian Cordillera transect of Lithoprobe. Canadian Journal of Earth Sciences, 32: 1483-1824. [12] Cook, F.A., Varsek, J.L., and Thurston, J.B. 1995. Tectonic significance of gravity and magnetic variations along the Lithoprobe Southern Canadian Cordillera Transect. Canadian Journal of Earth Sciences, 32: 1584-1610. [13] Cook, F.A., Varsek, J.L., Clowes, R.M., Kanasewich, E.R., Spencer, C.S., Parrish, R.R., Brown, R.L., Carr, S.D., Johnson, B.J., and Price, R.A. 1992. Lithoprobe crustal reflection cross section of the Southern Canadian Cordillera, 1, Foreland thrust and fold belt to Fraser River fault. Tectonics, 11: 12-35. [14] Cook, F.A., Green, A.G., Simony, P.S., Price, R.A., Parrish, R.R., Milkereit, B., Gordy, P.L., Brown, R.L., Coflin, K.C., and Patenaude, C. 1988. Lithoprobe Seismic reflection structure of the southeastern Canadian Cordillera: Initial results. Tectonics, 7: 157-180. [15] Coppens, F. 1985. First arrival picking on common-offset trace collections for automatic estimation of static corrections. Geophysical Prospecting, 33: 1212-1231. References 128 [16] Davis, G.A., Monger, J.W.H., and Burchfiel, B.C. 1978. Mesozoic construction of the Cordilleran "collage", central British Columbia to California. In Mesozoic paleography of the Western United States, Pacific Coast Paleography Symposium 2, edited by D.G. Howell and K.A. McDougall, Society of Economic Paleontologists and Mineralogists, Pacific Section, Los Angeles, California, pp. 1-32. [17] Dix, C H . 1955. Seismic velocities from surface measurements. Geophysics, 20: 6886. [18] Fookes, G. 1995. Karhunen-Loeve and spectral matrixfilteringin attenuation high amplitude peg-leg multiples. EAGE 57th Conference and Technical Exhibition, P004: 573-575. [19] Friedman, R.M., Mahoney, J.B., and Cui, Y. 1995. Magmatic evolution of the southern Coast Belt: constraints from Nd-Sr isotopic systematics and geochronology of the southern Coast Plutonic Complex. Canadian Journal of Earth Sciences, 32: 1681-1698. [20] Farrell, R.C, and Euwema, R.N. 1984. Refraction statics. Proceedings IEEE, 72: 1316-1330. [21] Freire, S. L.M., and Ulrych, T.J. 1988. Application of singular value decomposition to vertical seismic profiling. Geophysics, 53: 778-785. [22] Gardner, L.W. 1939. An areal plan of mapping subsurface structure by refraction shooting. Geophysics, 4: 259-274. [23] Gehrels, G.E., and Saleeby, J.B. 1985. Constraints and speculations on the displacement and accretionary history of the Alexander-Wrangellia-Peninsular superterrane. Geological Society of America, Abstracts with Programs, 17: 356. References 129 [24] Gelchinsky, B., and Shtivelman, V. 1983. Automatic picking of first arrivals and parameterization of traveltime curves. Geophysical Prospecting, 31: 915-928. [25] Gerbrands, J.J. 1981. On the relationship between SVD, KLT, and PCA. Pattern Recognition, 14: 375-381. [26] Gibson, B., and Larner, K. 1984. Predictive deconvolution and the zero-phase source. Geophysics, 49: 379-397. [27] Hagedoorn, G.J. 1959. The Plus-Minus method of interpreting seismic refraction sections. Geophysical Prospecting, 7: 158-182. [28] Hale, D. 1984. Dip-moveout by Fourier transform. Geophysics, 49: 741-757. [29] Hampson, D., and Russell, B. 1984. First-break interpretation using Generalized Linear Inversion. Journal of the Canadian Society of Exploration Geophysicists, 20: 40-54. [30] Hatherly, P.J. 1982. A computer method for determining seismic first arrival times. Geophysics, 47: 1431-1436. [31] Hyndman, R.D., and Lewis, T.J. 1995. Review: The thermal regime along the southern Canadian Cordillera Lithoprobe corridor. Canadian Journal of Earth Sciences, 32: 1611-1617. [32] Jones, A.G., and Gough, D.I. 1995. Electromagnetic images of crustal structures in southern and central Canadian Cordillera. Canadian Journal of Earth Sciences, 32: 1541-1563. [33] Jones, I.F., and Levy, S. 1987. Signal-to-noise ratio enhancement in multichannel References 130 seismic data via the Karhunen-Loeve transform. Geophysical Prospecting, 35: 1232. [34] Journeay, J.M. 1990. A progress report on the structural and tectonic framework of the southern Coast Belt, British Columbia. In Current Research, Part E, Geological Survey of Canada, Paper 90-1E: 183-195. [35] Journeay, J.M., and Friedman, R.M. 1993. The Coast Belt thrust system: evidence of late Cretaceous shortening in southwest British Columbia. Tectonics, 12: 756-775. [36] Kanasewich, E.R., Burianyk, M.J.A., Ellis, R.M., Clowes, R.M., White, D.T., Cote, T., Forsyth, D.A., Luetgert, J.H., and Spence, G.D. 1994. Crustal velocity structure of the Omineca Belt, southeastern Canadian Cordillera. Journal of Geophysical Research, 99: 2653-2670. [37] Lanczos, C. 1961. Linear differential operators, D. Van Nostrand Co. Ind., Ch. 3. [38] Lane, M. 1991. Coherency/semblancefiltering.LSPF Newsletter, 4(2): 19-26. [39] Lamer, K. L., Gibson, B.R., Chambers, R., and Wiggins, R.A. 1979. Simultaneous estimation of residual static and crossdip corrections. Geophysics, 44: 1175-1192. [40] Liner, C.L., and Bleistein, N. 1988. Comparative anatomy of common offset dip moveout. Expanded abstract of SEG 67th Conference and Technical Exhibition, S17.1: 1101-1105. [41] Marchisio, G., Pendrel, J.V., and Mattocks, B.W. 1988. Application of full and partial Karhunen-Loeve transformation in geophysical image enhancement, 58th Annual SEG meeting, Anaheim, California, Paper S 22.5: 1266-1269. References 131 [42] McGroeder, M.F. 1991. Reconciliation of two-sided thrusting, burial metamorphism and diachronous uplift in the Cascades of Washington and British Columbia. Geological Society of America Bulletin, v.103: 189-209. [43] Milkereit, B. 1989. Stacking charts: An effective way of handling survey, quality control and data processing information. Canadian Journal of Exploration Geophysicists, 25: 28-35. [44] Monger, J.W.H. 1985. Structural evolution of the southwestern Intermontane Belt, Ashcroft and Hope map-areas; British Columbia, In Current Research, part A, Geological Survey of Canada, Paper 85-1A: 349-358. [45] Monger, J.W.H., and Journeay, J.M. 1994. Guide to the Geology and Tectonic Evolution of the Southern Coast Mountains, May 1994. Geological Survey of Canada, Vancouver, B.C., open file 2490, 77 pp. [46] Monger, J.W.H., and Journeay, J.M. 1992. A field guide to accompany the Penrose Conference on the Tectonic Evolution of the Coast Mountains orogen, May 1992. Geological Survey of Canada, Vancouver, B.C., 97 pp. [47] Monger, J.W.H., and Price, R.A. 1979. Geodynamic evolution of the Canadian Cordillera - progress and problems. Canadian Journal of Earth Sciences, 16: 770-791. [48] Monger, J.W.H., Price, R.A., and Tempelman-Kluit, D.J. 1982. Tectonic accretion and the origin of the two major metamorphic and plutonic welts in the Canadian Cordillera. Geology, 10: 70-75. [49] Monger, J.W.H., van der Heyden, P., Journeay, J.M., Evenchick, C.A., and Mahoney, J.B. 1994. Jurassic-Cretaceous basins along the Canadian Coast Belt: Their bearing on pre-mid-Cretaceous sinistral displacements. Geology, 22: 175-178. References 132 [50] Neidell, N.S., and Taner, M.T. 1971. Semblance and other coherency measures for multichannel data. Geophysics, 34: 482-497. [51] O'Leary, D.M., Clowes, R.M., and Ellis, R.M. 1993. Crustal velocity structure in the southern Coast Belt, British Columbia. Canadian Journal of Earth Sciences, 30: 2389-2403. [52] Oppenheim, A.V. 1978. Application of Digital Signal Processing. Prentice-Hall, Inc. [53] Perz, M.J. 1993. Characterization of the Fraser fault, southwestern British Columbia, and surrounding geology through reprocessing of seismic reflection data. University of British Columbia, M.Sc. thesis, 170 pp. [54] [55] Robinson, E.A., and Treitel, S. 1980. Geophysical Signal Analysis, Prentice-Hall, New Jersey, 466 pp. [56] Roddick, J.A. 1983. Geophysical review and composition of the Coast Plutonic Complex, south of latitude 55 o N. In Circum-Pacific Plutonic Terranes. Edited by J.A. Roddick. Geological Society of America, Memoir 159, pp. 195-211. [57] Rutty, R.J., and Jackson, G.M. 1992. Wavefield decomposition using spectral matrix techniques. Exploration Geophysics, 23: 293-298. [58] Sacchi, M. 1996. Aperture compensated Radon and Fourier transform. University of British Columbia, Ph.D. thesis, 96 pp. [59] van der Heyden, P. 1992. A middle Jurassic to early Tertiary Andean-Sierran arc model for the Coast Belt of British Columbia. Tectonics, 11: 82-97. References 133 [60] Varsek, J.L., Cook, F.A., Clowes, R.M., Journeay, M., Monger, J.W.H., Parrish, R.R., Kanasewich, E.R., and Spencer, C.S. 1993. Lithoprobe crustal reflection structure of the southern Canadian Cordillera 2: Coast mountains transect. Tectonics, 12: 334-360. [61] Wheeler, J.O., and McFeely, P. (comp.) 1991. Tectonic Assemblage Map of the Canadian Cordillera and adjacent parts of the United States of America, Geological Survey of Canada, Map 1712A, scale 1:2 000 000. [62] Wiggins, R.A., Lamer, K., and Wisecup, R.D. 1976. Residual statics analysis as a general linear inverse problem. Geophysics, 41: 922-938. [63] Yilmaz, O. 1987. Seismic Data Processing, Society of Exploration Geophysicists, Tulsa, OK, 526 pp. [64] Zelt, B.C., Ellis, R.M., and Clowes, R.M. 1993. Crustal Velocity Structure in the Eastern Insular and Southernmost Coast Belts, Canadian Cordillera. Canadian Journal of Earth Sciences, 30: 1014-1027. [65] Zelt, B.C., Ellis, R.M., Clowes, R.M., Kanasewich, E.R., Asudeh, I., Luetgert, J.H., Hajnal, Z., Ikami, A., Spence, G.D., and Hyndman, R.D. 1992. Crust and upper mantle velocity structure of the Intermontane Belt, southern Canadian Cordillera. Canadian Journal of Earth Sciences, 29:1530-1548. [66] Zelt, C.A., and White, D.J. 1996. Crustal structure and tectonics of the southeastern Canadian Cordillera. Journal of Geophysical Research, 100: 24,255-24,274. [67] Zhang, W., and Clowes, R.M. 1995. Reprocessing southern Cordillera reflection data: Lines 88-14 and 88-17. LSPF Newsletter, 8(2): 12-16. References 134 [68] Zhang, W., and Clowes, R.M. 1994. Structure of the Eastern and Western Coast Belts, southwestern Canadian Cordillera, from reprocessing of Lithoprobe reflection data. In Canadian Geophysical Union Annual Meeting, Banff, Alberta, Programs and Abstracts, Abstract 12. Appendix A Dip-moveout by Fourier transform DMO, as the name implies, is a dip-dependent process. The two-dimensional (2-D) Fourier transform (over to and y) provides a particularly useful domain in which to apply a dip-dependent process, because all events having a particular slope in the (t , y)-domain 0 transforms to a single radial line in the (w ,fc)-domaingiven by 0 k u> Ato Ay 0 2sin9 v (A.l) where the 2-D Fourier transform of p (to,y,h) is defined by 0 P (^o, M ) = / dt e^ t0 0 0 J dye- P {t , y, h). iky 0 0 (A.2) This Fourier transform may be expressed in terms of NMO-corrected data by using the following change of integration variable, which was suggested by equation 4.4: ^ = ^ + 4^W^ ] 1 / 2 = + { ^y h 2 ? / 2 ( A 3 ) Defining and using equation (4.9) to replacePo(y^ + {Ato/Ay) h 2 2 2 y y, h) = p (t , n n y, h), the Fourier transform becomes Po(co ,M) = / dtnA-'e'^j 0 dye- yp (t ,y,h). ik n n (A.5) Using this integral transformation, which depends (through A) on the unknown zerooffset slope Ato/Ay, may seem even less practical than using the simpler equation (4.9). 135 136 Appendix A. Dip-moveout hy Fourier transform But equation (A.l) implies that the unknown slope may be interpreted as the ratio k/u> , 0 which suggests that A be defined as A = A(t ,w ,h,k) n 0 =(l+ ^y' 1 (A.6) 2 Finally, inverse Fourier transform Po(u>o,k,k) to obtain the desired "zero-offset" sections Po(t ,y,h), 0 i.e., Po{t , 0 yi )=^f h dwoe-^ J dke yp (u> , *, h). ik 0 0 (A.7) Equation A.5 resembles 2-D Fourier transform and is implemented digitally using fast Fourier transform (EFT) algorithm. The transform over y is, in fact, a Fourier transform and may be implemented digitally using any fast Fourier transform (FFT) algorithm. The transform over t , though not a Fourier transform, may be performed by numericn ally integrating over all t for each LO and k. The DMO by Fourier transform algorithm n 0 of equations A.5 to A.7 is, for constant velocity, accurate for all offsets and all dips. The differences between this algorithm andfinite-differenceDMO algorithm are analogous to the differences between frequency-wavenumber andfinite-differencealgorithm for migration.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Structure of the Eastern Coast Belt, southwestern Canadian...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Structure of the Eastern Coast Belt, southwestern Canadian Cordillera, from reproducing lithoprobe seismic… Zhang, Weimin 1997
pdf
Page Metadata
Item Metadata
Title | Structure of the Eastern Coast Belt, southwestern Canadian Cordillera, from reproducing lithoprobe seismic reflection lines 88-17 and 88-14 |
Creator |
Zhang, Weimin |
Date Issued | 1997 |
Description | As part of the Lithoprobe Southern Cordillera research program, seismic reflection data were acquired in 1988 along a 100 km profile across the Eastern and Central Coast Belts (lines 88-14 and 88-17) in a region of complex surface geology. A very good interpretation of the crustal section based on previous processing by industry has been made. However, the data were not processed specifically to enhance the features in the upper 5 s and provide a clear tie to surface geology, some of which is newly studied since the earlier interpretation. The present study involves reprocessing of these data to provide this enhancement. Severe crookedness of the line necessitated the implementation of unconventional processing techniques, including creation of several mini-profiles which cut through the midpoint scatter and application of a first-order correction for effects of reflector crossdip ( i.e., the crossline component of reflector dip). Refraction-based statics were computed using a two-step procedure consisting of initial identification and correction of systematic errors in the picked first break traveltimes, and subsequent application of a 2-D statics algorithm. Dip moveout (DMO) was applied to constant-offset gathers to effect prestack partial migration and reduce the dependence of normal moveout velocity on reflector dip. To reduce possible smearing in the upper few seconds, stacks were repeated with truncated spread lengths. The Karhunen-Loeve transform was applied efficiently and economically on the stacked seismic data to enhance coherent signals in all selected dips. The reprocessing has better delineated the subsurface geometry of the terranes and faults and expressed some of the structures which are not visible on previous processing sections, even though the 2-D seismic images cannot completely express the 3-D complexity inferred from the surface geology. The reprocessed sections are considered in light of the new geological mapping to provide a refined interpretation of subsurface structure and its implication for local tectonics. A new geological model is developed based on this refined interpretation and is being further tested by seismic forward modeling. |
Extent | 17579921 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-04-29 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0053000 |
URI | http://hdl.handle.net/2429/7693 |
Degree |
Master of Science - MSc |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1998-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1998-0077.pdf [ 16.77MB ]
- Metadata
- JSON: 831-1.0053000.json
- JSON-LD: 831-1.0053000-ld.json
- RDF/XML (Pretty): 831-1.0053000-rdf.xml
- RDF/JSON: 831-1.0053000-rdf.json
- Turtle: 831-1.0053000-turtle.txt
- N-Triples: 831-1.0053000-rdf-ntriples.txt
- Original Record: 831-1.0053000-source.json
- Full Text
- 831-1.0053000-fulltext.txt
- Citation
- 831-1.0053000.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0053000/manifest