Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Enhancement and modelling of signals from natural piezoelectric targets Butler, Karl Edward 1991

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1991_A6_7 B87.pdf [ 7.91MB ]
Metadata
JSON: 831-1.0052923.json
JSON-LD: 831-1.0052923-ld.json
RDF/XML (Pretty): 831-1.0052923-rdf.xml
RDF/JSON: 831-1.0052923-rdf.json
Turtle: 831-1.0052923-turtle.txt
N-Triples: 831-1.0052923-rdf-ntriples.txt
Original Record: 831-1.0052923-source.json
Full Text
831-1.0052923-fulltext.txt
Citation
831-1.0052923.ris

Full Text

ENHANCEMENT AND MODELLING OF SIGNALS FROM NATURAL PIEZOELECTRIC TARGETS by KARL EDWARD BUTLER B.Sc, Queen's University, 1987 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES (Department of Geophysics and Astronomy) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA October 1991 © Karl Edward Butler, 1991 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. Department of Aggfl/VSltf vrud fijTJJD^Oft^ The University of. British Columbia Vancouver, Canada Date CtfoftetV W , [CJFW DE-6 (2/88) ii A B S T R A C T Electromagnetic signals from piezoelectric targets in the earth can be generated using seismic sources and measured with electric or magnetic field receivers. The signals are typically small compared to the ambient electromagnetic noise and are difficult to identify in unprocessed records. Three data processing algorithms involving stacking, low-pass filtering, and sinusoid subtraction have been developed to enhance the signal-to-noise ratio of piezoelectric data acquired during field experiments. In addition, an analytic modelling technique has been developed to investigate the relationship between seismic and piezoelectric signals. The stacking technique, designed for use with repetitive seismic sources, employs a robust triggering algorithm that enables it to be used effectively even when the trigger signal is poor. High frequency noise is attenuated using a zero phase frequency domain low-pass filter with variable cut-off frequency and slope. The sinusoid subtraction technique is used to remove powerline noise which occurs at frequencies of 60 Hz and its harmonics. The amplitude and phase of each harmonic are estimated by calculating the Fourier series coefficients for that frequency. A sinusoid having the estimated amplitude and phase is then subtracted from the data to remove the harmonic. Remarkable improvements in the signal-to-noise ratio have been achieved by sinusoid subtraction as powerline noise levels typically exceed piezoelectric signal amplitudes by factors ranging from five to a few hundred times. For purposes of analytic modelling, a piezoelectric target is represented by a number of spheres each of which become independently polarized during the passage of a spherical elastic wave. The electric potential at a point in a uniform conductive medium surrounding the target is estimated by summing the potentials due to each of the polarized spheres. It is shown that the form of the electric potential time series generated by a uniformly polarized target is approximately proportional to the first derivative of the elastic wave particle velocity. In contrast, the potentials generated by non-uniformly polarized targets, representing quartz veins, are approximately proportional to the particle velocity itself. A distinct pulse of electric potential is also generated when the elastic wave reaches the outer limits of a vein-like target; such signals could yield important information on target dimensions if they could be identified in real piezoelectric data. iii TABLE OF CONTENTS Abstract ii Table of Contents iii List of Tables iv List of Figures v Acknowledgements vfi Chapter 1: Introduction 1 1.1 Piezoelectricity 2 1.2 Piezoelectricity in Rocks 7 1.3 The Piezoelectric Exploration Method 11 1.4 History of the Piezoelectric Method 14 (a) Development of the PEM in the USSR 14 (b) Development of the PEM outside of the USSR 18 1.6 Other Seismoelectric Effects 23 Chapter 2: Data Analysis 27 2.1 Introduction 27 2.2 Characteristics of Signals and Noise in Data Acquired by the PEM 27 (a) Characteristics of Piezoelectric Signals 27 (b) Sources of Noise 29 (c) Example: Piezoelectric Data from the Paymaster Quartz Vein Site..34 2.3 Processing Techniques for Improving the Signal to Noise Ratio 40 (a) Stacking using a Robust Triggering Algorithm 40 (b) Sinusoid Subtraction 53 (c) Low-pass Filtering in the Frequency Domain 63 Chapter 3: Modelling 68 3.1 Introduction 68 3.2 Theoretical Basis 69 (a) Electric Potential due to a Polarized Sphere 69 (b) A More General Model 72 (c) A Model for the Spherical Elastic Wave 74 (d) A Model for the Piezoelectric Fabric 81 3.3 Simulations 83 Chapter 4: Conclusions 98 References 102 Appendix A 106 Appendix B 135 L IST OF TABLES 3-1 Constants used to specify the spherical elastic wave LIST OF FIGURES V 1-1 Left and Right-handed Quartz Crystals 6 1- 2 Schematic Illustration of the Piezoelectric Exploration Method 13 2- 1 Plan View of the Paymaster Site 35 2-2 PAY60: A Typical Piezoelectric Record from the Paymaster Site 36 2-3 Demonstration of Signal Moveout at Paymaster 39 2-4 B37 Pre-stack Data 41 2-5 Plan View of the Layout used to Acquire data set B37 at Fairview 49 2-6 B37 Trigger Signals 50 2-7 B37 Mean Stacks 51 2-8 B37 Stacks after Sinusoid Subtraction using first 50 ms to Estimate Powerline noise.52 2-9 PAY59: A Raw (unprocessed) data set from Paymaster 58 2-10 PAY59: After Sinusoid Subtraction 59 2-11 Successive sinusoid subtractions on Dipole 1 60 2-12 Powerline Harmonics subtracted from Dipole 1 61 2-13 B37 Stacks after Standard Sinusoid Subtraction 62 2-14 Low-pass Filter Transfer Functions 66 2- 15 PAY60 after Sinusoid Subtraction and Low-pass filtering 67 3- 1 Spherical Model used to Represent a Piezoelectric Body 70 3-2 Model used to estimate Potential Due to a Piezoelectric Body of Arbitrary Shape 70 3-3 Model used to estimate the particle displacement generated by an explosive source...80 3-4 Symmetry of the Piezoelectric Fabric °°m 80 3-5 Geometry used for Piezoelectric Modelling 3-6 Pressure, Particle Displacement, and Particle Velocity 3-7 Comparison of Particle Velocity to Electric Potential 3-8 Electric Potential along Line Xj=18m in the x,x2 Plane 3-9 Decay of Normalized Electric Potential with Distance B-l Spherical and Cartesian Coordinate Systems with a Common Origin vii ACKNOWLEDGEMENTS I have had the good fortune of sharing the excitement and frustrations of seismoelectric research with a dedicated and diverse group of people. Special thanks are due to my supervisor Dr. R. Don Russell for his guidance and open-minded support. I express thanks to Mike Maxwell, Anton Kepic, Garfield Mellema, and Barry Narod for the genuine interest they took in my work, for their suggestions and of course for their expertise in the collection of field data. Dr. Tad Ulrych provided advice on data processing topics and kindly reviewed my thesis on short notice. I thank Ms. Rebecca Mills for her selfless support, and encouragement. This research was supported by Natural Sciences and Engineering Research Council (NSERC) Grants CRD-0094237 and A-0720 to R.D. Russell. The author was supported by an NSERC Postgraduate Scholarship. 1 1 I N T R O D U C T I O N The piezoelectric exploration method (PEM) involves the use of a source to generate seismic waves in the subsurface and receivers to measure the electromagnetic signals that are induced by the passage of the seismic wave through a piezoelectric body. The method can be used to search directly for deposits of piezoelectric minerals such as quartz, tourmaline, and sphalerite and indirectly for other valuable minerals that may be associated with these deposits. Exploration for quartz veins is expected to be the primary application as quartz veins are a common host for gold but are difficult to detect with conventional geophysical techniques. The principal objective of this thesis was to develop methods of data analysis and modelling that would aid in the understanding and development of the piezoelectric exploration method. The work was motivated by and complements field, laboratory, and theoretical studies of the method that researchers at the University of British Columbia (UBC) have been conducting since the late 1970's. The four chapters of the thesis are titled Introduction, Data Analysis, Modelling, and Conclusions. Piezoelectricity in crystals and rocks, the PEM, and other seismoelectric phenomena are reviewed in the Introduction. In Chapter 2 (Data Analysis) the characteristics of signals and noise recorded by the piezoelectric method are presented. This is followed by descriptions of processing techniques that have been developed to enhance the signal-to-noise ratio. Real piezoelectric field data is used to illustrate signal characteristics and the performance of the processing methods. In Chapter 3 (Modelling) the theoretical basis of an analytic model is presented followed by some modelling results. Synthetic time series of electric potential are generated for cases in which a spherically spreading, compressional seismic wave is incident upon a target approximated by an assemblage of piezoelectric spheres. The synthetic records 2 generated when a spherically spreading seismic wave passes through a planar piezoelectric target are of particular interest. All field data appearing in the figures of this thesis have been acquired during experiments conducted by the UBC seismoelectric research group (also known as the Geophysical Instrumentation Group). The group is directed by Dr. R.D. Russell. Its other members include Dr. M. Maxwell, A.W. Kepic, G.R. Mellema, Dr. B.B. Narod, and the author. 1.1 Piezoelectricity A piezoelectric material becomes electrically polarized when it is subjected to mechanical stress. The polarization vector P; is related to the stress O j k by the tensor equation Pi = dijk Ojk (i-1,2,3; j=l,2,3; k=l,2,3) (1.1) where the d i j k are constants called the piezoelectric moduli of the material. Note that each component of the polarization vector is linearly related to all components of the stress tensor. In cases where the material has a spontaneous polarization, P; must be interpreted as the change in polarization caused by the stress. Equation 1.1 expresses the direct piezoelectric effect. All piezoelectric materials also exhibit the converse piezoelectric effect which states that the application of an electric field E ; to the material causes it to change its shape according to the tensor equation e^diikEi (1.2) 3 where the components of the strain tensor describe the change in shape. The coefficients djjk relating polarization to stress in the direct effect are the same as those relating strain to the electric field in the converse effect. In general there are 27 components in a third rank tensor such as cL,^ . However, since the strain tensor is symmetrical (ie. = ekj) it follows from equation (1.2) that cLjk is also symmetrical in j and k. This reduces the number of independent d^ to a maximum of 18 and allows a more concise matrix notation to be used in place of the full tensor notation. In matrix notation the piezoelectric moduli are identified by new labels dy (i=l,2,3 ; 1=1,2,. 6). The corresponding matrix notation label for modulus dyk is given by the following scheme (Nye, 1957, p. 113): forj=k, d i j k = d i j , f4 forjk=23or32 for j * k, dijk = -T- where 1 = \ 5 for jk=13 or 31 16 forjk=12or21 Note that the first subscript is the same in both notations. In the absence of body torques (torques proportional to volume) the stress tensor Oj k is also symmetric (Nye, 1957) and its components may be assigned labels that are consistent with those used in the piezoelectric matrix by making the following change in notation: o"ll °"l2 ° 1 3 0"l ° 6 0"5 ° 2 1 <*22 ° 2 3 => 0"6 a 2 <*4 _ a 3 1 ° 3 2 a 3 3 _ -0"5 0-4 ° 3 In accordance with the convention used by Nye (1957) the longitudinal stress components (or the normal components of stress) <5y, o2, and o"3 are positive for tensile stresses and negative for compressive stresses. 4 In matrix notation, the direct piezoelectric effect given by equation (1.1) may be rewritten P-aVo - j (i=l,2,3;j=l,2,...6). (1.3) The principles of piezoelectricity in crystals have been studied extensively and are well documented (Cady, 1946; Nye, 1957). The number of independent moduli d^  in the piezoelectric matrix for a given crystal is dependent upon the crystal's symmetry. Centro-symmetric crystals cannot be piezoelectric (all 18 moduli are zero) since a centre of symmetry will not allow charges of opposite polarities to appear on opposite faces. Crystals with other symmetry properties have a variable number of non-zero piezoelectric moduli some of which are equal in magnitude. A summary of the forms of the piezoelectric matrices for all crystal classes is provided by Nye (1957, p. 124). According to Parkhomenko (1971, p. 36), the most extensive review of piezoelectric minerals is that of Bond (1943) who examined the symmetry of 830 minerals and determined that 70 of them should be piezoelectric. Piezoelectric moduli have been determined quantitatively for at least 22 minerals (Parkhomenko, 1971). Quartz is the most common piezoelectric mineral in rocks; nepheline, sphalerite, and tourmaline - although significantly less common than quartz - are among the more abundant of the remaining piezoelectric minerals (Sobolev andDemin, 1980, p. 18) Quartz is a primary rock-forming mineral second only to the feldspars with respect to abundance in the earth's crust (Dietrich and Skinner, 1979, p. 32). It is an enantiomorphic mineral meaning that it may crystallize in either right-handed or left-handed forms that are mirror images of each other. It is also polymorphic having four forms that are stable at different temperature ranges at a given pressure (Parkhomenko, 1971, p. 38). The low 5 temperature form - stable at temperatures below 573 degrees Celsius at atmospheric pressure - is called a-quartz. A member of symmetry group '32', oc-quartz has one three-fold axis (the c-axis) and three two-fold axes (the a-axes) of symmetry. The a-axes lie 120 degrees apart in a plane perpendicular to the c-axis as shown in Figure l-l(a). The piezoelectric matrix (cLj) for a-quartz has five non-zero piezoelectric moduli two of which are independent. The matrix is presented below along with the corresponding components of polarization and stress. <*1 <*2 <*3 <*4 °"5 0"6 Pi f d n 0 d 1 4 0 0 " p 2 0 0 0 0 -d 1 4 -2d„ p 3 . 0 0 0 0 0 0 This matrix is defined with reference to an orthogonal coordinate system xyz in which one of the a-axes coincides with * and the c-axis with z as shown in Figure 1-1 (b). The fact that all moduli on the bottom row of the matrix are zero illustrates that the c-axis can never be polarized regardless of the nature of the stress. The polarization vector must therefore lie in the same plane as the a-axes. The measured values of the piezoelectric moduli for right-handed a-quartz (Nye, 1957, p. 126) are d n = -2.3 pC/N and d 1 4 = -0.67 pC/N where the abbreviation pC/N means picocoulombs/newton. For a left-handed crystal, these values are multiplied by -1 (ie. d n = 2.3 pC/N and d 1 4 = 0.67 pC/N). At the temperatures that are encountered near the earth's surface where piezoelectric exploration may be utilized, a-quartz is the only stable quartz polymorph. The terms quartz and a-quartz will therefore be used interchangeably throughout the rest of this dissertation. 6 Left-handed quartz Right-handed quartz x = t a x-1-a z z Figure l-l (adapted from Bishop, 1978): (a) Views of left and right-handed quartz crystals showing the crystallographic axes. The c-axis (three-fold axis of rotation) is directed out of the page. The three a-axes (two-fold axes of rotation) lie in a plane parallel to the page, (b) Views of left and right-handed quartz showing the orientation of an orthogonal coordinate system according to the IRE (1949) conventions. 7 1.2 Piezoelectricity in Rocks The piezoelectric exploration method is intended to identify piezoelectric bodies in the subsurface. Such bodies are aggregates of piezoelectric crystals often mixed with non-piezoelectric minerals. Although each piezoelectric crystal in the aggregate would become polarized by the stress associated with the passage of a seismic wave, there is no guarantee that the body as a whole would exhibit a macroscopic polarization. If, for example, the piezoelectric axes of all crystals in a large body of rock were randomly oriented, the polarizations of individual crystals would be expected to sum to zero leaving the body with no net polarization. On the other hand, if the piezoelectric axes of all crystals were aligned in such a way that the polarizations of individual crystals tended to add constructively, the body would exhibit a net polarization and would make an ideal target for the PEM. If the piezoelectric axes of individual crystals (a-axes in the case of quartz) in an aggregate have a polar preferred orientation, the aggregate is said to have a piezoelectric fabric. Unlike piezoelectric effects in crystals, piezoelectric fabrics in natural aggregates (such as rocks) have received relatively little attention and their existence is supported by only a few experimental studies. A concise review of previous investigations of piezoelectric fabrics can be found in a recent paper by Ghomshei and Templeton (1989). A more comprehensive review of studies carried out prior to 1978 can be found in the Ph.D. thesis of J.R. Bishop (Bishop, 1978). A summary of research into piezoelectric fabrics -somewhat less comprehensive than that provided in either of the above references - is given below. According to Bishop (1978,1981) piezoelectric aggregates of fine quartz crystals were produced by Meissner in 1929 but the first detailed theoretical study of piezoelectric 8 fabrics in aggregates was conducted by the Soviet scientist Shubnikov (1946). Shubnikov extended the concepts of crystallographic symmetry to determine the piezoelectric moduli for aggregates with various symmetries (Bishop, 1978, p. 55). The first theoretical and laboratory investigations of piezoelectric effects in rocks were carried out by Parkhomenko and Volarovich beginning in 1953. Their studies were summarized in a book by Parkhomenko (1971). Laboratory measurements of piezoelectric effects in rocks are typically made by clamping a cubic or rectangular rock sample in a vise-like apparatus, applying a uniaxial stress to the sample, and measuring the component of polarization that appears across opposite faces of the sample. The polarization that is measured between the two clamped sample faces (ie. the component of the polarization vector parallel to the direction of uniaxial stress) is indicative of the strength of the longitudinal piezoelectric effect in the sample. Rocks bearing piezoelectric minerals may demonstrate a piezoelectric effect due to either a statistical effect or a piezoelectric fabric. The statistical effect refers to a statistical deviation from a random distribution for like ends of all piezoelectric crystal axes in the sample. Assuming a constant grain size, the effect is proportional to N 1 / 2 where N is the number of grains in the sample (Tuck et al. 1977; Bishop 1981). A statistical piezoelectric effect can also be caused by the presence of a single large uncompensated crystal in a sample otherwise composed of much finer grained crystals (Tuck et al., 1977). When relatively small piezoelectric effects are measured in rock samples it can be difficult to determine whether the effect is due to a piezoelectric fabric or merely statistical. This is particularly true for coarse-grained samples which contain relatively few crystals. Most of the rocks tested in the laboratory by Parkhomenko (1971) were quartz-bearing. The piezoelectric effects measured in about 200 samples of vein quartz and 9 pegmatitic quartz from several locations in the Soviet Union varied from 0.06% to 20% of the effect in monocrystalline quartz. Coarse grained samples tended to yield the larger effects while fine-grained samples exhibited the smallest effects. Parkhomenko found that vein and pegmatitic quartz tended to exhibit significantly larger piezoelectric effects than other quartz-bearing rocks including quartzites, granites, gneisses, and sandstones. In sandstone samples the effect was generally at least two orders of magnitude weaker than it was for typical samples of vein quartz. Parkhomenko attributed the large piezoelectric effects in vein quartz samples to the existence of piezoelectric fabrics. Tuck et al. (1977) measured the piezoelectric effects in samples of granite and in a highly aligned quartzite and found that their measurements were no larger than the expected statistical effect. The absence of a fabric-related piezoelectric effect in the quartzite sample was particularly informative since the c-axes and a-axes of the quartz grains in the quartzite were known to be highly aligned. On the basis of their results Tuck et al. concluded that natural piezoelectric fabrics probably do not exist. They surmised that the strong piezoelectric effects reported by Parkhomenko (1971) for vein and pegmatitic quartz were 'statistical effects arising from the dominance of one or a few quartz crystals or from an assembly of crystals produced by fragmentation of a larger one' (Tuck et al., 1977). Bishop (1978,1981) measured piezoelectric effects in several types of quartz-rich rocks and found that some of them - one quartz mylonite and two quartzites - exhibited a true piezoelectric effect due to fabric. No samples of pegmatitic or vein quartz were tested. The following three criteria were used to distinguish between statistical and fabric-related effects: (i) The measured effects were compared to the expected statistical effect for each sample. 10 (ii) Multiple specimens with the same orientation were cut from a single rock sample. Corresponding piezoelectric measurements on the different specimens were expected to be similar for rocks with a piezoelectric fabric. (iii) The optically observed c-axis orientation distribution was compared with that predicted by inversion of the piezoelectric data. The maximum effects measured by Bishop were approximately 1% of the effect due to a single quartz crystal. Ghomshei and Templeton (1989) measured longitudinal piezoelectric effects in many cubic samples from separated points along a quartz vein from the Fairview mining district in British Columbia. The samples were cut parallel to the vein wall and their orientations were preserved so that corresponding measurements in the different samples could be compared. The results suggested the existence of a piezoelectric fabric that was weakly coherent over distances of several tens of metres. The measured longitudinal piezoelectric coefficients averaged 2-3% of the longitudinal modulus d n for a single quartz crystal while the expected statistical effect was less than 0.5% of d n . Independent supporting evidence for a quasi-homogeneous piezoelectric fabric was provided by neutron diffraction analyses which were used to determine the orientation distribution (but not the polarities) of the a-axes in the vein samples. These analyses demonstrated the presence of a quasi-homogeneous a-axes fabric. Another form of independent, although indirect, evidence for the existence of a quasi-homogeneous piezoelectric fabric is provided by field experiments that UBC researchers have carried out at the vein. Electrical signals that appeared to be piezoelectric were observed in the field when a seismic stress was applied to the vein (Butler et al., 1991). In the absence of any preferential alignment of the piezoelectric axes in the vein, one would not expect to observe piezoelectric signals at the field scale. 11 Studies of a single spherical sample from the same quartz vein (Ghomshei et al., 1988; Ghomshei and Templeton, 1989) revealed the possible form of the piezoelectric fabric. Ghomshei et al. (1988) measured the longitudinal piezoelectric effect for 150 different orientations of the spherical sample and showed that the a-axes orientation distribution suggested by the piezoelectric data was in excellent agreement with that determined by neutron diffraction analyses of the same sample. The fabric suggested by the piezoelectric data could be approximated, to the first order, by the ideal quartz aggregate piezoelectric fabric 'infinity-m'. This fabric requires that there be equal proportions of right and left-handed quartz, that the c-axes be randomly oriented in a plane, and that one a-axis from each crystal be oriented polar parallel1 to a common direction. More information on the form of the 'infinity-m' fabric are presented in Chapter 4 where it is used in modelling the piezoelectric effect of a quartz vein. 1.3 The Piezoelectric Exploration Method (PEM) As shown in Figure 1-2 the PEM involves the use of a source to generate seismic waves in the subsurface and receivers to measure the electromagnetic signals that are induced by the passage of the seismic wave through a piezoelectric body. The sources and receivers may be placed on the surface, in boreholes, or inside mines. Since electromagnetic signals travel a few orders of magnitude faster than seismic waves, the time interval between the moment of source initiation and the reception of a piezoelectric response is essentially equal to the seismic traveltime from the shotpoint to the piezoelectric target. The distance from the shotpoint to the target can be calculated from the traveltime provided that a seismic velocity model is available for the survey area. Explosives appear to provide the most effective seismic source although weight drops and stacks of multiple jackhammer blows have been 1 Two piezoelectric axes are said to be polar parallel if they become electrically polarized in the same direction when subjected to the same stress. 12 used in experiments with varying degrees of success. In theory, it should be possible to measure the piezoelectric signals with both magnetic and electric field receivers. However, at the relatively low frequencies (typically lower than one kilohertz) used in the PEM it is expected (Russell and Barker, 1991) that the magnetic component of the signal will be much more difficult to detect than the electric component. The UBC research group therefore tends to use mainly electric field receivers such as grounded dipoles and capacitively coupled long wire antennas. Soviet scientists report that they have used both electric and magnetic field receivers to measure piezoelectric responses (Sobolev and Demin, 1980). 13 Figure 1-2 (from Sobolev et al., 1984): Schematic illustration of the piezoelectric exploration method. A seismic wave produced by an explosion or other source propagates through the earth and stresses a piezoelectric target (active zone). The target converts some of the seismic energy to an electromagnetic field which is detected by electric or magnetic field sensors linked to the recording apparatus. The shot moment circuit provides a signal indicating the time of detonation. 14 1.4 History of the Piezoelectric Method 1.4(a) Development of the PEM in the USSR Encouraged by their ability to observe the piezoelectric effect in rocks at the laboratory scale, Soviet scientists including M.P. Volarovich, E.I. Parkhomenko and G. A. Sobolev performed the first field tests of the piezoelectric method in the late 1950's. These tests were carried out at various surface sites in the USSR in areas of massive quartz-rich rocks (Volarovich et al., 1959 in Volarovich et al., 1962) and in the vicinity of pegmatite and quartz veins (Volarovich et al., 1962). Descriptions of the experiments were also provided by Sobolev and Demin (1980, p. 21-29). The seismic source consisted of a 100 kg weight dropped from a height of four metres to the bedrock surface. An eleven channel recorder incorporating amplifiers, bandpass filters, and galvanometers was used to record the responses of up to five seismometers and five electrical sensors as well as a signal marking the time of source impact. The electrical signals were received on outcrops by brass electrodes in the form of plates measuring 15x10x0.4 cm, and in sediments by brass stakes 70 cm long and 2 cm in diameter. Volarovich et al. (1962) implied (but did not state explicitly) that the electrical signals they recorded were from dipole sensors which measured the potential difference between two electrodes. The electrode separation was determined experimentally and varied from 1 to 5 m. A seismometer was typically deployed at the same location as each dipole to facilitate comparison of the arrival times and magnitudes of the seismic and piezoelectric signals. Volarovich et al. (1962) reported that piezoelectric responses from quartz-rich gneisses with magnitudes between 5 and 25 p:V were discernable at distances up to 25 m from the impact point. These signals appeared to travel at approximately the same velocity (2000 m/s) as the seismic surface waves but tended to precede the surface wave arrivals by 1 15 to 2 ms. The dominant frequencies of both the piezoelectric and the seismic arrivals were in the range 200-300 Hz. The physical process thought to be responsible for host rock piezoelectric signals was described by Sobolev and Demin (1980, p. 26). They state that an electric polarization is generated at each point in the rock as it is stressed by a passing elastic wave. The resultant electromagnetic field is propagated ahead of the elastic wave but is so weak and/or so highly attenuated with distance that it is only detectable a short distance in advance of the elastic wave. The quartz veins investigated by Volarovich et al. (1962) were coarse grained and 0.3 - 1 m wide. Piezoelectric signals varying in magnitude from 20 to 60 pV were measured by dipoles located directly above an exposed quartz vein and 14 m from the weight drop impact point (Volarovich et al., 1962, Fig. 6). Simultaneous piezoelectric arrivals were recorded on dipoles at different locations within a few metres of a vein; the common arrival time was found to correspond to the expected time required for a seismic wave to travel from the impact point to the vein. Dipoles located at greater distances from the vein recorded piezoelectric signals from the quartz-rich host rock rather than the vein. Measurements made along a vein that tapered from 1 m to 0.05 cm in thickness indicated that the piezoelectric response diminished as the thickness decreased. Volarovich et al. (1962) also reported that they were able to record signals from veins covered by 1 to 1.5 m of overburden provided that the weight was dropped directly onto bedrock. Finally, it was reported that the dominant frequencies of piezoelectric signals from quartz veins were in the range 500 - 600 Hz while those of the seismic signals were only 200 - 300 Hz. This discrepancy in dominant frequency was not observed for piezoelectric signals emanating from the quartz-rich host rock. Volarovich et al. (1962) also observed piezoelectric signals from a pegmatite vein. The vein consisted of coarse-grained feldspar, quartz, and mica, and inclusions of quartz-16 biotite mica, and varied in thickness from 2 to 12 m. It displayed large scale heterogeneities which were correlated with variations in the magnitude of the piezoelectric response at different points along its extent. Piezoelectric signals from this thick pegmatite vein were generally comparable in magnitude to those from quartz veins but displayed no increase in the dominant frequency relative to that of the seismic arrivals. In 1955 investigators at the All-Union Scientific Research Institute for Methodology and Technology of Exploration (VITR) began to study the possibility that seismoelectric phenomena could be used in exploration geophysics. The results of theoretical, laboratory, and field investigations of the piezoelectric method performed at the VITR were summarized by Neishtadt et al. (1972). In the laboratory investigations, physical models composed of aggregates of piezoelectric crystals were used to try to simulate the responses of natural piezoelectric targets. The field studies resulted in the development of instrumentation, and of procedures for the acquisition and interpretation of data. Neishtadt et al. (1972, p. 42-51) also outlined several case studies involving use of the piezoelectric method at various mineral deposits within the Soviet Union. In most cases the method was used to locate ores associated with quartz veins or quartz-pegmatite deposits. However, the method was also used to locate ore zones rich in the piezoelectric minerals sphalerite and chalcopyrite in polymetallic sulphide mines (ibid, p. 48-49). During the 1950's and 1960's research and development of the PEM proceeded concurrently at the VITR and at the Earth Physics Institute of the Academy of Sciences of the USSR. M. P. Volarovich, E. I. Parkhomenko, and G. A. Sobolev were among the principal investigators at the Academy of Sciences. Several more organizations became involved with the development of the method after it was introduced as a production tool for exploration during the period 1962-63 (Kondrashev, 1980, p. 3; Neishtadt et al., 1972, p. 2). Several designs of instrumentation for the PEM were developed by various groups for use on 17 the surface, in boreholes, and in mines (Neishtadt et al., 1972, p. 21). By 1978 there were over 100 PEM stations in operation in various regions of the USSR on deposits of gold, tin, mica, rock crystal, and other ores (Volarovich, 1978). The most recent general reviews of the piezoelectric method in the USSR are the monographs of Sobolev and Demin (1980), and Kondrashev (1980). Recently completed English translations of part of the former and all of the latter were commissioned by the Geophysical Instrumentation Group (GIG) at the University of British Columbia. Both of these references include an introduction to the piezoelectric effect in rocks, descriptions of instrumentation and field technique, examples of application of the PEM, as well as descriptions of data interpretation methods and theoretical modelling. Theoretical models were treated in more detail by Sobolev and Demin while Kondrashev devoted more attention to data interpretation. According to Kondrashev (1980, p. 145-147), explosives of various types are the standard seismic source for piezoelectric exploration. Charges weighing up to 10 kg may be used in areas where elastic waves are subject to poor coupling and rapid attenuation (eg. in soft, loose rocks with compressional wave velocities vp < 1000 m/s) but smaller charges are normally used. In underground conditions, typical charge sizes are 0.15-2 kg. Piezoelectric signals are normally sensed using electrodes to measure electric potentials. The amplitudes of piezoelectric signals observed during field studies in the USSR range from tens of microvolts (Volarovich et al., 1962) to a maximum of ten millivolts (Kondrashev, 1980, p. 42). These amplitudes refer to the potential differences that existed in the earth between two electrodes separated by various distances normally ranging from a few metres to several tens of metres. Sobolev and Demin (1980, p. 35-37, p. 141) described some experiments for which the receivers included electric and magnetic antennas but Kondrashev (1980) made no mention of such antennas. The PEM is most effective in its underground variants (in mines and boreholes) where the maximum range of investigation is typically 70-100 m (Kondrashev, 1980, p. 3; Sobolev and Demin, 1980, p. 129). However, a 18 limited number of underground experiments have suggested that the range of investigation for sphalerite-rich orebodies may be 300-500 m when large charges of up to 10 kg are used (Sobolev and Demin, 1980, p. 142). In the surface variant, the depth of exploration varies from only 3-5 m in areas where the quartz veins are situated in porous deposits with low seismic velocities (vp=300-500 m/s) to 30-40 m in more competent rocks characterized by relatively high velocities (eg. vp=3000 m/s) (Kondrashev, 1980, p. 4). 1.4(b) Development of the PEM outside of the USSR Outside of the USSR, the PEM has received relatively little attention and the author is not aware of any cases in which it has been used as a production tool for exploration. Ignorance of the method in the west is largely due to the fact that the Soviet experience has been poorly communicated. Some reports on the PEM have been translated from the original Russian but few are widely available. Furthermore, descriptions of the instrumentation are vague, while explanations of field techniques and plots of survey results typically lack clarity in the Soviet literature. Nonetheless, the claims of Soviet researchers regarding the usefulness of the PEM have encouraged some western researchers to investigate the method experimentally. The author is a member of a research group at the University of British Columbia (UBC) that is investigating the viability of the PEM with the support of a three year grant from NSERC and five industrial partners. Interest in the piezoelectric method at UBC was initiated during the late 1970's when Cominco Inc. supported laboratory and field tests conducted by B. Jose (Jose, 1979). The piezoelectric effect was measured in samples of quartz vein and host rock from the Con Mine in the N.W.T. and two field tests were conducted at the same mine. Jose (1979) claimed that some of his field results showed piezoelectric responses but his studies were not sufficiently comprehensive to be convincing. 19 In 1983 the Soviet scientists G. A. Sobolev and V. M. Demin, hosted by researchers from UBC, came to Canada to demonstrate the piezoelectric and PRRER2 exploration methods (Sobolev et al., 1980). Field tests were conducted at the Giant Yellowknife Mine, N.W.T. - where gold is associated with quartz lenses - and at the Sullivan iron-lead-zinc sulphide mine near Kimberley, B.C.. In the Giant Yellowknife Mine, piezoelectric signals from quartz lenses were measured using a single magnetic antenna with a bandwidth of 1-20 kHz. These signals implied that the magnitude of the electric field in rock in the vicinity of the antenna was about 0.2 mV/m for a typical layout having the shotpoint and antenna at distances of 30 m and 300 m respectively from the quartz lens. Two other antennas with bandwidths of 20 kHz to 3MHz were used to measure high frequency magnetic signals. These were interpreted as PRRER signals generated by thin layers of pyrite, arsenopyrite, and other natural semiconductors found within and around the quartz lenses. In the Sullivan mine strong PRRER signals were observed on the high frequency antennas after explosives were detonated at distances of about 30 m from an orebody containing beds of pyrrhotite, sphalerite, and galena. The strongest signals appeared to originate in the orebody. Despite the presence of the piezoelectric mineral sphalerite, no significant signals were sensed by the low frequency antenna. During the period June, 1985 to December, 1989 the UBC seismoelectric research group under the direction of R. D. Russell carried out seven field tests - four on the surface and three in mines - to investigate the piezoelectric and PRERR exploration techniques (Maxwell, 1988; Russell et al. 1990 NSERC Report). Some plausible piezoelectric and PRRER signals were observed. However the tests proved most useful for identifying sources 2 PRRER is an acronym for pulsed radio-wave range of electromagnetic radiation. See section 1.6 for a description of PRRER signals. of noise (particularly electrical noise associated with explosive detonation) and for guiding the development of instrumentation and field techniques. In 1987, A. Boyle, a research scientist with CRA Exploration Pty. Ltd., recorded piezoelectric signals at the Humboldt quartz vein site located 80 km northwest of Melbourne, Australia (Boyle, 1989, personal communication). The quartz occurs as numerous veins and stringers within sandstones, siltstones, and shales. Boyle's technique was novel in that he used a pneumatic breaker (similar to a jackhammer) for a seismic source and long wire electric field antennas for receivers. Each long wire antenna consisted of a 50 m length of insulated wire laid out along the ground and attached to one input of an amplifier; the other amplifier input was connected to a ground stake driven into the earth near one end of the insulated wire. The wire is capacitively coupled to electric potentials in the earth and the amplifier output is proportional to the average potential below the wire, the average being weighted according to the capacitance per unit length at each point along the wire (Russell and Maxwell, 1990 NSERC Report). Long wire antennas were used rather than grounded dipoles because the former were found to be more immune to shaking or contact phenomena during the passage of a seismic wave. The signal-to-noise ratio of the data was enhanced by stacking the responses from 500 breaker blows at each shotpoint. In 1990 samples of quartz from the Humboldt site were tested in the laboratory of the seismoelectric research group at UBC. The samples demonstrated very large piezoelectric effects indicative of longitudinal piezoelectric coefficients that were up to 12% of the piezoelectric modulus d n for an ideal quartz crystal (Maxwell et al., in preparation). Thus, on the basis of both field and laboratory studies the Humboldt site seemed to be worthy of further investigation. In April, 1990, members of the UBC group together with A. Boyle conducted a second piezoelectric experiment at the Humboldt site which reproduced and extended the results obtained during the 1987 test. Piezoelectric signals arrived 21 simultaneously on four long wire antennas at different locations. As the shotpoint was moved away from quartz targets the delay between the time of source impact and the reception of piezoelectric signal increased. The maximum shotpoint-to-target distance for which piezoelectric signals were observed is not known as the more distant shotpoints were located on overburden of unknown thickness which may have covered quartz veins that were not mapped. However, the maximum delay between the time of source impact and the arrival of a suspected piezoelectric signal was about 3 ms which suggests the maximum distance to the target was about 6 m given the typical seismic velocities measured at the site. The long wire antennas were located at distances ranging from 5 m to 30-40 m from the targets. Piezoelectric signals having peak-to-peak magnitudes of up to 25 U.V were recorded indicating that the maximum electric fields below the antennas were of the order of 1 pV/m. Modestly successful piezoelectric experiments were carried out by the UBC research group at the Fairview, B.C. quartz vein site in August, 1990 and April 1991. Piezoelectric signals were recorded by grounded dipole and quadrupole receivers while operating a pneumatic tamper source directly upon the vein. Quadrupole receivers - formed by subtracting one dipole from another - were used to reduce the level of powerline and telluric noise. Stacking was used to enhance the signal-to-noise ratio, with stacks of 1000 tamper blows being typical. The level of powerline noise was so high however that even stacked quadrupole signals acquired using 60 Hz notch filters were highly contaminated by harmonics of 60 Hz. A processing technique (described in Chapter 2) was therefore developed which removed powerline harmonics by subtracting sinusoids of the appropriate frequency, amplitude, and phase from the electrical sensor data. After processing, signals which appeared to be piezoelectric were clearly visible. The peak-to-peak magnitude of the signal on a 10 m dipole sensor located 5 m from the vein was approximately 7 pV. Therefore, as at Humboldt, the magnitude of the piezoelectric field was close to 1 pV/m. 22 In June, 1991, researchers from the UBC seismoelectric group performed their most convincing piezoelectric experiment to date at a surface site near the Paymaster mine in Timmins, Ontario (Russell et al., 1991 NSERC Report). The target was a steeply dipping quartz-ankerite vein, 2-3 m thick, located within a shear zone in basalt host rock.. The vein is part of a quartz ankerite system that hosts significant gold deposits. Using 0.23 kg sticks of high velocity explosive in shallow drill holes at distances up to 20 m from the vein, clear simultaneous piezoelectric arrivals were recorded on multiple grounded dipole and quadrupole receivers. The receivers were positioned at different locations 4 to 8 m from the vein. For shots located close to the vein, the piezoelectric signals recorded on 5 m dipoles had peak-to-peak magnitudes of about 5 mV representing electric fields of 1 mV/m. For the more distant shots electric field strengths of 0.1 - 0.2 mV/m were observed. Powerline noise was largely removed through the use of 60 Hz notch filters and the application of the sinusoid subtraction processing technique. During the experiment it was observed that A M radio broadcasts were being demodulated by the receiver preamplifiers and were appearing as noise in the frequency band of interest. A M radio pick-up was essentially eliminated and the signal-to-noise ratio of the data was greatly improved by putting capacitors across the amplifier inputs. The elimination of A M radio noise was an important step as this noise source had also been observed during previous experiments conducted by UBC researchers at Fairview and Humboldt. Plots of data obtained during the Paymaster experiment are presented in Chapter 2. Outside of the USSR there have been few attempts to investigate the piezoelectric method by means of theoretical or physical models. However, one useful theoretical model was developed by Russell and Barker (1991). They used it to make order of magnitude estimates of the strengths of electric and magnetic fields that could be generated by targets during piezoelectric surveys. The model consists of an insulating spherical body that becomes uniformly polarized within a homogeneous, isotropic, conductive whole space. 23 Using quasi-static approximations, Russell and Barker showed that, the strength of the electric field was proportional to the resistivity of the whole space and that, at the relatively low frequencies dominant in piezoelectric exploration, the electric field should be easier to detect than the magnetic field. They estimated that the electric and magnetic fields at a distance of 50 m from a piezoelectric body excited by a 3 kg explosive charge detonated at a distance of 10 m would be approximately 10 u,V/m and 0.4 pT respectively plus or minus an order of magnitude. In their calculation, they assumed that the resistivity was 1000 Qm and that the piezoelectric coefficient was 10% of the longitudinal modulus d n for a single quartz crystal. Given the natural variability of parameters such as the resistivity, the piezoelectric coefficient, and seismic wave coupling, the electric field estimate is in agreement with the magnitudes observed during both Soviet and western field experiments. The smaller signals obtained at Humboldt and Fairview are consistent with the fact that low energy seismic sources were used at those sites. In Chapter 3. the model of Russell and Barker is discussed in more detail and is used to simulate electric potentials generated when a spherical seismic wave polarizes a planar piezoelectric body. 1.5 Other Seismoelectric Effects Seismoelectric or mechanoelectric phenomena are mechanisms by which the production or transmission of seismic waves causes the production of electromagnetic waves. There are several phenomena that fall into the general category of seismoelectric or mechanoelectric effects. These include the following: (i) the piezoelectric effect (ii) the seismoelectric effect of the second kind or seismoelectric effect E (iii) pulsed radio-wave range electromagnetic radiation (PRRER) (iv) the seismic-electric effect or seismoelectric effect J 24 (v) blast-associated electromagnetic radiation (blast-associated EM) (vi) receiver or electrode shaking and contact phenomena The first four effects in the list above are generated by the interaction of seismic waves with various natural materials in the earth. Blast-associated EM signals are generated by a variety of mechanisms associated with explosions. Signals due to receiver shaking or contact phenomena may be generated when seismic waves pass by sensors of electric or magnetic fields. Each of these effects may or may not be observed during a piezoelectric field survey and it is important to be able to distinguish between piezoelectric signals and those caused by the other mechanisms. The piezoelectric effect has been described above. Brief descriptions of the other effects are given below. The seismoelectric effect of the second kind or E effect was discovered by Ivanov (1939). It is an electrofiltration effect observed in porous rocks that are at least partially saturated with liquid. In these rocks electric double layers are present meaning that charges of opposite sign exist on opposite sides of the rock-liquid interface. Under the influence of an elastic wave, ions of dominantly one polarity in the liquid phase move faster than those of the other polarity in the solid. The relative movement of ions is equivalent to an electric current which may be detected using electrodes or other receivers. If present, the E effect is normally observed at an electrode 2-10 ms before the seismic arrival. The frequency spectra of the seismic and electric E-effect signals are practically the same at any given point. More information on the E effect can be found in Parkhomenko (1971), Neishtadt (1972), and Kondrashev (1980). The seismic-electric effect was originally reported by Blau and Statham (1936). They injected a dc current into the earth and observed that the current was modulated during the passage of seismic waves. Thompson (1936,1939) suggested that the modulation was 25 caused by the variation of earth resistivity with elastic deformation. In the Soviet literature the seismic-electric effect is known as the seismoelectric effect J (Neishtadt, 1972, p.l). PRRER signals were discovered while carrying out geophysical surveys by the piezoelectric method in the vicinity of polymetallic (sulphide) ore bodies containing piezoelectric minerals (Sobolev et al., 1982). The dominant frequencies in the signals were 10-100 times higher than those in the elastic waves inducing them. At distances on the order of 100 m from the orebody the magnitude of a PRRER signal typically exceeded that of a piezoelectric response from the same target by a factor of 100-1000. Sobolev et al. (1982) hypothesized that PRRER signals arise in polymetallic orebodies as a result of the amplification of electrical impulses generated by the cracking of rock. Piezoelectric minerals within the orebody are polarized by the elastic wave and produce an electric field which combines with semiconductor minerals in the orebody to generate the required amplification. When a receiver such as an electrode, electric field antenna or magnetic field antenna is shaken, the magnetic flux (due to the earth's magnetic field) passing through it may change in such a way as to induce eddy currents within the receiver. Signals of this kind are called receiver or electrode shaking effects and may be generated when a seismic wave or air wave passes a receiver. In addition the changing stresses that occur at the contact between an electrode and the earth during the passage of a seismic wave may cause changes in the electrochemical conditions at the contact thereby causing fluctuations in the measured electric potential (Thompson, 1939); signals generated in this way are said to be due to contact phenomena. As described by Kondrashev (1980, p. 52-58), blast-associated E M signals are produced by several mechanisms; some of these are associated directly with the blast while 26 others are generated by the interaction of seismic waves with materials surrounding the explosive. Mechanisms associated directly with a blast include the expansion of thermally ionized explosion products, and (for explosions in an unconfined air space) the adiabatic compression of air at the front of the shock wave. For explosions in boreholes the expansion of thermally ionized products may be delayed from the time of the blast if the borehole is stemmed with materials that prevent the immediate escape of the products. Experiments described by Kondrashev (1980) have indicated that this delay can be up to 5-10 ms depending upon the length of the stemming, the degree to which it is packed, and the strength of its adherence to the bore hole walls. As seismic waves pass through the material surrounding the explosive they can produce electromagnetic signals by fracturing the rock (the triboelectric effect), by deforming electrical double layers and, in the case of an explosion within piezoelectrically active rock, by the piezoelectric effect. Blast-associated E M signals typically last 2-5 ms or longer and can be many times larger than the piezoelectric signals of interest (Neishtadt et al., 1972, p.32). Experiments have shown that the magnitude and duration of blast-associated EM radiation from charges in boreholes may be reduced by using certain types of explosives and stemming materials, by decreasing charge size, and by increasing the depth of the charge (Kondrashev, 1980, p. 52-58). Experimental results also suggest that the magnitude of an explosion-related signal is approximately proportional to R 1 / 3 where R is the distance between the explosion and the point of measurement (Sobolev and Demin, 1980, p. 41). Neishtadt et al. (1972, p. 32) reported that the dominant frequencies of explosive impulses were higher than those of piezoelectric responses and were typically in the 500 - 2000 Hz range. Blast-associated EM signals recorded in 1991 by the UBC seismoelectric group in two sulphide mines have been dominated by frequencies on the order of 105 Hz or higher. The higher frequencies recorded by the UBC group may indicate that different mechanisms dominated the explosive impulses observed by Neishtadt et al. or may simply be due to the wider bandwidth of the UBC records. 27 2 DATA ANALYSIS 2.1 Introduction This chapter is devoted to the first steps in the analysis of PEM data. These steps are the identification of piezoelectric signals and the enhancement of the signal to noise ratio by data processing. The characteristics of piezoelectric signals and various types of noise are presented first. Three processing methods - one involving stacking, the others involving filtering of certain frequencies - are then described. The stacking method uses a robust triggering algorithm to align piezoelectric records from multiple seismic impacts before they are stacked. The two filtering methods are a sinusoid subtraction technique used for the suppression of powerline noise, and a simple but flexible frequency domain low-pass filter. Real data examples are used to illustrate the characteristics of piezoelectric signals and the effectiveness of the processing techniques. 2.2 Characteristics of Signals and Noise in Data acquired bv the PEM 2.2(a) Characteristics of Piezoelectric Signals A piezoelectric signal that is recorded in the field generally displays the following characteristics: (1) It arrives simultaneously on receivers at different locations. (2) Its form or shape may differ on receivers having different locations or orientations relative to the target. (3) The delay between the time of source impact and the arrival of the signal corresponds to the minimum time required for a seismic wave to travel from the shotpoint to the target. 28 (4) Its amplitude spectrum is comparable to that of the seismic wave that induced it. Piezoelectric signals tend to arrive essentially simultaneously on receivers at different locations because they are electromagnetic waves travelling at very high velocities. The time required for an electromagnetic signal to travel between different receivers (typically separated by tens of metres or less) is too small to be detected given the accuracy with which arrival times can be determined. Thus, as stated by point (3) above, any delay between the time of source impact and a piezoelectric arrival can be attributed to the seismic traveltime between the source and the piezoelectric target. The rule of simultaneous arrivals may be violated if the receivers are located at different points in a weakly piezoelectric body of great extent. In this case, the electromagnetic signal at any point in the body is so weak that it can only be detected above the ambient noise when the elastic wave inducing it is close to a receiver. As a result, piezoelectric signals appear to travel between receivers with a speed that is roughly equivalent to the seismic velocity. Volarovich et al. (1962) observed piezoelectric arrivals of this type during an experiment conducted over quartz-rich gneissic bedrock. The similarity between the amplitude spectra of elastic and piezoelectric signals may be attributed to the linear nature of the piezoelectric effect. As indicated by equation (1.1), the polarization acquired at any point in a piezoelectric medium is linearly related to the seismic stress at that point. It is of interest to note, however, that Sobolev and Demin (1980) claim that piezoelectric responses from some targets are enhanced at certain frequencies corresponding to the natural frequencies of oscillation of the target as a whole and of individual sub-units within the target. In such cases the amplitude spectra of piezoelectric signals have resonant peaks that are not present in the seismic spectra. In a typical PEM survey, seismic signals recorded by geophones are dominated by frequencies below 1 kHz. 29 If natural frequencies of oscillation are significant, resonance of the target as a whole may enhance the piezoelectric spectra at certain frequencies of up to several kHz while resonances of individual sub-units are strongest at frequencies below 25 kHz4 (Sobolev and Demin, 1980, p. 43-49). The transfer function relating a piezoelectric signal to the elastic wave that induces it is examined in more detail in Chapter 3. 2.2(b) Sources of Noise There are many types of electrical noise that can obscure or completely overwhelm piezoelectric signals. The main sources of noise that may be encountered during field experiments are listed below: (1) powerline harmonics (2) A M radio transmissions (3) leakage, crosstalk, or radiation from geophones or other sensors (4) telluric noise (5) other seismoelectric effects including - blast-associated EM - receiver shaking and electrode contact phenomena - the seismoelectric effect E - PRRER signals - the seismic-electric effect (modulation of earth currents) During most of the piezoelectric field experiments conducted by UBC researchers, powerline harmonics have been the noise of highest amplitude. It is not uncommon to find that the ambient electric field in the earth at 60 Hz is in the range of 5-10 mV/m peak-to-4 Communicated by R.D. Russell and M. Maxwell of the UBC seismoelectric research group who spoke with V.M. Demin in Moscow in August, 1991. 30 peak. Such amplitudes are 5 to 1000 times larger than the amplitudes typically expected of piezoelectric signals. Harmonics of 60 Hz are weaker than the fundamental but several of them are typically larger than or comparable in magnitude to piezoelectric responses. Some mines in Ontario use 25 Hz power in addition to 60 Hz with the result that both fundamental frequencies and their harmonics can contaminate electrical recordings. During data acquisition, powerline harmonics can be attenuated by the use of notch filters and quadrupole receivers. Soviet investigators tend to record using high-pass filters with corner frequencies in the range of 200-500 Hz (Neishtadt, 1972; Sobolev and Demin, 1980; Kondrashev, 1980) presumably because powerline noise is particularly strong at 50 Hz (the fundamental frequency in the USSR) and its first few harmonics. Data processing methods such as the sinusoid subtraction technique described in section 2.3(b) may also be used to suppress powerline noise. During recent piezoelectric field experiments carried out by UBC researchers, 60 Hz notch filters have been used to remove the largest noise component and the remaining harmonics have been removed in processing by the sinusoid subtraction method. The notch filters prevent 60 Hz noise from limiting the dynamic range of the instrumentation. The processing technique eliminates the need for notch filters at many different frequencies and effectively suppresses powerline harmonics without attenuating or distorting other frequencies. UBC researchers have observed that the preamplifiers used for grounded dipole and long wire antenna receivers are capable of demodulating A M radio transmissions; the demodulated transmissions appear as noise some of which lies in the same frequency band as piezoelectric signals. Provided that it is not obscured by powerline harmonics, A M radio noise can be identified by feeding the amplified channel into headphones and listening for voices, music or other indications of broadcasts. The UBC group found that the reception of A M radio by grounded dipoles could be essentially eliminated by putting a small capacitor across the inputs of the preamplifier. The size of the capacitor was chosen so that it would 31 act aa shunt for A M carrier frequencies (hundreds of kilohertz) but would not adversely affect the performance of the preamplifier in the frequency range of piezoelectric signals. No references to A M radio noise have been found in the Soviet literature. Electrical noise due to leakage, radiation or crosstalk from seismic channels is of particular concern because it could be mistaken for a piezoelectric signal. The form of such noise is similar to that of the seismic signal, although minor differences in arrival time and signal shape may be observed due to differences in the frequency response of seismic and electrical channels. Crosstalk noise is typically characterized by its identical form on all affected channels. It is caused by induction or by currents leaking from one channel to another within multi-channel instrumentation such as amplifiers, filters, and digitizers. Electromagnetic radiation from seismic sensors and the leakage of seismic sensor currents into the earth differ from crosstalk in that noise from these sources is received by the electrical sensors. The amplitude of radiation or leakage from a seismic sensor decrease as the distance between the electrical and seismic detectors is increased. Seismic signals may also radiate from poorly shielded cables. The standard test used by the UBC seismoelectric group to check for radiation, leakage, and crosstalk involves tapping the geophone (or other seismic sensor) after it and all other sensors have been laid out for an experiment. If any signals similar in character to that of the geophone channel appear on the electrical channels, geophone leakage, radiation, or crosstalk is suspected. Telluric noise is caused by natural currents in the earth, and by natural fluctuations in the atmospheric magnetic field due to lightning, magnetic storms, polar aurora, etc. (Neishtadt et al., 1972, p.35; Kondrashev, 1980, p.61). It occurs at random, may appear as a sudden impulse or as pulse packets lasting for 30-50 ms and longer, and is typically stronger in the latter half of the day (Neishtadt et al., 1972, p. 35). Telluric noise may be 32 distinguished from piezoelectric signals by its shape or by its failure to appear at the same delay time when an experiment is repeated. Although other seismoelectric effects can provide information about the earth, they can also interfere with or look similar to piezoelectric effects and are therefore sources of noise for the PEM. The other effects have been described in section 1.6. The ways in which they can be distinguished from the piezoelectric effect during PEM surveys are outlined here. When blast-associated E M signals are observed, they usually occur coincident with the time of the blast and last 2-5 ms (Neishtadt et al., 1972, p. 32). They tend to be dominated by higher frequencies than piezoelectric signals. Neishtadt et al. (ibid) claimed that the dominant frequencies of explosive impulses were in the range 500-2000 Hz. During recent tests conducted by the UBC seismoelectric group in sulphide mines, the blast-associated EM signals were dominated by frequencies of the order of 10 kHz or higher. Noise due to receiver shaking or electrode contact phenomena will not be observed before the seismic wave arrives at the receiver. Unlike piezoelectric signals, therefore, this type of noise will not arrive simultaneously on receivers at different distances (or more correctly at different traveltimes) from the shotpoint. Furthermore, a signal that appears on a receiver before the seismic arrival at that receiver cannot be attributed to receiver shaking or contact phenomena. The E effect may be observed in the vicinity of wet and porous rocks or sediments. E effect signals differ from piezoelectric responses in that they do not appear simultaneously on receivers at different locations. Rather, they are typically detected by a receiver 2-10 ms before the arrival of the seismic wave and appear to travel between receivers at speeds comparable to seismic velocities (Neishtadt et al., 1972, p. 34; Kondrashev, 1980). When 33 grounded dipole receivers are used, the E effect is generally maximized by aligning the two electrodes with the shotpoint and minimized by placing the two electrodes equidistant from the shotpoint (Neishtadt et al., 1972, p.34). The magnitudes of piezoelectric signals are not correlated in this way with dipole orientation relative to the shotpoint. Finally, Neishtadt et al. (ibid) also claim that the dominant frequencies in E effect signals are lower than those in piezoelectric responses. PRRER signals are of much higher frequency and magnitude than signals of piezoelectric origin. They are dominated by frequencies that range from tens of kilohertz to a few megahertz and their magnitude typically exceeds that of a piezoelectric response from the same target by a factor of 100-1000 (Sobolev et al., 1982; Sobolev et al., 1984). PRRER signals also tend to have a pulse-like character with abrupt onsets while piezoelectric responses are quasi-sinusoidal in form (Sobolev et al., 1984). The seismic-electric effect described by Blau and Statham (1936) and by Thompson (1936, 1939) was generated by using a seismic source to modulate a dc current that was being actively injected into the earth. If natural telluric currents in the earth are strong enough, it is possible that their modulation by seismic waves could also be detected and could be a troublesome source of noise for the PEM. Like piezoelectric signals, seismic-electric responses would arrive simultaneously on antennas at different locations and would be dominated by seismic frequencies. However, the author has not found any investigations in the literature of telluric current modulation by seismic waves. Neither Neishtadt et al. (1972) nor Sobolev and Demin (1980) nor Kondrashev (1980) list the seismic-electric effect (or the seismoelectric effect J as they call it) as a source of noise for the piezoelectric method. 34 2.2(c) Example: Piezoelectric Data from the Paymaster Quartz Vein Site A description of the Paymaster piezoelectric experiment performed by the UBC seismoelectric group in June, 1991 can be found in section 1.5(b). The target was a steeply dipping, 2-3 m thick quartz-ankerite vein. Figure 2-1 is a plan view of the most successful experimental layout used at the site. Four grounded dipoles were located at distances of 4-8 m south of the southern edge of the vein. Dipoles 1,2, and 3 were each 5 m long. Dipole 4 was 4 m long and oriented at right angles to the others. All shotpoints were located within or north of the vein with the most distant being SP7 at a distance of 20 m from the vein's northern edge. A shotpoint geophone planted 2 m from the top of each shothole was used to trigger the recording of each piezoelectric record. The seismic traveltime between each shotpoint and its shotpoint geophone was measured in a separate test and was used to calculate the actual time of the explosion (to within 1 or 2 ms) on the piezoelectric record. Figure 2-2 shows the record acquired at SP3.5E approximately 7.0 m north of the vein. The traces labelled SPG and GEO come from geophones located 2 m from the shotpoint and at the eastern electrode for Dipole 1, respectively. Those labelled D1-D4 are from Dipoles 1-4 located as shown in Figure 2-1. All dipole channels were acquired using a 60 Hz notch filter and have been processed by sinusoid subtraction to remove powerline noise at 25 Hz and at harmonics of 60 Hz. The bandwidth of the electrical channels was 16 to 1000 Hz. The frequency response drops off at 6 dB/octave in the region 1-2.9 kHz and at approximately 18 dB/octave above 2.9 kHz. P A Y M A S T E R S I T E + 7 N + 6 + 6E 5 m + 4 + 4.5E + 3A + 2.5 + 3.5E QUARTZ-ANKERITE DIPl -x+GEO * DIP2 DIP3 DIP4 F i g u r e 2-1. P l a n v i e w o f the Paymaster site near T i m m i n s , O n t a r i o s h o w i n g the locations o f f o u r g r o u n d e d d i p o l e s ( D I P l - 4 ) a n d thirteen shotpoints. A shotpoint geophone was placed two metres away f r o m each shot hole and a second g e o p h o n e was p l a c e d at the eastern g r o u n d stake for D i p o l e 1. T h e quartz-ankcritc v e i n d i p p e d steeply to the north and had been e x p o s e d by trenching. PAY60 SP3.5E SOURCE IMPACT 45.54 mS 5.0 V 400. mV 2.0 mV 2.0 mV 2.0 mV 2.0 mV o. 26. 52. 78. 104. 130. TIME (ms) Figure 2-2. A typical piezoelectric record from the Paymaster experiment acquired at SP3.5E seven metres north of the vein using four grounded dipole sensors (DIP1-4) laid out as shown in Figure 2-1. The traces labelled SPG and GEO come from geophones located two metres from the shotpoint and at the eastern electrode for Dipole 1 respectively. A piezoelectric signal appears simultaneously on the four dipoles approximately 2.9 ms after the time of explosion indicated by the dashed line. A scale bar to the left of each dipole trace indicates the voltage that existed in the ground between two electrodes spaced four to five metres apart. Powerline harmonics have been removed from the dipole traces using the sinusoid subtraction technique. 37 The explosive source was detonated at 45.5 ms. Approximately 2.9 ms later a piezoelectric signal arrived simultaneously on all four dipoles. Note that the phase of the signal's first arrival is not the same on all four dipoles. This is likely due to the different location and orientations of the four dipoles relative to the part of the vein that was initially polarized by the seismic wave. The dominant frequencies in the signal are 200-400 Hz and are comparable to the dominant frequencies in the early part of the seismic signal labelled GEO. The lower frequencies that dominate the later part of this seismic trace may be related to surface wave arrivals. If the part of the vein closest to the shotpoint is responsible for the first part of the piezoelectric signal then the seismic wave must have travelled the 7 m between the shotpoint and the vein in 2.9 ms. This corresponds to an average seismic velocity of 2400 m/s which is not unreasonable for travel through the near surface. The average seismic velocity over the 13.6 m between the shotpoint and the geophone labelled GEO was comparable at approximately 2100 m/s. The magnitude of the piezoelectric response was strongest on Dipoles 3 and 4 where it reached approximately 2.5 mV peak-to-peak. The magnitude of the electric field was therefore approximately 0.5 mV/m at these sensors. As expected, the electric field strength at Dipole 2 - the sensor most distant from the vein - was considerably less at about 0.1 mV/m. The suspected piezoelectric signal cannot be attributed to electrode shaking or contact phenomena because it arrives simultaneously on the four dipoles and in advance of the seismic arrival at the eastern electrode for Dipole 1 (the trace labelled GEO). The simultaneity of the arrivals also indicates that they are not due to the seismoelectric effect E. It is unlikely that the signal is caused by geophone crosstalk, radiation, or leakage because its arrival time does not match either of the seismic arrivals and its shape is not consistently similar to either of the geophone signals. Furthermore, geophone tap tests of the kind 38 described in section 2.2(b) were performed at Paymaster and no signs of crosstalk, radiation, or leakage were observed. Blast-associated EM noise was not observed on any of the records obtained at Paymaster. Figure 2-3 shows the piezoelectric signals received by Dipole 1 for twelve shotpoints at various distances from the vein. Time zero is the moment of explosion and each trace is plotted at the distance corresponding to its north-south displacement relative to SPO (a shotpoint within the vein). As expected, the delay between the time of explosion and the arrival of a piezoelectric response increased regularly as the shotpoint was moved farther away from the vein. Also, with the exception of shotpoints 1 and IE which were within the vein, the signal magnitudes generally decreased as the shotpoint to vein distance was increased. All of the traces in Figure 2-3 were acquired using a 60 Hz notch filter and the bandwidth described for Figure 2-2. Harmonics of 60 Hz, and the 25 Hz fundamental were removed using the sinusoid subtraction technique and an 800 Hz low-pass filter with a slope of 36 dB/octave was applied in processing. Figure 2-3. (The figure, adapted from Russell et al., 1991, is on the following page.) Demonstration of signal moveout at the Paymaster experiment. The piezoelectric arrivals sensed by Dipole 1 are shown for twelve shotpoints located as shown in Figure 2-1. Each trace is plotted at the distance corresponding to its north-south displacement relative to SPO. The shotpoint number for each trace is superimposed on the distance axis. Source detonation occurs at time zero. Note that the delay between the time of detonation and the arrival of the piezoelectric response increases regularly with the distance between the shotpoint and the vein. 40 2.3 Processing Techniques for Improving the Signal-to-Noise Ratio of PEM Data 2.3(a) Stacking Using a Robust Triggering Algorithm Stacking and Triggering Considerations The UBC seismoelectric group has used low energy, repeatable seismic sources such as tampers during some piezoelectric experiments; the responses from multiple (up to 3000) tamper blows at a single shotpoint have been stacked to improve the signal-to-noise ratio (S/N). In the presence of random noise, the stacking of n common signals would be expected to improve S/N by a factor of n1'2 (Sheriff and Geldart, 1982, p. 126). For best results, the responses from individual tamper blows must be properly aligned so that they add constructively during stacking. This may be accomplished by ensuring that the time of seismic source impact is the same on all pre-stack data traces. A trigger signal, typically provided by a geophone or accelerometer located close to the shotpoint, is used to indicate the time of impact. The effectiveness of stacking is therefore dependent on the quality of the trigger signal and the reliability of the triggering algorithm used to detect its onset. Figure 2-4 shows the first 750 ms of pre-stack data acquired during a stacking experiment at the Fairview, B.C. quartz vein site in August, 1990. The uppermost trace in Figure 2-4 is from a trigger geophone and shows signals from six separate tamper blows. The geophone was positioned seven metres away from a point on the quartz vein where the tamper was being operated. The other traces are from two grounded quadrupole sensors that were used to measure piezoelectric responses. Powerline harmonic noise dominates the pre-stack quadrupole data. In order to stack responses from separate tamper blows a triggering algorithm must be used to extract a window of data (on all three channels) around each B37 SP1: PRE-STACK DATA TIME (ms). . Figure 2-4. The first 750 ms of data recorded on three channels during a stacking experiment at the Fairview, BC quartz vein site. A tamper was used to pound directly upon the vein approximately 10 times per second. The top channel is from a trigger geophone located 7 m away from the shotpoint and shows signals from six tamper blows. The other channels are from two grounded quadrupole sensors and are dominated by powerline harmonics. Scale bars indicate the size of the signals received by the sensors prior to any amplification. Figure 2-5 shows the sensor layout. ^ 42 trigger signal. The onset of the trigger signal should be the same within each window to ensure optimal stacking. Most of the piezoelectric data collected by the UBC seismoelectric group using tamper-like sources has been stacked using a computer program that is directly linked to the digital data acquisition system. The program, written by R.D. Russell, employs a simple trigger algorithm, based on signal level and slope, that allows it to be implemented in real time. Each record from an individual tamper blow is added to the stack immediately after it is acquired and no pre-stack data is saved. The trigger algorithm extracts a data record whenever the output of the trigger geophone rises above a certain positive voltage (or alternatively whenever the output drops below a certain negative voltage). The program has proven to be tremendously effective for stacking when the trigger geophone is located within about 0.5 to 2 m of the shotpoint. At these distances trigger signals tend to be well-defined with clear arrivals, relatively fast rise times, and no more than one or two large amplitude cycles. However, errors in triggering can still occur for a few reasons. For example, the amplitude of the trigger signal can vary considerably from one source impact to the next; the weaker the signal the longer it is likely to take to rise above the trigger level. Furthermore, some trigger signals may be so weak that the first break is missed entirely and triggering could occur on the second or third cycle of the signal. High amplitude noise may also cause triggering to occur at the wrong time. Another source of concern is the fact that jackhammers and tampers typically operate at 10-20 impacts per second while the real time stacking program can stack only about three records per second (for typical record lengths of about 100 ms). Thus, there is a chance that the algorithm will finish stacking one trace and begin looking for the next trigger point just after a source impact. In that case, the first break would be missed and triggering could occur on a later cycle in the signal. 43 When the seismic source is located very close to or on top of a piezoelectric target it is sometimes desirable to set the trigger geophone several metres away from the shotpoint. In this way a piezoelectric signal can be recorded on the electrical sensor channels before the arrival of the seismic signal on the trigger channel. Such experiments - called remote trigger tests- are done to demonstrate that the piezoelectric signal cannot be attributed to leakage, radiation, or crosstalk from a geophone channel. Unfortunately, as the distance between the shotpoint and trigger geophone is increased the trigger signal deteriorates; it becomes noisier, longer and more ringy and the sources of trigger error mentioned above typically become overwhelming for the real-time stacking program. The data in Figure 2-4 was acquired using a remote trigger geophone at a distance of seven metres from the shotpoint. Interest in effectively stacking remote trigger data was the major motive for development of the robust trigger algorithm. The algorithm was originally developed by the author as part of a course project for Geophysics 514 (Advanced Time Series Analysis) taught by Dr. T J . Ulrych. A description of this project can be found in a report written by the author (Butler, 1990). The Robust Trigger Algorithm The robust trigger algorithm determines how to align traces for stacking by cross-correlating the envelope of trigger signal k with the envelope of a reference trigger signal. Signal k is then shifted forward or backward in time to the position where the correlation is highest. A FORTRAN program has been written to implement the algorithm and is listed in Appendix A. Due to its computational complexity, the program is not fast enough to stack data in real-time. Pre-stack data such as that shown in Figure 2-4 must be stored in a file, and stacked at a later time. A step by step description of the program follows. Triggering and stacking is controlled by several user-specified parameters. First of all, the user must specify the location within the data file of the signal to be used as the 44 reference or pilot trigger trace. This requires the user to plot the first few signals on the trigger channel and choose one that looks typical. Next, a trigger level (LEVEL) must be specified; the level should be high enough that it is above the noise but low enough that the reference trigger signal rises above it. The total number of samples wanted in each trace (NFT) and the number wanted in the pre-trigger (PRE) must also be specified. The pre-trigger is the portion of the trigger signal that precedes the trigger point (the point where the signal first exceeds the trigger level). After calculating the envelope of the reference trigger signal, the program begins searching for trigger points in the pre-stack data. A 'rough' trigger point is located where the data on the trigger channel first rises above the specified LEVEL. A trace consisting of NFT samples including a pre-trigger of PRE samples is extracted about this 'rough' trigger point. The next step is to refine the initial trigger point selection by cross-correlating the reference envelope with the envelope of the extracted trace. The trigger point is shifted to align the two envelopes as closely as possible (ie. to obtain the maximum correlation between the envelopes). Traces on all data channels are then extracted about the refined trigger point and are added to their respective stacks. The program skips over a user-specified number of samples (SKIP) following the refined trigger point before it begins searching for the next trigger. This process of picking a rough trigger point, refining it by cross-correlation, extracting traces and stacking them is repeated over and over until the program has searched through as much of the data file as specified by the user. The final stacks on each channel are then normalized by dividing by the number of traces stacked and written to an output file. In practice the pre-stack input data file is usually too large to be held all at once in the memory of a personal computer. For example, a typical file contains enough data for a stack of 1500 traces on four channels. Typical sample rates and source repetition rates are 0.0003 45 seconds and twelve impacts per second respectively and all samples are stored as two byte integers. The size of such a data file would be 1500 samples .12x0.0003 channel x 4 channels x 2 bytes = 3.33 x 106 bytes. sample The program therefore reads and operates on only one swath of data at a time. Stacks of all traces within a given swath of data - called intermediate stacks - are also written to an output file. The speed of the triggering algorithm could be increased by omitting the calculation of envelopes and simply cross-correlating the trigger signals themselves. However, two very ringy trigger signals (such as those shown in Figure 2-4) can appear to be well correlated even when they are a cycle or two out of alignment. The envelope of a ringy signal is much smoother than the signal itself. As a result the lag or shift for which two envelopes are best correlated is better defined. Envelopes are therefore used in order to enhance the algorithm's ability to deal with poor trigger signals. The envelope E(t) of a real function g(t) is defined as the magnitude of the associated analytic signal (g(t) - /GHi(t)} where i= (-1)1/2 and GHi(t) is the Hilbert transform of g(t). The analytic signal can be obtained by taking the Fourier transform of g(t), zeroing the negative frequency components, doubling the positive frequency components, leaving the DC component unchanged, and transforming back to the time domain (Bracewell, 1986). A detailed derivation of this result can be found in a report written by the author (Butler, 1990). The stacking program uses a FFT (Fast Fourier Transform) and the procedure described above to find the analytic signal and calculate its (time varying) magnitude E(t) whenever the envelope of a trace is required. The cross-correlation of two functions g(t) and h(t) is defined as 46 corr(g,h)(t) = J g ( x + t)h(x)dt . -OO The cross-correlation of two real functions may also be obtained from the relation corr(g,h)(t) = J-i{G(co)H*(o))} where G(co) is the Fourier transform of g(t) H*((o) is the complex conjugate of the Fourier transform of h(t) and Jx{ } denotes the inverse Fourier transform. When two time series are to be cross-correlated at many lags it is more efficient to make use of the FFT and perform the calculations in the frequency domain. Trigger signal envelopes are cross-correlated in this way using subroutine CORREL from the Numerical Recipes library (Press et al., 1989). In using the FFT to facilitate cross-correlation in the frequency domain it is implicitly assumed that the envelopes are periodic with a period equal to the length of the envelope time series. The assumption is false and will introduce error into the cross-correlation at all lags other than zero. Thus, if correlations for lags as large as ±n sample intervals are desired then the last n samples in each envelope time series are set equal to zero. In this way, relatively little error is introduced into the cross-correlation provided that the samples near the beginning and end of the envelope time series are already tending toward zero. In the stacking program, n is a user-specified parameter called MAXLAG. When choosing this parameter the user should keep in mind that the last MAXLAG samples in all envelope traces are replaced with zeros. If MAXLAG is made too large , important information regarding the shape of each envelope may be lost. On the other hand, if the maximum correlation between two envelopes occurs at a lag greater than MAXLAG, no trace shifting is done; rather the program discards the current trace and begins searching for the next trigger point. Thus, if MAXLAG is made too small, valuable traces may be discarded. 47 Real Data Example The performance of the robust trigger algorithm will be illustrated using the data set shown (in part) in Figure 2-4. As mentioned previously, this data was acquired using a tamper source directly upon the quartz vein at Fairview, B.C.. A plan view of the experimental layout is shown in Figure 2-5. The vein is about 1.5 m wide and dips steeply NNE. The trigger geophone was located seven metres away from the shotpoint. Both grounded quadrupole receivers had dimensions of 10 x 10 m. The distances from the shotpoints to the nearest receiver electrodes were 4.5 m and 9.5 m for Quadrupoles 1 and 2 respectively. Figure 2-6 shows the reference trigger signal (top), followed by the first six signals on the trigger channel as they were aligned for stacking. The lowermost trace is the final average stack of 2504 signals from the trigger channel. It is clear that the trigger algorithm has done a good job of aligning all but one of the pre-stack signals despite their poor signal-to-noise ratio, varying amplitudes, and ringy character. The fact that the final stack exhibits little high frequency attenuation is further evidence of good pre-stack signal alignment. The misalignment of the first break on the trace labelled T3 can be attributed to its anomalous shape. The final average stacks of 2509 traces on all three channels are plotted in Figure 2-7. The vertical dotted line indicates the approximate time of source impact (± 1-2 ms). This was established by measuring the seismic traveltime from the shotpoint to the trigger geophone and subtracting it from the arrival time of the trigger signal. The stacked quadrupole signals are still dominated by powerline noise. However, stacking has reduced 48 the noise level - as indicated by the scale bars in Figures 2-4 and 2-7 - by factors of approximately 40 and 60 on Quadrupoles 1 and 2 respectively. In order to remove the powerline harmonics, sinusoid subtraction processing was applied to the stacked quadrupole traces. The processed data is shown in Figure2-8. A signal that appears to be piezoelectric is now clearly visible at the time of source impact. The signal arrives simultaneously on both sensors. It cannot be attributed to geophone crosstalk, radiation, or leakage because it arrives well in advance of the trigger signal. Electrode shaking and/or contact phenomena may be responsible for some of the large amplitude excursions that occur several milliseconds after the time of source impact. The magnitudes of the suspected piezoelectric signal are about 12 and 6 pV peak-to-peak on Quadrupoles 1 and 2 respectively. In contrast, the noise levels in the pre-stack data (as shown in Figure 2-4) acquired by Quadrupoles 1 and 2 were approximately 1.2 and 3 mV peak-to-peak. The noise in the pre-stack data was therefore approximately 100 times larger than the signal on Quadrupole 1 and 500 times larger than the signal on Quadrupole 2. 49 10 T R I G G E R G E O P H O N J E QUAD 2 QUAD 1 Figure 2-5. Plan view of the layout used to acquire data set B37 at the Fairview, BC quartz vein site. A tamper was used to pound directly upon the vein at the location marked SP1. A remote geophone located 7 m away from SP1 provided the trigger signal for stacking. Two grounded quadrupole receivers with dimensions of 10 x 10 m measured the piezoelectric response. The vein dips steeply NNE. B37 TRIGGER SIGNALS 0. 34. 68. 102. 136. 170. TIME (ms) Figure 2-6. Trigger signals extracted from data set B37 and aligned for stacking by the robust triggering algorithm. The reference signal and final stack of 2504 signals are the uppermost and lowermost traces respectively. The traces labelled TI - T6 are the first six trigger signals in the data set and can also be seen, prior to extraction, in Figure 2-4. B37 SP1: RAW STACKS OF 2504 SOURCE IMPACT 55.40 mS TIME (ms) Figure 2-7. Mean stacks of seismic and electric signals generated in response to 2504 tamper blows upon a quartz vein at the Fairview, BC site. The experimental layout is shown in Figure 2-5. The uppermost trace is from the trigger geophone and the others are from two grounded quadrupole receivers. Powerline noise dominates the quadrupole data and obscures any piezoelectric signal that may be present. The dotted line at 55.4 ms indicates the time of source impact. B37 SP1 (AFTER SINUSOID SUBTRACTION) SOURCE IMPACT 55.40 mS 0. 34. 68. 102. 136. 170. TIME (ms) Figure 2-8. The result of subtracting powerline harmonics from the stacked quadrupole traces shown in Figure 2-7 is illustrated in this plot. A signal that appears to be piezoelectric is now clearly visible on both quadrupole channels at the time of source impact. The signal cannot be attributed to geophone crosstalk as it appears before the seismic arrival on the trigger channel. The amplitudes and phases of the powerline harmonics were estimated using the first 50 ms of each quadrupole trace. 53 2.3 (¥) Sinusoid Subtraction A processing technique involving the subtraction of sinusoids from piezoelectric records has been developed for the suppression of powerline noise. This type of noise is known to occur at the fundamental frequency of power transmission (typically 60 Hz in North America) and at frequencies that are harmonics of the fundamental (eg. 120 Hz, 180 Hz, 540 Hz, etc.). The amplitude and phase of the sinusoid to be subtracted at each powerline harmonic frequency are estimated by calculating its Fourier series coefficients. The success of the method is dependent upon the amplitude and phase of each powerline harmonic remaining constant over the record length. This condition is usually satisfied and the technique has proven to be highly effective in practice. The concepts of the Fourier series are presented in detail in many readily available mathematics texts (eg. Stremler, 1982). A brief review of the concepts relevant for sinusoid subtraction processing is presented below. Any real-valued, periodic function f(t) having finite energy per period may be represented by a trigonometric Fourier series. The representation is given by OO OO f(t) = ao + J^cos nG}ol + 5X s i n ncoGt n=l n=l 1 T where a* = ~ J f(t) dt 0 0 0 T is the period of f(t) 2TC and co0 = ~ . 54 Note that f(t) is represented by its mean value (aQ) plus an infinite sum of sines and cosines with frequencies that are positive integer multiples of 0 ) o . The fundamental frequency co0 is controlled by the period T. The sine and cosine terms for each frequency may be combined into a single cosine term of amplitude cn = yj&n2 + b n 2 and phase (|>n = arctant-bn/aj. This yields a more compact notation for the trigonometric Fourier series, oo f(t) = X C n C O S ( n 0 V + 4>n) • n=0 Suppose a piezoelectric record y(t) is contaminated by powerline harmonics having frequencies mcop (m = 1,2, 3,...). If the Fourier series representation is calculated for a portion y*(t) of the record that is dominated by powerline noise, the amplitudes and phases of the powerline harmonics may be approximated by the Fourier coefficients cn and <(>„ at frequencies mcop. In order to generate Fourier coefficients at these frequencies, the length of the extracted portion y*(t) must be a positive integer multiple of 27i/cop. In other words, y*(t) must contain an integral number of cycles of the frequency cop.5 Furthermore, it must be assumed that y*(t) is periodic with a period equal to its length (this assumption would be true if y*(t) was composed wholly of powerline harmonics). The presence of other types of noise or of piezoelectric signals in y*(t) influences the value computed for the Fourier coefficients. This limits the accuracy with which the amplitudes and phases of powerline harmonics can be estimated. A computer program has been written to implement the sinusoid subtraction technique on a DOS-based personal computer. The FORTRAN source code - called SLNSUB.FOR - and program documentation can be found in Appendix A. Data files 5 R.D. Russell has pointed out that this restriction would not apply if the values a and b for the cosine and sine components of a given frequency were determined using a general least squares formulation. The Fourier series coefficients are in fact the least squares solutions obtained for a and b in the special case for which the length of the data record is an integral number of cycles of the frequency of interest. 55 containing up to eight channels, with a maximum of 4096 samples per channel, may be processed. The user must specify which harmonics of 60 Hz are to be subtracted and may indicate that certain channels (eg. seismic channels) are to be written to the output file without undergoing sinusoid subtraction. The program is fast enough to be used in the field to process data sets immediately after they are acquired. The time required for a typical run in which 15 harmonics are subtracted from 4 channels of 4096 samples each is about 3 minutes on a 80286 laptop computer with a math coprocessor. A portion y*(t) of each piezoelectric record is extracted by the program for purposes of calculating Fourier coefficients at the frequencies specified by the user. This portion contains as many full cycles of the frequency cop as possible. If the total length of the time series is L seconds then y*(t) consists of the first T seconds where _ 2nk _ k % " fp fp is the fundamental frequency of power transmission (assumed to be 60 Hz) and k is the largest positive integer for which T < L. Note that the length of the piezoelectric record must be at least 27t/cop (0.01666... seconds for fp = 60 Hz) to permit calculation of the Fourier coefficient for Cup. The fundamental frequency u)0 in the Fourier series representation of y*(t) is given by co0 = 2n/T = u>p/k. The amplitude and phase of the powerline harmonic mcop are therefore approximated by the Fourier coefficients cn and (|>n where n = kxm. Once cn and <])n have been determined, the sinusoid cncos(nu)0t + <j>n) (or equivalently cncos(mu)pt + (|>n)) is subtracted from the full piezoelectric record thereby stripping off the powerline interference at frequency mtop. Real Data Examples Data set PAY59 in Figure 2-9 was acquired during the Paymaster piezoelectric experiment at SP4.5E using four grounded dipole sensors laid out as shown in Figure 2-1. 56 The four dipole traces are highly contaminated with powerline noise despite the fact that a 60 Hz notch filter was used during acquisition. This is not surprising since the site was located within approximately 200 m of two transmission lines supplying 60 Hz and 25 Hz to the Paymaster mine. The recording bandwidth was the same as that described in section 2.2(c) for the dipole traces in Figure 2-2. The result of processing data set PAY59 using the sinusoid subtraction program is shown in Figure 2-10. A piezoelectric signal, arriving about four milliseconds after the time of explosion is now clearly visible in the data from Dipoles 1 and 4. The signal is weaker but still discernable in the traces for Dipoles 2 and 3. A total of 26 powerline harmonics - 25 Hz and 25 harmonics of 60 Hz - were subtracted from each of the dipole traces. Note that the amplitude of the noise on Dipole 2 prior to sinusoid subtraction - at 13 mV peak-to-peak - was approximately 50 times larger than the underlying piezoelectric signal. Figure 2-11 illustrates the effect of successive sinusoid subtractions on the trace for Dipole 1. The original record and the final record remaining after the subtraction of 26 powerline harmonics are plotted at the top and bottom of the figure respectively. The intermediate traces show how the signal to noise ratio of the data improved steadily as harmonics were subtracted. Most of the improvement was obtained by subtracting the nine largest harmonics of 60 Hz. The amplitudes of the harmonics subtracted from Dipole 1 are plotted in Figure 2-12. The amplitude of the largest harmonic (540 Hz) is four times larger than the amplitude of the piezoelectric signal on Dipole 1. A second example of the performance of the sinusoid subtraction technique is provided by the stacked data set B37 discussed in section 2.3(a). Figure 2-7 shows the original stacked data. Figure 2-8 shows the data after the subtraction of 17 harmonics of 60 Hz from the quadrupole traces. The largest harmonics occurred at 60 and 180 Hz. The 57 standard sinusoid subtraction program SINSUB.EXE was modified in this case so that it used only the first 50 ms of the quadrupole traces to estimate the amplitudes and phases of the powerline harmonics. The data in this time window was recorded before the time of source impact and was dominated by powerline noise to a greater extent than that part of the trace recorded after source impact. However, as mentioned in the program description above, the standard sinusoid subtraction program was written to automatically extract as long a portion of the trace as possible for purposes of estimating the powerline noise; for data set B37, the first 166.7 ms (10 full cycles of 60 Hz noise) would be extracted. The result of using this wider time window to estimate the harmonics is illustrated in Figure 2-13. The suspected piezoelectric signal occurring at the time of source impact is visible but the signal to noise ratio is clearly better in Figure 2-8. Evidently the piezoelectric signal and the ringy noise (possibly caused by electrode shaking) visible after the time of source impact had an undesirably large impact on the amplitudes and/or phases estimated for powerline harmonics. This example demonstrates that, for purposes of estimating powerline harmonics it is desirable to record a significant length of data before the time of source impact or detonation. Furthermore, it would be beneficial for the user of the sinusoid subtraction program to have the option of specifying which part of the record should be used for estimating powerline harmonics. This option will be incorporated in a future version of the program. PAY59 SP4.5E (BEFORE SINUSOID SUBTRACTION SOURCE IMPACT 44.58 mS 0. 26. 52. • 78. 104. 130. TIME (ms) Figure 2-9. A raw (unprocessed) data set from the Paymaster experiment acquired at SP4.5E approximately ten metres north of the quartz vein. Four grounded dipole receivers (DIP1 - DIP4) were located as shown in Figure 2-1. The traces labelled SPG and GEO come from geophones located two metres from the shotpoint and at the eastern electrode for Dipole 1 respectively. The time of explosion is indicated by the dashed line. Any piezoelectric signals that may be present in the dipole data are hidden by powerline noise. A scale bar to the left of each dipole trace indicates the voltage that existed in the earth across two electrodes spaced four to five metres apart. PAY59 SP4.5E (AFTER SINUSOID SUBTRACTION) SOURCE IMPACT 44 .58 mS . 4.0 V 100. mV 1.0 mV 1.0 mV 1.0 mV 1.0 mV 26. 52. 104. 130. TIME (ms) Figure 2-10. The piezoelectric data set from SP4.5E at Paymaster after the suppression of powerline harmonics by sinusoid subtraction. A piezoelectric signal arriving approximately four milliseconds after the time of explosion can be seen on all four dipole receivers. Compare this plot to Figure 2-9 in which the same data set is plotted prior to sinusoid subtraction. SUCCESSIVE SINUSOID SUBTRACTIONS ON DIPOLE1 8.0 mV 4.0 mV 4.0 mV 4.0 mV 4.0 mV 4.0 mV 4.0 mV o > tr-1 o i n o 26. 52. 104, 130. TIME (ms) Figure 2-11. The effect of successive sinusoid subtractions on the trace recorded by Dipole 1 of data set PAY59. The original record is plotted at the top of the figure. One or more powerline harmonics were subtracted from each trace in order to generate the trace immediately below it. A total of 26 harmonics were subtracted to yield the final trace at the bottom of the plot. Note that the original record is plotted at half the scale of the others. 2.0 HARMONICS SUBTRACTED FROM DIPOLE 1 1.5 > CD zs 1.0 < 25 Hz .0 T t I T T 0 10 15 20 25 60 Hz Harmonic 30 35 Figure 2-12. Amplitudes of the 26 powerline harmonics subtracted from the data recorded by Dipole 1 of data set PAY59. All but one of the frequencies (25 Hz) are harmonics of 60 Hz. The amplitude of 60 Hz itself is small because a 60 Hz notch filter was used during acquisition. ON B37 SP1 (AFTER SINUSOID SUBTRACTION) SOURCE IMPACT 55.40 mS 0. 34. 68. 102. 136. 170. TIME (ms) Figure 2-13. Stacked data set B37 after the application of standard sinusoid subtraction processing. The amplitudes and phases of the powerline harmonics were estimated using most (the first 166.7 ms) of each quadrupole trace. A signal that appears to be piezoelectric is visible on the two quadrupoles at the time of source impact. However, the signal to noise ratio is not as high as that obtained by using only the first 50 ms of each trace to estimate powerline noise (see Figure 2-8). 63 2.3(c) Low-pass Filtering in the Frequency Domain Low-pass filtering is used to attenuate high frequency noise. Noise at frequencies above those present in a piezoelectric signal may be generated by telluric or atmospheric sources, by power lines, A M radio stations, recording instrumentation and other sources. Some low-pass filtering is performed during data acquisition due to the limited frequency response of amplifiers and recording instruments. In addition, analog time domain low-pass filters are normally used during acquisition to eliminate noise at frequencies well above those expected in the piezoelectric data and to prevent aliasing during digitization. Following acquisition it is sometimes desirable to enhance the signal-to-noise ratio by further low-pass filtering. The frequency domain filter described in this section is a flexible and effective tool for this purpose. A time series is filtered in the frequency domain by multiplying its Fourier transform by a transfer function and transforming the product back to the time domain. The general expression for frequency domain filtering of a time series g(t) is gs(t) = J1{J{g(t)}xS(o))} gs(t) = J-i{G(co)S(co)} (2.1) where gs(t) is the filtered time series to is radian frequency !F[ } and T1 { } denote the Fourier and inverse Fourier transforms G(co) is the frequency domain representation (Fourier transform) of g(t) and S(co) is the frequency domain transfer function. In general, S(cu) is a complex function of frequency and may be expressed as S(to) = a(co) + i'b(co) or S(o)) = | S(co) | e®«») 64 where |s(co)| =>/a(a))2 + b(to)2 is the magnitude of the filter at frequency ca and 6(00) = arctan(b(co)/a(o))) is the phase of the filter at frequency to. Equation (2.1) may therefore be re-written gs(t) = ^1{G(co)|s(co)U'^ )} (2.2). Equation (2.2) shows that, in general, a filter changes both the amplitude and the phase of all frequencies in the Fourier transform of the original time series. The application of a negative phase shift to a component in the frequency domain is equivalent to delaying that sinusoidal component in the time domain. Phase-shifting is therefore undesirable when the objective is simply to attenuate certain frequencies without otherwise changing the shape of a signal. Analog, time domain filters always introduce phase shifts; the shifts are greatest near the filter cut-off frequencies and outside the pass bands. Real-valued frequency domain filters however may be applied without affecting the phase of any frequencies. If S(co) is real-valued then b(co)=0 and the filter phase 6(u)) is zero for all frequencies. The lack of phase shifting is an advantage of filtering in the frequency domain rather than the time domain. A low-pass filter computer program has been written for use in processing digital seismoelectric data on a DOS-based personal computer. The FORTRAN source code, called LPF.FOR , and documentation explaining how to use the program are listed in Appendix A. A FFT algorithm (called REALFT) from the Numerical Recipes library (Press et al., 1989) is used to perform the necessary transformations between the time and frequency domains. Up to eight channels of data, each containing a maximum of 4096 samples, may be filtered. If desired, certain channels may be written to the output file without being filtered. The cut-off frequency fc (in Hz) and filter slope s (in dB/octave) are provided by the user and determine the filter transfer function S(f) according to the equation 65 S(f) = "—T~JT: where f represents temporal frequency in Hz. Since 1 + ( f / f c r / 0 ; S(f) is a real-valued function, the filter does not introduce phase shifts at any frequency. Note that S(fc) = 1/2 for any filter slope s. The cut-off frequency is therefore the 6 dB down point of the amplitude response of the filter. The transfer function is plotted for several different values of slope in Figure 2-14. It can be shown that as the filter slope is increased causing the amplitude response to roll off more sharply, the impulse response of the filter becomes more oscillatory or 'ringy'. In the limiting case of a response that drops from unity to zero at the cut-off frequency (ie. a boxcar filter), the impulse response is given by the sine function sin(27tfct)~ s(,)=2f<h^r which exhibits damped oscillation at the cut-off frequency fc (Stremler, 1982, p.l 15). Since the filtering of a time series corresponds to convolving it with the filter impulse response, strongly oscillatory impulse responses are undesirable. Thus, in the selection of filter slope a compromise is sought between large values that provide rapid roll off of the amplitude response with frequency and a low values that yield less oscillatory impulse responses. The result of processing the data shown in Figure 2-2 using the low-pass filter program is displayed in Figure 2-15. The cut-off frequency fc and slope s were 800 Hz and 24 dB/octave respectively. Note that the filtering has improved the signal to noise ratio of the dipole data without affecting the general shape or phase of the piezoelectric and seismic signals. L O W - P A S S TRANSFER FUNCTIONS CD--a 0 10 - 2 0 - 3 0 40 - 5 0 - 6 0 - \ V \ \ WW -- \ \ \ CP \ \ \ \ \ °~ \ \ \^ \ \ \ \ \ \ I l l \ - J \ \ \ \ t O \ \ \ l l I I I I \ l \ l \ I I \ I- I I I 10 -1 10 o 10 Relative Frequency ( f / f c ) Figure 2-14. Transfer function of the frequency domain low-pass filter for five different slopes ranging from 12 to 72 dB/octave. The phase response is zero (ie. no phase shift) for all frequencies. Note that the cut-off, fc, is the frequency at which the amplitude response is 6 dB down. ON ON PAY60 SP3.5E (AFTER SIN. SUB. and 800/24 LPF) SOURCE IMPACT 45.54 mS 5.0 V 400. mV 2.0 mV 2.0 mV 2!0 mV 2.0 mV "0 o o I—( 0. 26. 52. 78. 104. 130. TIME (ms) Figure 2-15. Data set PAY60 from the Paymaster experiment after processing by the sinusoid subtraction and low-pass filter programs. The same data is shown in Figure 2-2 prior to low-pass filtering. The filter cut-off frequency (6 dB down point) and slope were 800 Hz and 24 dB/octave respectively. 5 68 3 MODELLING 3.1 Introduction The analytic model described in this chapter estimates the electric potential generated in a conductive whole space by the passage of a spherical elastic wave through a piezoelectric body. Synthetic time series are generated in order to show how the electric potential varies as a function of both time and position. The main objective is to investigate the transfer function relating the elastic wave to the piezoelectric signal it induces. The model can also be used to make quantitative estimates of the magnitude and direction of the electric field created by a target possessing a specific piezoelectric fabric. The theoretical basis of piezoelectric modeling is described in section 3.2. The model developed by Russell and Barker (1991) for calculating the electric potential due to a polarized sphere is reviewed first. It provides a useful and intuitive approximation to the piezoelectric signal generated by a uniformly polarized target in conductive media. In order to examine the effect of non-uniform polarization, the target is next represented by many spheres each of which is independently polarized due to the stress associated with a spherical elastic wave. The expression for the electric potential generated by such a model is rather complicated as it depends on many factors including, for example, the curvature of the seismic wavefront, the dimensions of the target, and the form of the piezoelectric tensor. Consequently, a computer has been used to calculate synthetic records of the electric potential generated by 'typical' models thought to be reasonably representative of natural piezoelectric targets and actual piezoelectric field methods. For purposes of modelling it is necessary to specify the form of the spherical elastic wave and the piezoelectric tensor. A pulse-like particle displacement function, based on Sharpe's (1942) model, is used for the 69 elastic wave and the ideal fabric °°-m is used to define the piezoelectric tensor dj^ . Synthetic data generated for three different targets is presented in section 3.3. 3.2 Theoretical Basis 3.2 (a) Electric Potential due to a Polarized Sphere A first approximation to the transfer function relating seismic and piezoelectric signals was developed by Russell and Barker (1991). Their primary objective was to obtain order of magnitude estimates of the strengths of electric and magnetic fields that could be generated by natural piezoelectric targets. These estimates are described briefly in section 1.5(b) of this thesis. The target was approximated by a spherical, insulating body which acquired a uniform, time-varying electrical polarization P when stressed by a seismic wave. The sphere was surrounded by a conductive whole space of resistivity p. A spherical coordinate system oriented as shown in Figure 3-1 was used for the analysis. Using quasi-static approximations, Russell and Barker showed that the electric potential at a position r outside the sphere was given by the frequency domain expression ¥ ( r , c o ) = ¥ ( r , c o ) = 47ten 4rcen [/co + 2/(3peo)J m L/co + 2/(3peo)_ P(co) cos9 "p*(co) -7 (3.1) where V is the volume of the sphere e0 is the permittivity of free space and j = ypl. For typical resistivities of 105 Qm or less and for the frequencies of a few kilohertz or lower used in piezoelectric exploration, co « (2/3p£o) and equation 3.1 is well approximated by 3pVcos6 ¥(r,CO) » K 8 m 2 ;coP(co) (3.2a). 70 (r,e) Figure 3-1 (from Russell and Barker, 1991). Illustration of the spherical model used by Russell and Barker (1991) to represent a polarized piezoelectric body imbedded in a conductive medium. A spherical coordinate system is oriented with its z-axis parallel to the direction of polarization and its origin at the centre of the model. Figure 3-2. Diagram of the model used to estimate the electric potential induced at point (xj, x2, x3) by the interaction of a spherical elastic wave with a piezoelectric target of arbitrary shape. The target is approximated by a group of spheres. It is assumed that the total potential can be obtained by a linear superposition of the potentials generated by each polarized sphere. The origin of a rectangular coordinate system OxjX 2x 3 is located at the shotpoint. 71 Expressing this result in the time domain we have v 3pVcos6 8P dt (3.2b). Thus, the electric potential outside the sphere is proportional to the first time derivative of the sphere's polarization. Recall that polarization is linearly related to seismic stress by the piezoelectric tensor d^ according to the equation P; = d^a^. For plane waves, Russell and Barker point out that seismic stress is proportional to the product of particle velocity D and acoustic impedance. Thus, P = ku where k is a constant and the electric potential is given by 3pVcos0 ¥(r,<o) « k K 8 ? c r 2 yea v(co) (3.3a) or , » , 3pVcos9 ch) dt (3.3b). Equation (3.3a) is a simple transfer function relating seismic particle velocity to the electric potential expected outside a piezoelectric target. In the time domain, the electric potential is expected to be a scaled version of the first time derivative of the particle velocity at the sphere. This is a useful starting point for modelling the relationship between seismic and piezoelectric waveforms. However, the model is somewhat restrictive in that it assumes the target is spherical and that all parts of the target experience the same seismic stress as a function of time. Quartz veins, which tend to be planar in shape, area common target for piezoelectric exploration. The seismic wave will take longer to reach more distant parts of a target such as a vein and the resulting polarization will not be uniform across the target. Furthermore, shotpoints are usually located within several tens of metres of the target; at these distances the seismic wavefront may not be well approximated by a planar wave. Russell and Barker (1991) recognized the restrictions of their model and suggested that more realistic targets could be simulated by linear superposition of the potentials of 72 many spheres. Their suggestion has been followed in this chapter in order to simulate the electric potential due to a planar quartz vein. In addition, a piezoelectric fabric has been assigned to the vein, and a spherical elastic wave has been used to induce the polarization. 3.2(b) A More General Model A model for the PEM that incorporates spherical seismic waves and a target of arbitrary shape having a non-uniform polarization is illustrated in Figure 3-2. The seismic source is located at the origin of a rectangular coordinate system and n piezoelectric spheres are used to approximate the target. The spheres lie at different distances from the seismic source and acquire different time-varying polarizations. The electric potential generated by each polarized sphere is calculated using equation (3.1). In using this equation, it is implicitly assumed that the potential due to each sphere is not significantly altered by the presence of the other insulating spheres making up the target. This approximation is expected to be poorest immediately adjacent to an insulating target (where the current paths are likely to be most significantly altered) but to improve as the distance from the target is increased. The potential at point (xls x2, x3) is then obtained by summing the contributions from all n spheres in the target according to the equation n ¥(x! , x2 , x3, co) =-'47160 [;co + 2/(3p6o) 1 /co V/P /(C0)T' f (3.4) where r';= ((x, - xt/), (x2 - x2;), (x3 -x3/)). Note that equation (3.4) is a frequency domain expression. The polarization of sphere / is calculated using equation (1.1) in the frequency domain: (P c^o)), = d^fV, . (3.5) 73 The seismic stress cjjk is dependent on frequency (or time in the time domain), the distance r from the seismic source, and the type of particle displacement in the seismic wave. In order to reduce the complexity of the model, only compressional waves (p-waves) will be considered. The particle displacement for a spherical compressional wave is back and forth in the radial direction. It may be represented in the spherical coordinate system by u(r,t) in the time domain or U(r, co) in the frequency domain. For an isotropic medium, the seismic stress may be calculated from the particle displacement as shown in Appendix B. The resulting expressions for the normal stresses Ga and shear stresses o^  (i * j) are f ^ U 2X + 2p 1 - -5- + (X + 2\i)xi^ + X au au •> 9XJ kdx kj , i * j * k (3.6) u aa(r,a>)= = 3U 3U X i dXj + xJ8X i 2xiXjU where X and p are the elastic constants of the medium. It is assumed that the elastic constants and densities of the whole space and piezoelectric target are similar so that reflections and refractions of the seismic wave at the target boundaries may be ignored. The seismic stress is calculated at the centre of each sphere and is assumed to be uniform over the sphere's volume. This assumption is reasonable provided the diameter of the sphere is significantly less than one half seismic wavelength. Equations (3.4), (3.5), and (3.6) define a transfer function relating the particle displacement U(r,co) at n spheres to the electric potential at any point in the conductive whole space surrounding the target. However, the function is complicated and it is not easy to visualize the form of the electrical potential. In general, different spheres will become polarized at different times and it is expected that some will contribute constructively and others destructively to the potential at a given point. The relative importance of any given sphere in determining the potential at a given point depends not only upon its distance r' 74 from the point but also on the strength and direction of its polarization. The polarization is in turn dependent upon the stress at the sphere and the piezoelectric tensor d^. A computer program has been written to calculate the electric potential using the transfer function described above for a target composed of many piezoelectric spheres. The program of course requires that the particle displacement function u(r,t) and the piezoelectric tensor cLjk be specified. The calculations are performed in the frequency domain and time series of electric potential are obtained using an inverse FFT. It is therefore desirable to use a function u(r,t) that has an analytic Fourier transform and is reasonably representative of seismic signals generated by real sources (such as explosions). The particle displacement function and piezoelectric tensor chosen for modelling are described below. 3.2 (c) A Model for the Spherical Elastic Wave The particle displacement function chosen for use in piezoelectric simulations is based on a model of the elastic wave produced by an explosive source. The model was developed by J.A. Sharpe (1942). Sharpe suggested that the explosion could be approximated by the application of a pressure pulse to the interior surface of a spherical cavity in an ideally elastic medium. The model is illustrated schematically in Figure 3-3. The radius of the cavity is a and the pressure pulse is denoted as q(t). The medium surrounding the cavity is assumed to be homogeneous and infinite (i.e., a whole space) so that no reflections or refractions are generated. Sharpe also simplified the model by assuming the elastic constants X and \i were equal. He claimed that the assumption was not required but resulted in a considerable simplification of the analysis and was "substantially correct" for material below the low velocity layer. Since the pressure is applied uniformly 75 over the interior surface of the cavity, the particle displacement is purely compressional and in the radial direction. Starting with the elastic equation of motion and appropriate boundary conditions, Sharpe derived the particle displacement u in the medium resulting from the application of a step change in the pressure to the interior of the cavity. For a step change in pressure of the form JO t<0 q(t) = H(t) q 0 where H(t) = U t > Q the resulting particle displacement is u,(r,X) = H(X) aqo 4*1 ( f - V I -co, rjtA/2 sin(coox + arctan\/2) + >/2 a] -co0TA/2 . Lrj (3.7) where T = t — r-a I v j 2y[2v C°o = - 3 a a + 2ti /3u. . , . , and v =\ I — = "W ^ is the compressional wave velocity. The particle displacement given by equation (3.7) is made up of three terms, two of which are oscillatory. The first two terms depend on the distance from the source as (a/r)2 while the third is proportional to a/r. As a result, at distances of more than a few times the cavity radius, the particle displacement may be approximated by the exponentially damped sinusoid a 2 Oo -co„xA/2 us(r,T) = H ( T ) ^ - - - e sinco0T. (3.8) 76 The frequency to0 of the sinusoid is directly proportional to the p-wave velocity and inversely proportional to the cavity radius. The damping is high. For example, if the frequency of the sinusoid is 500 Hz, the amplitude of the signal has dropped to lie of its peak value within >/2/(27t500) ~ 0.5 milliseconds or one quarter period of the sinusoid. The particle displacement resulting from the application of an arbitrary pressure function q(t) can be calculated by convolving q(t) with the particle displacement impulse response. The impulse response is obtained from the response to a unit step in pressure (i.e., qQ=l) by differentiation as follows 9 uI(r,T)=— us(r,T) . Sharpe (1942) used this approach to compute the particle displacements generated by several pressure functions. Alternatively, the particle displacement can be calculated in the frequency domain as the product of the Fourier transforms of the pressure function and the impulse response. This approach is more useful for piezoelectric modelling since the calculations required to calculate electric potential are most easily performed in the frequency domain (using equation 3.4). At this point it is convenient to re-write the approximate expression for the unit step response in terms of time t as follows a2 -co0(t-t0)/>/2 us(r,t) = H(t-t0) e ° ' V sinco0(t-t0) (3.9) r-a where t0 = — and the value 1 has been substituted for q0. 77 The particle displacement produced by an arbitrary pressure function q(t) is then given by d u(r,t) = q(t) <8> ^ us(r,t) where <8> denotes convolution. Taking the Fourier transform of both sides of the above equation we obtain U(r,co) = 7 | q ( t ) ® ^ u s ( r , t ) = {^q(t)} x ^ u 3 ( r , t ) U(r,to)=;coQ(co)Us(r,co) The Fourier transform of the approximate unit step response given by equation (3.9) is derived in Appendix B. The resulting expression is Us(r,co) = _a* -;'cot0 2>/2p: <f2 Therefore, the Fourier transform of the particle displacement generated by an arbitrary pressure function q(t) is U(r,co) = 2>/2p: -Mo ;coco0Q(co) j(0 + V"2 (3.10) for r » a where Q(co) = J{q(t)}. 78 The components of the stress tensor for a spherical seismic wave depend upon U(r,co) and 3U/9XJ as shown by equation (3.6). The spatial derivative, calculated keeping in mind that r and tQ are both dependent on Xj since r = -^Xj 2 + x 2 2 + x 3 2 and tQ = (r-a)/v, is given by dU(r,co) Xj - 3 — = - U ( r , c o ) 7 I v r. (3.11) Substituting equations (3.10) and (3.11) into (3.6), and applying the assumption X= p we obtain (after some algebra) the following expressions for the stress tensor aii(x1,x2,x3,Q)) = pU(r,co) 3 f_x£ ja r r 3 v 2 x :2 + 1 -2pXixkU(r,co) fjco 2I . , o^x^x^.to) = ~2 j — + ~ | , l * k (3.12) where r = •\Jx^+ x 2 2 + x 3 2 and U(r,ca) is given by equation (3.10). Equations (3.12) are based on the approximate step response given by equation (3.9) and are therefore good approximations to the true stress at distances r of more than a few times the cavity radius a. It is interesting to note that as r -»<*>, the maximum normal stress ci{, which occurs for r=x;, is given by (aii)max = -M-U(r,coy 2X: 2 1 7 7 + 1 ;coU(r,co) 79 but v =A/ — = A / p (since we have assumed A, = p.) (<*«)«« = -M/Sip U(r,co). Recall that, for plane waves, acoustic impedance Z = pv = "\J3\ip for X, = p.. (a a)M X = -ya)ZU(r,<o) or, in the time domain 3u (° i i )max = - z 9 7 • Thus, as r-»°°, and the spherical wave becomes more like a plane wave, the maximum normal stress is given by the product of acoustic impedance and particle velocity. This approximation was used by Russell and Barker (1991) to estimate stress. In the simulations presented in this thesis however, the full stress tensor for a spherical elastic wave is calculated using equations (3.12). Figure 3-3 (adapted from Sharpe, 1942). Illustration of the model used by Sharpe (1942) to estimate the elastic wave particle displacement generated by an explosive source. The explosion is approximated by the application of a pressure pulse q(t) to the interior of a spherical cavity in an ideally elastic medium. Figure 3-4 (from Bishop, 1981). Symmetry of the crystallographic axes within a quartz aggregate having the piezoelectric fabric °° m. The aggregate contains equal proportions of left and right-handed quartz. One a-axis from each crystal is polar parallel to an a-axis from every other crystal. The c-axes are randomly oriented in a plane. A rectangular coordinate system is shown with its x 3 axis parallel to the infinity axis of the fabric. 81 3.2 fd) A Model for the Piezoelectric Fabric The concept of a piezoelectric fabric was introduced in section 1.3 of the Introduction. For purposes of piezoelectric modelling, the fabric of a target is important because it determines the form of the piezoelectric tensor d^. It is also very difficult to measure partly because the fabric of a natural piezoelectric aggregate is expected to be somewhat heterogeneous (Bishop, 1981; Ghomshei and Templeton, 1989). There are six ideal piezoelectric fabrics for quartz aggregates (Bishop, 1981). The fabrics are ideal in the sense that the piezoelectric axes of all crystals in the aggregate must conform to one of six symmetry groups. Only two of the ideal fabrics are possible for aggregates containing equal proportions of left- and right-handed quartz. Natural quartz aggregates are usually composed of equal volumes of the two hands although aggregates with a significant excess of one hand have also been reported (Bishop, 1981; Ghomshei and Templeton, 1989). A natural quartz aggregate may or may not possess a piezoelectric fabric. Furthermore, if present, the fabric will not in general correspond exactly to any of the six ideal fabrics. For purposes of modelling, however, we will assume that the target has the ideal fabric °°m. This is one of the two ideal fabrics that can exist in aggregates having equal proportions of left- and right-handed quartz. Furthermore, the results of a laboratory study of samples from the Fairview, B.C. quartz vein (Ghomshei and Templeton, 1989) suggested that the fabric of the vein was similar to °°-m. The alignment of piezoelectric axes for the «»-m symmetry is depicted in Figure 3-4. One a-axis from each crystal in the quartz aggregate is aligned polar parallel to a common direction (the x3 direction of the coordinate system shown in Figure 3-4). The other a-axes 82 are oriented over a cone and the c-axes are randomly oriented in a plane. Note that the polar parallel a-axes lie along a symmetry axis of infinite fold. That is, the piezoelectric properties of the aggregate remain unchanged for any rotation about this axis. With the infinity axis oriented parallel to the x3 coordinate axis, the piezoelectric moduli cLj of the °°-m fabric are (dij) 0 0 0 0 -a„ 0 0 0 0 A 0 0 2 -an 2 an 0 0 0 (3.13) where 3H= -2.3 x l O 1 2 C/N is the longitudinal piezoelectric modulus for a right-handed quartz crystal (Bishop, 1981). After the type of fabric to be used for modelling has been established, there remains the question of how the fabric should be oriented within the target. Ghomshei and Templeton (1989) found that the fabric of the Fairview, B.C. vein could be approximated by an °°-m fabric for which the <*> axis was roughly vertical and within the plane of the steeply dipping vein. The effect of fabric orientation on the electric potential generated by a vein-like target is investigated in the simulations that follow. 83 3.3 Simulations The steps involved in generating synthetic records of the electric potential due to a target composed of many piezoelectric spheres may be summarized as follows: (i) Calculate the seismic stress at each sphere in the target. (ii) Multiply the stress tensor by the piezoelectric tensor of the sphere to calculate its polarization. (iii) Calculate the potential generated by each polarized sphere. (iv) Calculate the sum of the potentials due to all spheres in order to approximate the potential due to the target as a whole. The models and equations used to implement these four steps have been described in section 3.2. A brief summary of the relevant equations is provided below. Using Sharpe's (1942) model for a spherical compressional elastic wave produced by explosive pressures, the stress tensor is given by equations (3.12) aii(x1,x2,x3,co) = pU(r,co) + 1 (3.12) where the particle displacement U(r,co) is defined by equation (3.10) 84 Note that equations (3.10) and (3.12) are far-field approximations to the stress tensor and particle displacement valid for r » a . Once the type and orientation of the piezoelectric fabric has been chosen, the polarization of each sphere is calculated using equation (1.1) Pi(x1,x2,x3,Q)) = dijko^(x1,x2,x3,a)) (1.1) or, in matrix notation Pi(xi,x2,X3,co) = dijoj(x1,x2,x3,co) (1.3) Note that targets having heterogeneous piezoelectric fabrics could be modelled by making clj a function of position. This would be implemented by assigning a different cLj to each sphere. For the simulations presented below, however, the d^  is assumed to be identical for all spheres in the target. The electric potential generated by each polarized sphere is calculated using the model of Russell and Barker (1991). The total electric potential at a point is obtained using equation (3.4) ^ ( x „ x2 , x3, co) -j _ [ yto l X ^ V / P ? ^ - ^ 47ceo [/co + 2/(3pe0)J / j ( r ^ /=l (3.4) where r^ is the displacement from sphere / to the point (xj,x2,x3) at which potential is being calculated. 85 Note that all calculations are performed in the frequency domain. Time domain records of the electric potential \|/(xi,x2,x3,t) are obtained by applying an inverse FFT to Y(X 1,X 2,X 3,CD). The FORTRAN source code, called PZM.FOR written to carry out the calculations, is listed in Appendix A. The program was run on the MTS mainframe computer system at UBC One part of the model that has been ignored up to this point is the form of the pressure function Q(0)) (q(t) in the time domain) needed to specify the particle displacement in equation 3.10. One function that provides a plausible approximation to the pressure generated by an explosions and can be expressed analytically in both the time and frequency domains is q(t) = H(t)bt2e (3.14a) Q(oo) = 2b 1 a + (3.14b) where a and b are constants given by cc-2G>0=4V2v/3a and b = qmax(ae)2/4 . The pressure q(t) has the form of a pulse beginning at time t=0 and reaching a maximum q ^ at time t= 2/a. This pressure function is used for all of the simulations presented below. Figure 3-5 illustrates the geometry used for modelling. The shotpoint was located at the origin of a rectangular coordinate system Ox,x2x3 eight metres away from the piezoelectric target. One spherical and two planar targets were considered. The planar Figure 3-5. Geometry used for piezoelectric modelling. The shotpoint was located at the origin of a rectangular coordinate system and targets were centred at point (8m ,0,0). Planar targets, representing quartz veins, were located in the plane xl = 8m. The potential was calculated at 70 points located along the dotted lines and at the intersections of the grid shown in the x 1x2 plane. 87 bodies - intended to represent quartz veins - were located in the plane x = 8m. If the x3 direction is considered 'up', then the planar bodies were vertically oriented. The electric potential was calculated at 70 points (potential points) at the intersections of the grid shown in the horizontal (x,x2) plane and along the dotted horizontal and vertical lines. The distances between the target and the potential points varied from 5 to 100 m. The values assigned to the constants needed to specify the particle displacement u are listed in Table 3-1. The maximum pressure q ^ was chosen in order to obtain a peak particle velocity x> of approximately 3.75 x 103 m/s at a distance of 10 m from the shotpoint. Russell and Baker (1991) determined that a particle velocity of this magnitude could reasonably be generated by a 3 kg forcite explosion or a 50 kg weight dropped from a height of 2 m. The estimate is also comparable to particle velocities measured by Sharpe (1942b) during a study of elastic waves produced by explosions; he measured velocities ranging from 2 x 10-3 to 1 x 10"5 m/s using calibrated geophones in boreholes at distances of 5 to 100 m from 1/8 pound of 60 percent gelatine dynamite. The value of 2 x 1010 N/m2 assigned to the elastic constants X and p. (recall that it is assumed X = u.) is characteristic of the shear modulus reported by Birch (1966) for samples of granite. The value implied for Poisson's ratio (o = X,/2(A.+u.)) is 0.25 which lies midway between the values of 0.05 and 0.45, characteristic of very hard, rigid rocks and soft, poorly consolidated materials respectively (Telford et al., 1976). The radius a of the spherical cavity was chosen to obtain a reasonable dominant frequency for the elastic wave. Recall that the dominant frequency of the particle displacement step response decreases as a is increased according to the relationship coo = 2y[2 v/3a. Also note that a is more than five times smaller than the maximum distance of eight metres between the shotpoint and target. Thus, the far-field approximation that has been used to calculate the particle displacement at the target is reasonable. 88 Table 3-1. Constants used to specify the spherical elastic wave according to the model of Sharpe (1942). Description Symbol Value maximum pressure Imax 250x106 N elastic constants (k = p) n 2.0xlO , 0 N /m 2 radius of cavity a 1.5 m seismic velocity V 3000 m/s Figure 3-6 illustrates the pressure function q(f) and the resulting particle displacement u (r,t) and particle velocity u (r,t) = du/dt at a distance of eight metres from the shotpoint. The particle displacement and velocity were obtained in the time domain by taking the inverse FFT of U (r,co) and jcoU (r,co) respectively. The explosive pressure is a single pulse which peaks at 0.5 ms and has a total duration of approximately 2 ms. It gives rise to a quasi-sinusoidal particle displacement. The dominant frequencies of u and v are approximately 300 and 400 Hz respectively. The difference in dominant frequency is to be expected since velocity is related to displacement by differentiation in the time domain, or by multiplication by jco in the frequency domain. The particle displacement and velocity signals arrive at time t=2.2 ms. This corresponds to the time required for the elastic wave to travel the 6.5 m from the surface of the spherical cavity (at r=a=1.5m) to r=8m at a velocity of 3000 m/s. The constants needed to specify the electrical properties of the model include the piezoelectric moduli dy of the target and the resistivity of the whole space surrounding it. All simulations were carried out using a resistivity of 1000 £1 m - a value that is typical of 89 EXPLOSIVE P R E S S U R E q( t ) 2 . 5 0 E + 0 8 -O.OOE+OO -5 10 15 Time" ( m s ) PARTICLE DISPLACEMENT (at r = 8 m ) 5 10 15 T ime ( m s ) PARTICLE VELOCITY (at r = 8 m ) 2 0 0 . 0 0 4 0 . 0 0 2 -E 0 . 0 0 0 - 0 . 0 0 2 -- 0 . 0 0 4 10 T ime ( m s ) Figure 3-6. Generation of a spherical elastic wave using the model of explosive pressures developed by Sharpe (1942). The explosion was approximated by the application of a pressure pulse (top) to the interior of a spherical cavity in an ideally elastic medium. The particle displacement and particle velocity of the resulting compressional wave are shown at a distance of eight metres from the shotpoint. 90 the host rock resistivity at the Fairview, B.C. quartz vein site (Russell and Barker, 1991). The ideal piezoelectric fabric °°- m was assigned to each of the three simulated targets but with two different orientations. The shape, dimensions, and piezoelectric moduli of the three targets are described below with reference to the rectangular coordinate system shown in Figure 3-5. Target 1 consists of a single piezoelectric sphere one metre in diameter centred at point (8m, 0, 0). The °° axis of the piezoelectric fabric is oriented in the x2 direction. Note that the moduli given for the « • m fabric by (3.13) apply for the case in which the °° axis is parallel to x3. By expressing (3.13) in its full tensor form (dijlc) and applying standard tensor transformation principles, it can be shown that the appropriate moduli for Model 1 are 0 0 0 0 0 A 2 2 3n 0 0 0 0 0 0 -an 0 0 (3.15) where dn = -2.3 x l O 1 2 C/N is the longitudinal piezoelectric modules for a single right-handed quartz crystal. Target 2 consists of a single layer of piezoelectric spheres lying in the plane x, = 8, with its centre at point (8m, 0, 0) and lateral dimensions of 44 x 44 m. The target's orientation is shown in Figure 3-5 (with reduced dimensions). A total of 1936 spheres, each one metre in diameter, were stacked side by side to form the layer. The °° axis of the piezoelectric fabric lies in the plane of the layer and is oriented in the x2 direction. The piezoelectric moduli are therefore given by (3.15). This target is intended to approximate a quartz vein. Its piezoelectric fabric approximates, in both form and orientation, the fabric 91 observed by Ghomshei and Templeton (1989) in samples from the Fairview, B.C. quartz vein. Target 3 is identical to Target 2 except that the °° axis of the piezoelectric fabric is oriented in the Xj direction perpendicular to the plane of the simulated vein. The form of the piezoelectric matrix (djj) for this fabric orientation is -an 0 0 0 0 0 0 0 0 A (3.16) . 0 0 0 0 - 3 M 0 J Note that the targets each have a thickness of just one sphere. Since individual spheres are modelled as uniformly polarized bodies, the effects of the target thickness have not been modelled using these targets. However, it is expected that the effects of non-uniform polarization in thick piezoelectric bodies will be similar to the effects of non-uniform polarization in laterally extensive bodies such as Targets 2 and 3. Figure 3-7 illustrates the relationship between particle velocity and the electric potential generated at a point by the three targets. The particle velocity x> (r,f) was calculated for r = 8m which is the minimum distance between the shotpoint and each target. The electric potential was calculated for the point with coordinates (18, -10,0) (in metres) located ten metres away from the vein-like target as shown in Figure 3-5. There is a significant difference between the potential generated by a single sphere (Target 1) and that generated by the simulated veins (Targets 2 and 3). The potential due to a single sphere is approximately proportional to the first time derivative of particle velocity 8u/9t. As discussed in section 3.2(a), Russell and Barker (1991) predicted that such a relationship would apply for the case of a plane seismic wave. In contrast, the potentials generated by the PARTICLE VELOCITY and ELECTRIC POTENTIALS PARTICLE VELOCITY AT r=8 m TARGET 1 POTENTIAL TARGET 2 POTENTIAL TARGET 3 POTENTIAL I I i I i i i I i i ;—I 0 2 4 6 8 10 12 14 16 18 20 TIME (ms) Figure 3-7. Comparison of the particle velocity at a distance of eight metres from the shotpoint to synthetic records of electric potential calculated for the three piezoelectric targets. The potential at point (18m, -10m, 0) is shown. This point is located on the grid shown in Figure 3-5. The records generated by the two planar targets (Targets 2 and 3) resemble the particle velocity. The record generated by a single sphere (Target 1) appears to resemble the first derivative of particle velocity. two vein-like targets bear a striking resemblance to the particle velocity itself. Thus, for Targets 2 and 3, the transfer function relating particle velocity and electric potential is well approximated by xP(x1,x2,x3,co) = k!F{x>(T,t)} or, in the time domain, by \|/(x1,x2,x3,t) = k D(r,t) where k is a constant. The value of this constant is not important if the objective is simply to estimate the shape of the electric potential signal. In Figure 3-7 it is clear that the signals generated by Targets 2 and 3 are of lower frequency than the signal due to the single sphere. This reflects the fact that the spheres making up the vein-like targets become polarized at different times and do not contribute in perfect phase to the electric potential. The potential due to Target 2 displays a trailing pulse which arrives approximately 5 ms after the onset of the initial signal. The late pulse originates from spheres located near the edges of the square, planar target. There is little evidence of this trailing pulse in the synthetic potential for Target 3. This is attributed to the difference in fabric orientation. Near the outer limits of Targets 2 and 3 the components of normal stress would have been relatively strong in the x2 and x3 directions. The piezoelectric moduli for Target 2 (given by (3.15)) are more strongly coupled to normal stress in the x2 direction than are the moduli for Target 3 (given by (3.16)). The presence of the trailing pulse in the potential due to Target 2 suggests that the piezoelectric method may be capable of distinguishing the edges of targets having limited dimensions. More generally, pulses may be generated wherever a large scale heterogeneity in the piezoelectric fabric is encountered. Note that, with the exception of the trailing pulse, the potential due to Targets 2 and 3 decays to zero after one and a half quasi-sinusoidal cycles. It appears that the homogeneous piezoelectric fabrics gave rise to destructive 94 interference patterns that effectively cancelled the contributions to electric potential of spheres located beyond a certain distance from the shotpoint. At the edges of Target 2, where the homogeneous fabric was interrupted, the efficiency of the destructive interference waned, thereby giving rise to the trailing pulse. Figure 3-8 shows the synthetic records of electric potential calculated for the three targets at eight points along the line Xj = 18 in the XjX2 plane. The line is parallel to and ten metres away from the vein-like targets as shown in Figure 3-5. Note that the three sets of records exhibit significant differences with respect to amplitude, symmetry, and waveform. The waveforms of signals produced by the three targets have been compared above. However, it is of interest to note how the trailing pulse in the records for Target 2 grows as the edge of the target is approached. The pulse is largest at point (18, 20, 0) which is the closest to the target edge (recall that the target is 44 x 44 m and centred at point (8,0,0)). The potentials generated by Targets 1 and 2 have the same symmetry with respect to position but differ in magnitude by factors of 30 - 40. The difference in magnitude is relatively small considering that the volume of Target 2 is 1936 times larger. For a uniformly polarized target, one would expect the potential to be proportional to target volume as demonstrated by Russell and Barker (1991). The similar symmetries exhibited by the synthetic potentials for Targets 1 and 2 are attributed to the fact that they have the same piezoelectric moduli djj. The differences in the magnitudes and symmetries of synthetic records for the two vein-like targets are due to differences in the piezoelectric moduli of the two targets. Figure 3-8 shows that the potentials due to Targets 2 and 3 are symmetric and anti-symmetric respectively along line Xj = 18, about point (18,0, 0). Thus, the symmetric deployment of a 95 dipole sensor, about point (18, 0, 0) would yield zero potential difference for Target 3 but a significant difference for Target 2. Figure 3-9 shows how the electric potential \|f(xj,x2,X3,t) at a fixed time decays as the distance to the target is increased along line x2 = 2m in the x tx2 plane (see Figure 3-5). The electric potential at each point along the line was divided by the potential at the point closest to the target to obtain the normalized values used in the plot. The fixed times used for Targets 1, 2, and 3 were 2.27,2.50, and 2.50 ms respectively. These are the times for which the electric potential was maximum. The solid and dotted lines in Figure 3-9 have slopes corresponding to 1/r2 and 1/r3 respectively. Note that the potentials due to the two vein-like targets decay at significantly different rates. This must be attributed to the differing orientations of their piezoelectric fabrics since the targets were otherwise identical. Targets 1 (the sphere) and 2 have the same fabric orientation and their potentials decay at comparable rates. Intuitively, one might not have expected this result since Target 1 is essentially a point source while Target 2 appears to be a broad planar source. However, Target 2 was non-uniformly polarized, and at any given time, only part of it contributed constructively to the external potential. ELECTRIC POTENTIAL A L O N G LINE X j = 18 m IN THE X X 2 PLANE TARGET 1 (SINGLE SPHERE) 0.2 mV 5 10 15 T i m e ( m s ) 20 TARGET 2 10 mV 10 Time (ms) 15 20 TARGET 3 40 mV 10 Time (ms) 15 20 (18,-5,0) (18,-2,0) (18,-1,0) (18,0,0) (18,1,0) (18,2,0) (18,5,0) (18,10,0) (18,15,0) (18,20,0) (18,-5,0) (18,-2,0) (18,-1,0) (18,0,0) (18,1.0) (18,2,0) (18,5,0) (18,10,0) (18,15,0) (18,20,0) (18,-5,0) (18,-2,0) (18,-1,0) (18,0,0) (18,1,0) (18,2,0) (18,5,0) (18,10,0) (18,15,0) (18,20,0) Figure 3-8. Synthetic records of the electric potential generated at ten points along line xl=18m by three piezoelectric targets. The line lies in the xlx2 plane at a distance of ten metres from the planar vein-like targets (Targets 2 and 3). The point coordinates (in metres) relative to the coordinate system shown in Figure 3-5 are listed to the right of each trace. NORMALIZED POTENTIAL vs. DISTANCE 0.1 0.001 0.0001 i 8 8c 1 V v— V \ > * • Target 1 D Taraet 2 — > j] • • Tar get 3 s1 1 s s s—-r 1 • 1/(r V - | < • -s ] • - - " - 1/(rA3) N 1 1 } 1 — v - s—l, w 1 ,J—' i N • > S 1 10 100 DISTANCE FROM PLANE Xl=8m(m) Figure 3-9. Decay of normalized electric potential at a fixed time as the distance from plane Xj=8 m (the plane containing the piezoelectric targets) is increased along line x2=2 m. The electric potential has been normalized for each target by dividing the potential at each point by the potential at the point located closest to plane Xj=8 m. 98 4 C O N C L U S I O N S The overall objective of this thesis has been to further the development of the Piezoelectric Exploration Method. This objective has been approached in two ways: (1) Data processing techniques have been developed to enhance the signal-to-noise ratio of data acquired during piezoelectric field experiments. (2) An analytic modelling technique has been developed to investigate the nature of piezoelectric signals from natural targets. As a member of the UBC seismoelectric research group, the author has been involved in all stages - design, acquisition, processing, and interpretation - of several field experiments. Significant problems and questions that could be addressed by data processing and modelling were identified on the basis of this experience. The many sources of noise that can contaminate or completely dominate piezoelectric records are described in Chapter 2. Three processing techniques - including a stacking algorithm, sinusoid subtraction, and low-pass filtering - have been developed to improve the signal-to-noise ratio. The algorithms have been implemented as FORTRAN computer programs that can be run on the DOS-based personal computers used for data acquisition. The stacking technique employs a robust triggering algorithm that enables it to be used effectively even when the trigger signal is very poor. The development of the algorithm allowed stacking experiments to be carried out with the trigger geophone located several metres away from the seismic source. In this way it was demonstrated that suspected piezoelectric signals generated at the time of the source impact could not be attributed to cross-talk or radiation from the trigger geophone. 99 The frequency domain low-pass filter is a simple algorithm that calculates the FFT of a time series, multiplies it by transfer function having a variable cut-off frequency and slope, and transforms the product back to the time domain. Unlike analog filters applied during data acquisition, the frequency domain filter does not influence the phase at any frequency. More importantiy however, the availability of a simple low-pass program allows piezoelectric records to be acquired with a wide bandwidth and low-pass filtered later if necessary. Sinusoid subtraction has proven to be the most useful of the three processing techniques. Its purpose is to subtract powerline harmonics from piezoelectric records. During most of the piezoelectric experiments conducted by the UBC group, powerline harmonics have been the highest amplitude noise. The ambient electric field in the earth at 60 Hz is commonly one to three orders of magnitude larger than the amplitudes expected of piezoelectric signals. Consequently, dramatic improvements in the signal-to-noise ratio have been obtained by the use of the sinusoid subtraction technique. The amplitude and phase of each harmonic frequency is estimated by calculating its Fourier series coefficients. The presence of piezoelectric signals or other types of noise in the data influences these coefficients and limits the accuracy with which the powerline harmonics can be estimated. Experience has shown that it is sometimes preferable to base the estimate on a portion of the data recorded before the piezoelectric arrival rather than on the whole data record. The sinusoid subtraction program is now used routinely during experiments conducted by the UBC seismoelectric research group. The incentive for analytic modelling was provided by an interest in establishing the transfer function relating seismic and piezoelectric signals. Such a transfer could be applied to geophone signals recorded during seismoelectric experiments in order to estimate the form of any piezoelectric signal that might be generated. Ideally, the estimated signal could then 100 be used to deconvolve piezoelectric records. Russell and Barker's (1991) model of a uniformly polarized sphere in a conductive whole space was used as the starting point for the modelling presented in this thesis. A target having non-uniform polarization was approximated by a group of many piezoelectric spheres each of which became independendy polarized during the passage of a spherical elastic wave. Synthetic time series of the electric potential \|/(x1,x2,x3,t) generated by three targets located eight metres away from a shotpoint were computed using Sharpe's (1942) model to simulate the spherical elastic wave. The targets included a single uniformly polarized sphere, and two broad planar bodies intended to represent quartz veins. All three targets had the same piezoelectric fabric but the fabric orientation differed for one of the planar targets. The three sets of synthetic records displayed differences related to both non-uniform polarization and fabric orientation. A comparison of the elastic wave particle velocity to the electric potentials at a point showed that the relationship between seismic and piezoelectric signals can be significantly influenced by non-uniform polarization. The potential generated by single piezoelectric sphere was approximately proportional to the first derivative of particle velocity. In contrast, the potentials generated by the two non-uniformly polarized vein-like targets were lower frequency and resembled the particle velocity itself. A second effect associated with non-uniform polarization is the generation of 'trailing pulses'. After a short time (one and a half quasi-sinusoidal cycles), the electric potentials generated by the planar targets decayed to zero due to destructive interference caused by non-uniform polarization. However, as the elastic wave passed the edges of one planar target, a second pulse of electric potential was generated. If 'trailing pulses' of this kind could be identified in real piezoelectric data, they could yield important information on target dimensions. 101 The rate of decay of electric potential with distance from the piezoelectric targets was also investigated. The potentials due to the single, uniformly polarized sphere and the broad, planar target having the same piezoelectric tensor decayed at comparable rates. The similarity of the decay rates for a point source (the sphere) and a planar source were attributed to the non-uniform polarization of the planar target. The effects of piezoelectric fabric orientation were evaluated by comparing the electric potentials generated by two vein-like targets that were identical in all respects except for their fabric orientations. The electric potentials due to the two targets had different magnitudes, different symmetries with respect to position, and exhibited different rates of decay with distance from the target. 102 REFERENCES Bishop, J.R. 1978. An investigation into piezoelectric effects in rocks, unpublished Ph.D. thesis, Macquarie University. Bishop, J.R. 1981. Piezoelectric effects in quartz-rich rocks. Tectonophysics 11, 297-321. Blau, L.W. and Statham, L. 1936. Method and apparatus for seismic-electric prospecting. U.S. Patent No. 2 054 067, September 15, 1936. Bond, W.L. 1943. A mineral survey for piezoelectric materials. Bell Syst. Tech. J. 22, 145-152. Bracewell, R.N. 1986. The Fourier Transform and its Applications. McGraw-Hill, New York. Butler, K.E. 1990. Two approaches to stacking piezoelectric exploration data. Unpublished term report for Geophysics 514, University of British Columbia. Butler, K.E., Kepic, A.W., Russell, R.D. and Mellema, G.R. 1991. Analysis of piezoelectric signals from quartz veins. Paper presented at the CSEG National Convention in Calgary, Alberta, May 14-16. Cady, W.G. 1946. Piezoelectricity. McGraw-Hill, New York. Dietrich R.V. and Skinner B.J. 1979. Rocks and Rock Minerals. John Wiley and Sons, New York. Ghomshei, M.M. and Templeton, T.L. 1989. Piezoelectric and a-axes fabric along a quartz vein. Physics of the Earth and Planetary Interiors 55, 374-386. Ghomshei, M.M., Narod, B.B., Templeton, T.L., Arnott, A.S. and Russell, R.D. 1988. Piezoelectric pole figure of a vein quartz sample. Textures and Microstructures 1, 303-316. 103 Institute of Radio Engineers, 1945. Standards on piezoelectric crystals. Proc. Inst, of Radio Eng. 49, 1378-1395. Ivanov, I.G. 1939. Effect of electrostatic charging of slabs of earth by passing elastic waves through them. DAN SSSR 24, no. 1. Jose, B.F. 1979. Feasibility of the piezoelectric effect for quartz vein detection, unpublished M.Sc. thesis, University of British Columbia. Kondrashev, S.N. 1980. The Piezoelectric Method of Exploration. Nedra, Moscow (in Russian). An English translation by G.M. Volkoff is available at the University of British Columbia. Maxwell, M. 1988. Catalog of field experiments. In: Seismoelectric exploration review reports. Unpublished internal report of the Geophysical Instrumentation Group, University of British Columbia. Maxwell, M., Russell, R.D., and Narod B.B. 1991. Piezoelectric responses of natural quartz. Manuscript in preparation. Neishtadt, N.M., Mazanova, Z.V., Binevich, L.Ya. and Maiko, M.I. 1972. Piezoelectric method of exploration (methodological recommendations). ONTI VITR, Leningrad. Nye, J.F. 1957. Physical Properties of Crystals. Clarendon Press, Oxford. Parkhomenko, E.I. 1971. Electrification Phenomena in Rocks (translation by G.V. Keller). Plenum, New York. Press, H.P., Flannery, B.P., Teukolsky, S.A. and Vetterling, W.T. 1989. Numerical Recipes (FORTRAN version). Cambridge University Press, New York. Russell, R.D. and Maxwell, M. 1990, Development of a Seismoelectric Exploration Technique, unpublished annual report to NSERC for project no. CRD-0094237, July 31, 1990. 104 Russell, R.D., Maxwell, M., Butler, K.E., Kepic, A.K., and Mellema, G.R. 1991, Development of a Seismoelectric Exploration Technique, unpublished annual report to NSERC for project no. CRD-0094237, July 31,1991. Russell, R.D. and Barker, A.S., Jr. 1991. Seismo-electric exploration: expected signal amplitudes. Geophysical Prospecting^,105-118. Sharpe, J.A. 1942. The production of elastic waves by explosion pressures. 1. Theory and empirical field observations. Geophysics 7,144-154. Sharpe, J.A. 1942b. The production of elastic waves by explosion pressures. 2. Results of observations near an exploding charge. Geophysics 7, 311-321. Sheriff, R.E. and Geldart, L.P. 1982. Exploration Seismology Volumel: History, Theory, and Data Acquisition. Cambridge University Press, Cambridge. Shubnikov, A.V. 1946. Piezoelectric Textures. Akad. Nauk SSSR, Moscow (in Russian). Sobolev, G.A. and Demin, V.M. 1980. Mechanoelectric Phenomena in the Earth. Nauka, Moscow (in Russian). An unpublished partial English translation by G. M. Volkoff is available at the University of British Columbia. Sobolev, G.A., Demin, V.M., Los1 and Maybuk, Yu.Ya. 1982. Study of the electromagnetic radiation of rocks containing semiconductor and piezoelectric minerals. Izvestiya, Academy of Sciences USSR (English translation by A.G.U.) 18, 888-897. Sobolev, G.A., Demin, V.M., Narod, B.B. and Whaite, P. 1984. Tests of piezoelectric and pulsed-radio methods for quartz vein and base-metal sulfides prospecting at Giant Yellowknife Mine, N.W.T., and Sullivan Mine, Kimberley, Canada. Geophysics 49, 2178-2185. Stremler, F.G. 1982. Introduction to Communication Systems. Addison-Wesley, Reading. Telford, W.M., Geldart, L.P., Sheriff, R.E., and Keys, D.A. 1976. Applied Geophysics. Cambridge University Press, New York. 105 Thompson, R.R. 1936. The seismic-electric effect. Geophysics 1, 327-335. Thompson, R.R. 1939. A note on the seismic-electric effect. Geophysics 4,102-105. Tuck, G., Stacey, F.D. and Starkey, J. 1977. A search for the piezoelectric effect in quartz-bearing rocks. Tectonophysics 39, T7-T11. Volarovich, M.P., Parkhomenko, E.I. and Sobolev, G.A. 1959. A field investigation of the piezoelectric effect in quartziferous rock. Proc. (Doklady) Academy of Sciences USSR 128,525. Volarovich, M.P., Sobolev, G.A. and Parkhomenko, E.I. 1962. The piezoelectric effect in pegmatite and quartz veins. Izvestiya, Academy of Sciences USSR, Physics of the Solid Earth (English translation by A.G.U.), p. 145-152. Volarovich, M.P. 1978. Investigations of the physical properties of rocks and minerals at high pressures and temperatures. Izvestiya, Academy of Sciences USSR, Physics of the Solid Earth (English translation by A.G.U.) 14, 696-703. APPENDIX A FORTRAN source code, parameter files, and documentation are listed for the following four computer programs: (1) SINSUB.FOR - used to estimate and subtract harmonics of 60 Hz from piezoelectric records (2) LPF.FOR - a frequency domain low-pass filtering program (3) TRIG43EX.FOR - a stacking program that uses a cross-correlation technique choose trigger points (4) PZM.FOR - a piezoelectric modelling program 107 PROGRAM DOCUMENTATION Program: SINSUB.EXE Function: Sinusoid subtraction (subtraction of powerline harmonics) Written by Karl Butler, May, 1991, Last updated Sept. 10/91 Documentation written May 30,1991, Last updated Oct. 8/91 Source code: SINSUB.FOR Compiler: MS FORTRAN v. 5.0 Description: SINSUB.EXE is used to subtract powerline harmonics (assumed to be harmonics of 60 Hz) from up to eight traces each comprised of up to 4096 points. The magnitude and phase of the sinusoids to be subtracted are found by calculating the trigonometric Fourier Series coefficients for each frequency. The user specifies which harmonics of 60 Hz are to be removed and which channels are NOT to undergo the sinusoid subtraction process. The first operation performed on each channel is the subtraction of its mean value. The harmonic estimation and subtraction process is performed next. Sinusoid subtraction at a particular frequency is carried out on all channels (except those specificied by the user) before proceeding to estimate and remove the next harmonic. After all sinusoids have been subtracted the mean values of the processed traces are removed and all traces are written as REAL*4 numbers to a binary output file. If only one channel is read from the input file the trace is written to an ASCII output file following the subtraction of each harmonic. Note that this program may also be used to subtract the mean values from all channels without subtracting any sinusoids. This is done by specifying in the parameter file that each channel is to be excluded from the sinusoid subtraction process. The neccessary input and output files for program SENSUB.EXE are described below: Parameter file: SINSUB.PRM - contains the name of the binary input data file, the type (INTEGER*4 or REAL*4) of data in the input file, the binary file record length, the number of samples in each trace, the sample interval, the number of traces in the input data, information on which traces are NOT to undergo sinusoid subtraction, the number of harmonics to be subtracted, and a list of the harmonics themselves. See the attached sample parameter file for a more detailed desription. Input data file: The name of this file is specified by the user on the first line of the parameter file. It must be a BINARY file containing multiplexed data of type REAL*4 or INTEGER*4. For example, binary output files generated by the lowpass filter program LPF.EXE and *.DBL files created by stacking program RCED.COM can be used as input data files for this program. Output data file: SINSUB.BIN - a BINARY file containing multiplexed data of type REAL*4 Output text file: SINSUB.TXT - contains information on the parameters used and the amplitudes of the sinusoids subtracted from each trace SINSUB.PRM (October 8,1991 6:32pm) PAY 5 9.DBL (filnam) 32 (RLEN - record length i n input data f i l e ) 1*4 (TYPE - type of input data - 1*4 or R*4) 4096,0.000032 (M - max. of 4096 samples, CELT) 1 ,8 (STCHAN,LACHAN - each may range from 1 to 8) 4,1,2,7,8 (NRW) , (RWCHAN(I) ,1-1, NRW) 25 (NHARM - no. of harmonics may range from 0 to 30) 1 2 3 4 5 (HAR(l)) 6 7 8 9 11 12 13 15 17 19 21 23 25 27 29 31 33 34 35 37 INPUT PARAMETERS: filnami name of the input data f i l e (CHARACTER*40) RLEN: record length In the input data f i l e (INTEGER) TYPE: type of input data (1*4 for INTEGER*4 or R*4 for REAL*4) (CHARACTER*3) M: number of samples i n each trace i n the input data (maximum of 4096 sample S ) (INTEGER) DELT: sample i n t e r v a l (in seconds) of the traces i n the input data f i l e (REAL) STCHAN: f i r s t channel to be read from the input f i l e (INTEGER) LACHAN: l a s t channel to be read from the input f i l e (INTEGER) A l l channel s from STCHAN to LACHAN incl u s i v e are read. Both STCHAN and LACHAN may range between 1 and 8 but LACHAN must be greater than or equal to STCHAN. NRW: The number of channels that are to be read but are not to undergo the harmonic subtraction process. (INTEGER) RWCHAN(I),1-1,NRW : Channel numbers of the NRW channels that are to be read but not undergo harmonic subtraction.(INTEGER) NHARM: The number of d i f f e r e n t harmonics to be subtracted. This number may range from 0 to 30. (INTEGER) HAR(I),I-1,NHARM : the harmonics to be subtracted (INTEGER) o OO sinsub.for (October 8,1991 5:25pm) Q *********************************************************************** C SINSUB.FOR - Written by K. Butler, May, 1991 C Last updated Sept. 10/91 C C This program is used to subtract power line harmonics from up to C eight traces each comprised of up to 4096 points.The magnitude and C phase of the sinusoid to be subtracted at each harmonic frequency C are found by calculating the trigonometric Fourier Series coefficients C for each harmonic. C The user specifies which channels are to be read from the input C data f i l e and specifies which channels are NOT to undergo the C sinusoid subtraction process. The f i r s t operation performed on each C channel is the subtraction of its mean value. Following this, the C harmonic estimation and subtraction process is performed. Sinusoid C subtraction at a particular frequency is carried out on a l l channels C (except those specified by the user) before proceeding to estimate C and subtract the next harmonic. After a l l sinusoid subtractions have C been completed the mean values of the processed traces are removed and C a l l traces (including those from which harmonics were not subtracted) C are written as REAL*4 numbers to a binary output f i l e . If only one C channel is read from the input f i l e (ie. i f NCHAN-1), the trace i s C written to an ASCII ouput f i l e following the subtraction of each harmonic. C C NOTES: C (1) This program is similar to sinusoid subtraction programs C ADPINT.EXE and ADPREAL.EXE but i t has been improved in the C following ways: C (i) Execution time has been reduced by 30-40% by reducing C the number of calls to the trigonometric functions sine and C cosine. The trigonometric functions are evaluated once for each C harmonic frequency and are stored in arrays SI(J) and CO(J). C (ii) The data in the multiplexed binary input f i l e may be of C type INTEGER*4 OR REAL*4. C C (2) The ALLOCATE statement and the allocatable array, option used C for arrays IDATA and RDATA are MS FORTRAN v.5.0 extensions to C the ANSI FORTRAN77 full-language standard. C C C Explanation of INPUT parameters: C C FILNAM: name of the input data f i l e (CHARACTER*40) C RLEN: record length in the input data f i l e (INTEGER) C TYPE: type of input data (1*4 or R*4) (CHARACTER*3) C M: number of samples in each trace in the input data (INTEGER) C DELT: sample interval (in seconds) of the traces in the input C data f i l e (REAL) C STCHAN: fir s t channel to be read from the input f i l e (INTEGER) C LACHAN: last channel to be read from the input f i l e (INTEGER) C Al l channels from STCHAN to LACHAN inclusive are read. C NRW: The number of channels that are to be read but are not to C undergo the harmonic subtraction process. (INTEGER) C RWCHAN(I),1-1,NRW : Channel numbers of the NRW channels that are to C be read but are not to undergo harmonic C subtraction.(INTEGER) C NHARM: The number of different harmonics to be subtracted. (INTEGER) C HAR(I),I-1,NHARM ; the harmonics to be subtracted (INTEGER) C C ************************************************************ INTEGER I,J,K,M,L,NHARM,STCHAN,LACHAN,NCHAN,MP,N,RLEN,RLEN2,NRW INTEGER FLAG,TRACE,RWCHAN(8),INTDUM(7),INTVAR(8),IERR,KK,LIST(8) REAL*8 DATA[ALLOCATABLE, HUGE](:,:),SI[ALLOCATABLE](:) REAL*8 CO [ALLOCATABLE] (:) ,A (30, 8) ,B (30, 8) ,C (30, 8) , PHASE (30,8) REAL* 8 HAR(30) ,FREQ(30) ,AVHAR(30) ,SUM(8) ,MEAN (8), FUND,DELT, PI REAL*8 AVG,ARG1 REAL*4 RDUM(7) ,RVAR(8) CHARACTER*40 FILNAM CHARACTER*3 TYPE ALLOCATE (DATA(8,4096), CO(4096), SI(4096), STAT-IERR) IF (IERR .NE. 0) THEN WRITE(6,*)'Not enough storage for data; execution aborted.' STOP END IF OPEN (UNIT-2, FILE-' SIN SUB. PRM' , STATUS-'OLD' ) READ (2,'(A39)') FILNAM READ(2,*) RLEN READ(2,'(A)') TYPE OPEN (UNIT-1, FILE—FILNAM, STATUS—' OLD', ACCESS-' DIRECT', + FORM—' UNFORMATTED' ,RECL—RLEN) -READ(2,*) M.DEL? READ(2,*) STCHAN,LACHAN READ (2, *) NRW, (RWCHAN (I) , I—1,NRW) NCHAN—LACHAN-STCHAN+1 RLEN2—NCHAN* 4 OPEN(UNIT-3,FILE-'SINSUB.ASC ) OPEN (UNIT-4, FILE-'SINSUB. BIN' , STATUS-'UNKNOWN' , ACCESS-'DIRECT' , + FORM—' UNFORMATTED' , RECL-RLEN2) OPEN (UNIT-7, FILE—' SINSUB. TXT' ) C Read the harmonics to be subtracted. READ(2, *) NHARM DO 10 I—1, NHARM 10 READ (2, *) HAR(I) C Write a message to the screen to remind the user of the type C of data he has specified for the binary input data f i l e . IF (TYPE ,EQ. '1*4') THEN WRITE(6, 890) WRITE(6,891) ELSE IF (TYPE .EQ. 'R*4') THEN WRITE(6,890) WRITE(6, 892) ELSE WRITE(6,893) 893 FORMAK//' DATA TYPE ENTERED INCORRECTLY: EXECUTION ABORTED') STOP END IF END IF 890 FORMAT(//////////' NOTE: THE DATA IN THE BINARY INPUT FILE MUST') 891 FORMAT(' BE MULTIPLEXED AND OF TYPE INTEGER*4.'//) 892 FORMAT (' BE MULTIPLEXED AND OF TYPE REAL*4. ' / /) WRITE(6, 900) FILNAM WRITE (7, 900) FILNAM 900 FORMAT(/' Input data w i l l be read from ',A39) IP (NRW . GT. 0) THEN WRITE(7, 902) STCHAN,LACHAN 902 FORMAT(' Channels',12,' through', 12,' w i l l be read.') 903 FORMAT(' Sinusoids w i l l not be subtracted from channels:',812) WRITE(7,903) (RWCHAN(I),1-1,NRW) WRITE(7,905) NCHAN-NRW 905 FORMAT(' Sinusoids WILL be subtracted from the remaining',12, + ' channels.') ELSE WRITE(7,906) STCHAN,LACHAN 906 FORMAT(' Channels',12,' through', 12,' w i l l be f i l t e r e d . ' ) END IF C The lowest frequency harmonic (ie. the fundamental) has units of C Hz and i s stored i n the variable FUND. FUND-60.0D0 PI-3.1415927D0 C Calculate the number of times (N) that the period of fundamental C sinusoid (0.166666... sec for 60 Hz) divides into the length of the C data trace. Then c a l c u l a t e the number of points MP that w i l l C be used i n the truncated data trace. N-INT (FUND*DBLE (M) *DELT) MP-INT(DBLE(N)/(FUND*DELT) + 0.5D0) C Write some comments to an output f i l e . WRITE(7,925) DELT WRITE(6,926) DELT*1000000.D0 WRITE (7, 927) MP 925 FORMAT(' The sample i n t e r v a l is',F12.7,' seconds.') 926 FORMAT(' The sample i n t e r v a l is',F8.1,' microseconds.'/) 927 FORMAT(' There are',15,' samples i n the truncated traces.'/) WRITE(7,930) WRITE(7,931) 930 FORMAT!/' The Fourier components of the following frequencies') 931 FORMAT(' are to be estimated:') DO 20 I-1,NHARM 20 WRITE(7,*) HAR(I)*FUND WRITE (7,940) WRITE (7,941) 940 FORMAT(' and sinusoids at those frequencies are to be') 941 FORMAT(' subtracted from the input data.') WRITE(7,950) WRITE(7,951) WRITE (7, 960) 950 FORMAT!/' The frequencies nearest to those l i s t e d above whose') 951 FORMAT(' Fourier components can be EXACTLY calculated using the') 960 FORMAT)' given sample i n t e r v a l are:') DO 22 I—1,NHARM 22 . WRITE(7,*) (HAR (I) *DBLE (N) ) / (DBLE (MP) *DELT) C Read i n the input traces. DO 25 I-1,NCHAN 25 SUM(I)-0.0D0 IF (TYPE .EQ. '1*4') THEN DO 28 I-1,M READ(1,REC-I) (INTDUM(J),J-1,STCHAN-1),(INTVAR(K),K-1,NCHAN) DO 27 K-1,NCHAN DATA(K, I)-DBLE(INTVAR(K)) 27 SUM(K)-SUM(K)+DATA(K,I) 28 CONTINUE ELSE DO 30 I-1,M READ(1,REC-I) (RDUM(J),J-1,STCHAN-1),(RVAR(K),K-1,NCHAN) DO 29 K—1,NCHAN DATA (K, I) —DBLE (RVAR (K) ) 29 SUM(K)-SUM(K)+DATA(K, I) 30 CONTINUE END IF WRITE(6,*)'Finished reading data' WRITE(7,922) 922 FORMAT!//' SUMMARY OF SINUSOID SUBTRACTION:') KK-0 DO 100 TRACE—1,NCHAN C Subtract the mean values from a l l traces. MEAN (TRACE)—SUM! TRACE) /DBLE (M) DO 33 I-1,M 33 DATA (TRACE, I) —DATA (TRACE, I) -MEAN (TRACE) IF (NCHAN .EQ. 1) THEN DO 34 I-1,M 34 WRITE(3, 920) DATA(TRACE,I) END IF 920 FORMAT(F14.4) C Determine which traces are to undergo'sinusoid subtraction and C store t h e i r trace numbers i n array LIST(i). FLAG-0 DO 36 J-1,NRW IF (TRACE .EQ. RWCHAN(J)) THEN FLAG-1 END IF 36 CONTINUE IF (FLAG .EQ. 0) THEN KK-KK+1 LIST(KK)-TRACE END IF 100 CONTINUE C For each harmonic calculate the amplitude and phase present i n C each of the traces to be processed. DO 200 I—1,NHARM AVHAR(I)-0.D0 FREQ(I)-HAR (I) *FUND ARG1-2.D0*PI*FREQ(I)*DELT WRITE(6,993) FREQ(I) WRITE(6,994) 993 FORMAT(/' SUBTRACTING',F12.3,' Hz:') 994 FORMAT(' Trace Amplitude (levels)') DO 40 J-l.M SI(J)-DSIN(ARG1*DBLE(J-l) ) CO(J)-DCOS(ARG1*DBLE(J-l)) 40 CONTINUE C Calculate the amplitude and phase of the current harmonic C for each trace In turn. DO 70 K-1,KK L-LIST(K) A(I,L)-0.D0 B(I,L)-0.D0 DO 50 J-l,MP A(I,L)-A(I,L)+DATA(L, J)*CO(J) B(I,L)-B(I,L)+DATA(L, J)*SI(J) 50 CONTINUE A(I,L)-A(I,L)*2.D0/DBLE(MP) :• B(I,L)-B(I,L)*2.D0/DBLE(MP) ' C(I,L)-DSQRT(A(I,L)**2+B(I,L)**2) C Calculate the phase of the harmonic. Since the function C DATAN2 is not defined when both of its arguments are zero, C an IF statement i s used to arbitrarily set the phase to C zero i f both of the arguments are very small. IF (DABS(A(I,L)).LT.l.D-75 .AND. DABS(B(1,1)).LT.1.D-75) THEN PHASE(I,L)-0.D0 WRITE(6,*) 'PHASE ARBITRARILY SET EQUAL TO ZERO FOR TRACE',L ELSE PHASE (I, L) —DATAN2 (-B(I,L) ,A(I,L) ) END IF AVHAR (I)-AVHAR (I) +C (I, L) C Subtract the harmonic from the current trace. DO 60 J-1,M 60 DATA (L, J) - D A T A (L, J) -A (I, L) *CO (J) -B (I, L) *SI (J) WRITE(6,995) L,C(I,L) 995 FORMAT(2X,I2,5X,G14.6) IF (NCHAN .EQ. 1) THEN DO 65 J-1,M 65 WRITE(3,920) DATA(TRACE,J) END IF 70 CONTINUE AVHAR (I)-AVHAR (I) /DBLE(KK) 200 CONTINUE C Subtract the mean values from a l l traces that have been filtered. DO 220 K-1,KK L-LIST(K) AVG-0.D0 DO 205 I-1,M 205 AVG—AVG+DATA(L,I) AVG—AVG/DBLE (M) DO 210 I-1,M 210 DATA(L,I)—DATA(L,I)-AVG 220 CONTINUE C Write out the amplitudes and phases of the harmonics subtracted C - from each trace. DO 250 K-1,KK L-LIST(K) WRITE(7,982) L,STCHAN+L-1 982 FORMAT)/' TRACE',12,'(Input channel',12,'):') WRITE (7,984) WRITE(7,985) 984 FORMAT(' Cosines of the following frequency, phase, and') 985 FORMAT(' amplitude will be subtracted:') WRITE(7,986) 986 FORMAT (4X,' FREQUNCY', 7X,' PHASE' , 6X, 'AMPLITUDE', 7X, ' SIN COMP .' + ,4X,'COS COMP.') DO 230 1-1,NHARM 230 WRITE(7,990) FREQ(I),PRASE(I,L),C(I,L),A(I,L),B(I,L) 250 CONTINUE 990 FORMAT(F12.3,F12.3,4X,G14.6,F12.3,F12.3) C Write out the average amplitudes of the extracted harmonics. WRITE(6,997) WRITE(7,997) WRITE(6,998) WRITE (7,998) 997 FORMAT!/' AVERAGE AMPLITUDES OF THE EXTRACTED HARMONICS:') 998 FORMAT!' Subtraction Harmonic Frequency Amplitude') DO 255 1-1,NHARM WRITE(6,999) I,HAR(I) , FREQ (I) , AVHAR (I) WRITE(7,999) I,HAR(I),FREQ(I),AVHAR(I) 999 FORMAT(3X,I2,10X,F5.1,2X,F12.3,3X,G14.6) 255 CONTINUE C Write the processed data to a binary output f i l e . DO 260 I-1,M 260 WRITE(4,REC-I) (SNGL(DATA(K,I)),K-l,NCHAN) CLOSE (UNIT-1) CLOSE(UNIT—2) CLOSE(UNIT-3) CLOSE(UNIT-4) STOP END 112 PROGRAM DOCUMENTATION Program: LPF.EXE Function: Lowpass filtering Written by K.Butler, May, 1991 Documentation written May 30,1991 Source code: LPF.FOR Compiler: MS FORTAN v.5.0 Description: LPF.EXE applies a zero phase low pass filter to up to eight real-valued time series each comprised of up to 4096 samples. Filtering is performed in the frequency domain by applying a FFT to each trace, multiplying the spectrum by the filter transfer function, and applying an inverse FFT to the product. The transfer function T(f) of the filter is given by T ( f ) = 1 + (f/fc)(s/6) where f is frequency fc is the cutoff frequency (6dB down or half amplitude point) s is the filter slope in dB/octave The cutoff frequency and filter slope are specified by the user. Note that the FFT algorithm used in this program demands that the number of samples in each trace be a POWER OF TWO! The neccessary input and output files for program LPF.EXE are described below: Parameter file: LPF.PRM - contains the name of the binary input data file, the type (INTEGER*4 or REAL*4) of input data, the number of channels in the data, information describing which channels are NOT to be filtered, the the number of samples in each trace, the sample interval, the cut-off frequency (6 dB down point), and the filter slope (in dB/octave). See the attached sample parameter file for a more detailed description. Input data file: The name of this file is specified by the user on the first line of the parameter file. It must be a BINARY file containing multiplexed data of type REAL*4 or INTEGER*4. For example, binary output files generated by program SINSUB.EXE and the *.DBL files created by stacking program RCED.COM can be used as input data files for this program. Output data file: LPF.BIN- a BINARY file containing multiplexed data of type REAL*4 LPF.PRM (October 8,1991 6:23pm) PAY60.BIN (filnam) R*4 (TYPE of input data: 1*4 or R*4) 8 (NCHAN: maximum of 8) 2,7,8 (NRW, (RWCHAN(I) ,1-1,NRW) ) 4096 (NFT: must be a power of 2!, max. of 4096) 0.000032 (DELT) 800.0 (CUTOFF) 24.0 (SLOPE) Explanation of input parameters: filnam: name of binary input data f i l e (CHARACTER*40) TYPE: type of input data (1*4 for INTEGER*4 or R*4 for REAL*4) (CHARACTER*3) NCHAN: number of channels i n the input data f i l e (maximum of 8), (INTEGER) NRW: number of channels that are not to be f i l t e r e d , (INTEGER) RWCHAN(I): channel numbers of those channels that are not to be f i l t e r e d (INTEGER NFT: The no. of samples i n each trace This number MUST BE A POWER OF TWO and less than or equal to 4096. (INTEGER) DELT: sample Interval of data i n seconds, (REAL) CUTOFF: cu t o f f frequency (6dB down point) i n Hz, (REAL) SLOPE: f i l t e r slope above the cutoff frequency i n db/octave, (REAL) LPF.FOR (May 30,1991 4:27pm)  c #*****#*****•********#**+********************^  C Program LPF.FOR - written by K.Butler, Hay, 1991 C C Description: This program applies a rero phase low pass filter to up C to eight real-valued time series each comprised of up to 4096 samples. C The filter is applied in the frequency domain and all calculations are C performed using double precision (REAL*8) variables. The transfer C function T(f) of the filter is given by C C T(f) - l/(l+(f/fc)*Ms/6)) C where f - frequency C fc - cutoff frequency <6dB down point) C s - filter slope in dB/octave C Subroutine REALFT (slightly modified for use with REAL*8 data) from the C NUMERICAL RECIPES library is used to perform the neccessary forward and C inverse FFTs for each trace. NOTE that the number of samples (NFT) in C each time series MUST BE A POWER OF 2 (a restriction imposed by C subroutine REALFT). C Data in the binary input file may be of type REAL*4 or INTEGER**. Up to C eight channels may be present in the input file and it is assumed that C the data is stored in a multiplexed format (ie. each record is assumed C to contain one sample from each channel). The binary output file LPF.BIN C contains data of type REAL*4 and is written in the same (multiplexed) C format as the input file. C The following parameters must be specified by the user in the input C parameter file LPF.PRM: C (1) FILNAM: name of the binary input file (CHARACTER*40) C (2) TYPE: type of input data (1*4 or R*4) (CHARACTERS) C (3) NCHAN: number of channels to be read and filtered, (INTEGER) C (4) NRW: number of channels that are not to be filtered (INTEGER) C (5) RWCHAN: channel no.s of those channels that are not to be filtered (INT) C (6) NFT: no. of samples in each trace (MUST be a power of 2), (INTEGER) C (7) DELT: sample interval of the data in seconds, (REAL) C (8) CUTOFF: cutoff frequency (6dB down point) in Hi, (REAL) C (9) SLOPE: filter slope above the cutoff frequency in dB/octave, (REAL) C C Note: ALLOCATABLE arrays are used to reduce the size of the C executable file. This is an extension to the FORTRAN77 standard C found in the MS FORTRAN v.5.0 compiler. C C a************************************************** INTEGER I, J,RLEN, NFT,NCHAN, IERR,NRW, RWCHAN (8), INTVAR(8) REAL*8 DELT,CUTOFF,SLOPE REAL*8 Y[ALLOCATABLE, HUGE](:,:), FILT[ALLOCATABLE](:) REAL*4 RVAR(8) CHARACTER*40 FILNAM CHARACTER*3 TYPE ALLOCATE (Y(4096,8), FILT(4096), STAT-IERR) IF (IERR .NE. 0) THEN WRITE(6,*)'Not enough memory for data; execution aborted.' STOP END IF OPEN (UNIT-2, FILE-' LPF.PRM' .STATUS-' OLD' ) READ (2,'(A39)') FILNAM READ (2, ' (A) ' ) TYPE READ(2,*) NCHAN READ (2, *) NRW, (RWCHAN(I), I—1,NRW) RLEN—NCHAN* 4 OPEN (UNIT-1, FILE-FILNAM, STATUS—'OLD' , ACCESS-'DIRECT' , + FORM—'UNFORMATTED', RECL-RLEN) OPEN (UNIT-3, FILE-' LPF.BIN', STATUS-'UNKNOWN', ACCESS-'DIRECT', + FORM—'UNFORMATTED' , RECL-RLEN) OPEN (UNIT-7,FILE-'LPF.TXT' .STATUS-'UNKNOWN' ) C Write a message to the screen to remind the user of the type C of data he has specified for the binary input data file. IF (TYPE .EQ. '1*4') THEN WRITE(6,890) WRITE(6,891) ELSE IF (TYPE .EQ. 'R*4') THEN WRITE(6, 890) WRITE (6, 892) ELSE WRITE(6,893) 893 FORMAT!//' DATA TYPE ENTERED INCORRECTLY: EXECUTION ABORTED') STOP END IF END IF 890 FORMAT (//////' NOTE: THE DATA IN THE BINARY INPUT FILE MUST') 891 FORMAT (' BE OF TYPE INTEGER*4.' /) 892 FORMAT (' BE OF TYPE REAL*4.' /) WRITE (7, 895) 895 FORMAT(' PARAMETERS USED BY LOWPASS FILTER PROGRAM LPF.EXE:') WRITE(6, 900) FILNAM WRITE(7,900) FILNAM 900 FORMAT)//' Input data will be read from ',A39) C Read the number of points in the time series, the sample C interval (in seconds), the desired cutoff frequency (6 dB down C point), and the filter slope (in dB/octave). READ(2,*) NFT READ(2,*) DELT READ(2,*) CUTOFF READ(2,*) SLOPE WRITE(7,901) DELT WRITE(6,902) DELT*1000000.D0 901 FORMAT(' The sample interval is',F12.7,' seconds.') 902 FORMATC The sample interval is',F8.1,' microseconds.'/) WRITE(7,904) 1,NCHAN 904 FORMAT(' Channels',12,' through',12,' will be read.') IF (NRW .GT. 0) THEN WRITE (7, 907) 907 FORMAT(' The following channels WILL NOT be filtered:') WRITE(7,908) (RWCHAN(I),1-1,NRW) WRITE(7,909) NCHAN-NRW 908 FORMAT(1613) 909 FORMAT(' The remaining',12,' channels WILL be filtered.') ELSE WRITE(7,910) NCHAN 910 FORMAT(' All',13,' channels will be filtered.') END IF WRITE(7,920) WRITE(7,922) CUTOFF WRITE(7,924) SLOPE 920 FORMAT(/' Lowpass f i l t e r specifications:') 922 FORMAT(' Cutoff (6 dB down point):',F9.2,' Hz') 924 FORMAT(' Slope (in dB/octave) :',F9.2) C Read in the input traces. IF (TYPE .EQ. '1*4') THEN DO 28 I-1,NFT READ(1,REC-I) (INTVAR(K) ,K-1,NCHAN) DO 27 K-l,NCHAN 27 Y(I,K)-DBLE(INTVAR(K) ) 28 CONTINUE ELSE DO 30 I-1,NFT READ(1,REC-I) (RVAR(K) ,K-l,NCHAN) DO 29 K-l,NCHAN 29 Y(I,K)-DBLE(RVAR(K)) 30 CONTINUE END IF WRITE(6,980) 980 FORMAT(' Finished reading data',/) C Call a subroutine to construct the frequency domain lowpass C f i l t e r . The f i l t e r coefficients are stored in array FlLT(i) C in an order that allows the elements of FILT to be multipled C directly with the elements of the FFT array returned by subroutine C REALFT. CALL LP (FILT,DELT,CUTOFF, SLOPE,NFT) C Determine which traces are to be filtered and f i l t e r these C traces one by one using subroutine FRQFIL. DO 40 I-1,NCHAN FLAG—0 DO 20 J-l,NRW IF (I .EQ. RWCHAN(J)) THEN WRITE(6,983) I FORMAT( ' Channel',12,' will not be filtered.') FLAG—1 END IF - CONTINUE IF (FLAG .EQ. 0) THEN WRITE(6,984) I 984 FORMAT(' Beginning to fi l t e r channel', 13,'.') CALL FRQFIL(Y(1,1),FILT,NFT) END IF 40 CONTINUE C Write the processed data to a multiplexed binary f i l e of type REAL*4. DO 50 I-1,NFT 50 WRITE(3,REC-I) (SNGL (Y (I, J) ), J - l , NCHAN) CLOSE(UNIT-1) CLOSE(UNIT—2) CLOSE(UNIT-3) CLOSE(UNIT—7) STOP END C ******************************************************************** C Subroutine LP - a subroutine that returns an array FILT(i) C holding the f i l t e r coefficients for a frequency domain lowpass C f i l t e r . The coefficients are ordered in such a way the the elements C of FILT may be multiplied with the corresponding elements of the C FFT array returned by subroutine REALFT. C Q ******************************************************************** SUBROUTINE LP(FILT,DELT,CUTOFF, SLOPE,NFT) INTEGER NFT,I,J REAL*8 FILT(NFT),DELT,CUTOFF,SLOPE,P,DELF,FREQ P-SLOPE/6.0D0 DELF-1.0D0/(DELT*DBLE(NFT)) C Make the frequency domain f i l t e r in a form such that i t can be C multiplied directly with the FFT returned by subroutine REALFT. FILT(1)-1.0D0 FREQ—DBLE(NFT/2)*DELF FILT(2)-(1.0D0/(1.0D0+(FREQ/CUTOFF)**P)) DO 10 1-3,(NFT-1),2 FREQ—DBLE((1-1)/2)*DELF FILT(I)-(1.0D0/(1.0D0+(FREQ/CUTOFF)**P)) J-I+l FILT(J)-FILT(I) 10 CONTINUE RETURN END c ******************************************************* ************* C Subroutine FRQFIL - A subroutine that calculates the FFT of a real C signal, multiplies i t by a frequency domain fi l t e r , and performs an C inverse FFT on the product to return a filtered signal. The FFT and C inverse FFT are performed by subroutine REALFT from the NUMERICAL C RECIPES library. The trace to be filtered (DATA) must be real and C the number of points (NFT) in the trace must be a power of two C (these restrictions are imposed by subroutine REALFT). C Q ******************************************************************** SUBROUTINE FRQFIL(DATA,FILT,NFT) INTEGER NFT,N,I REAL*8 DATA(NFT),FILT(NFT) C Calculate the FFT of the real time series. N-NFT/2 CALL REALFT(DATA,N, 1) C Apply the frequency domain f i l t e r . DO 20 I-1,NFT 20 DATA(I)-DATA(I)*FILT(I) C Perform the inverse FFT and apply a scaling factor to return C the filtered data to the time domain. CALL REALFT(DATA,N,-1) DO 30 I-1,NFT 30 DATA (I)-DATA (I)/DBLE (N) RETURN END Q ************************************* ************* * ******** ********* SUBROUTINE REALFT (DATA, N, ISIGN) REAL*8 WR,WI,WPR,WPI,WTEMP,THETA REAL*8 DATA(*) THETA-6.28318530717959DO/2.0DO/DBLE(N) . Cl-0.5 IF (ISIGN.EQ.1) THEN C2—0.5 CALL FOURl(DATA,N,+l) ELSE C2-0.5 THETA—THETA ENDIF WPR— 2. 0D0*DSIN (0 . 5D0*THETA) **2 WPI-DSIN(THETA) WR-1.OD0+WPR WI-WPI N2P3-2*N+3 DO 11 1-2,N/2+1 11- 2*1-1 12- I1+1 13- N2P3-I2 14- I3+1 WRS-SNGL (WR) WIS-SNGL(WI) H1R-C1* (DATA (ID+DATA (13) ) H1IK:1* (DATA(I2)-DATA(I4) ) H2R—C2* (DATA (12)+DATA (14)) H2I-C2* (DATA(I1)-DATA(I3) ) DATA(II)—H1R+WRS*H2R-WIS*H2I - DATA(I2)-H1I+WRS*H2I+WIS*H2R DATA(13)-H1R-KRS*H2R+WIS*H2I DATA(14)—H1I+WRS*H2I+WIS*H2R WTEMP-WR WR-WR*WPR-WI*WPI+WR WI-WI *WPR+WTEMP *WP I+WI 11 CONTINUE IF (ISIGN.EQ.l) THEN HlROATA(l) DATA(1)—H1R+DATA(2) DATA(2)-H1R-DATA(2) ELSE HlR-OATA(l) DATA(1) -CI*(H1R+DATA(2)) DATA (2) -CI* (H1R-DATA (2) ) CALL FOURl(DATA,N,-l) ENDIF RETURN END Q ******************************************************************** SUBROUTINE FOUR1 (DATA, NN, ISIGN) REAL*6 WR,WI,WPR,WPI,WTEMP,THETA REAL*8 DATA(2*NN) N-2*NN J-l DO 11 I-1,N,2 IF(J.GT.I)THEN TEMPR-DATA(J) TEMPI-OATA(J+l) DATA(J)-DATA(I) DATA(J+1)-DATA(I+1) DATA(I)—TEMPR DATA(1+1)-TEMPI ENDIF M-N/2 1 IF ((M.GE.2) .AND. (J.GT.M) ) THEN J-J-M M-M/2 GO TO 1 ENDIF J-J+M 11 CONTINUE MMAX-2 2 IF (N.GT.MMAX) THEN ISTEP-2*MMAX THETA-6.28318530717959D0/(ISIGN*MMAX) WPR—2.D0*DSIN(0. 5D0*THETA) **2 WPIOSIN (THETA) HR-1.D0 WI-0.D0 DO 13 M-l.MMAX,2 DO 12 I-M,N,ISTEP ' J-I+MMAX TEMPR-SNGL(WR)*DATA(J)-SNGL(WI)*DATA(J+l) TEMPI-SNGL(WR)*DATA(J+l)+SNGL(WI)*DATA(J) DATA(J)-DATA(I)-TEMPR DATA(J+l)—DATA(1+1)-TEMPI DATA(I)-DATA(I)+TEMPR DATA(I+1)-DATA(1+1)+TEMPI 12 CONTINUE WTEMP-WR WR-WR*WPR-WI*WPI+WR WI-WI *WPR+WTEMP *WPI+WI 13 CONTINUE MMAX-ISTEP GO TO 2 ENDIF RETURN END trig4.3ex.for (February 15,1991 8:32pm) Q ******************************************************************** C TRIG43EX.FOR - written by K.E. Butler C The main program and subroutine STACK were o r i g i n a l l y C completed i n August, 1990 and l a s t modified Feb.9/91. C Subroutine LYTKB2 i s a modified version of subroutine LYT C obtained from T.J. Ulrych. Subroutines CORREL, TWOFFT, C REALFT, and F0UR1 are double p r e c i s i o n versions of the C Numerical Recipes (Press et al.,1989) subroutines of the C same name. C C DESCRIPTION -This i s a robust t r i g g e r i n g and stacking algorithm designed C for use on data sets that have a tri g g e r channel and C up to three a u x i l i a r y data channels. In p a r t i c u l a r , the C program has been designed to stack p i e z o e l e c t r i c data C sets that have been acquired using a r e p e t i t i v e source, C and a geophone signal for t r i g g e r i n g . Traces are C aligned for stacking by cross-correlating the envelope C of the i t h t r i g g e r (geophone) signal with the envelope C of a reference tri g g e r s i g n a l . Signal i i s then shi f t e d C forward or backward i n time to where the c o r r e l a t i o n C i s highest. C C ** - D i f f e r s from TRIG42EX.FOR only i n that the order i n which some C variables are declared has been changed, and some MS FORTRAN v.S.OO C extensions to the FORTRAN77 standard have been incorporated to C reduce the s i z e of the executable f i l e and increase speed of C execution. These extensions are: C (i) use of the ALLOCATABLE and HUGE attributes and the ALLOCATE C statement to allocate the storage space for arrays S(I), and E(I) C at run time rather than at compile time. This reduced the s i z e of C the *.EXE f i l e by well over 100Kb. C ( i i ) use of the NEAR at t r i b u t e and INTERFACE TO statement to instruct C the compiler to place subroutines LYTKB2 and STACK i n the same C memory segment as the main program. This appears to have increased C the speed of execution by about 10% (when the /Gt4096 option i s C used at compilation). C ** NOTE: The MS FORTRAN documentation warns that the NEAR and HUGE C attributes are meaningful only on processors, such as the C 8086, which have a segmented architecture. Therefore, t h i s C source code i s NOT AS PORTABLE AS TRIG42EX.FOR. c c **** . DIFFERS FROM TRIG41EX.FOR i n that the input parameter f i l e has C been renamed TRIG42EX.PRM and i n that the format of some of the C output f i l e s has been changed. The output f i l e s generated by C t h i s program are explained i n comments below. C c .«.* . DIFFERS FROM TRIG40EX.FOR i n that subroutine LYTKB (renamed C LYTKB2 i n t h i s program) has been made more e f f i c i e n t ; i t has C - been changed to c a l l REALFT rather than FOUR1 i n order to C calculate the FFT of each extracted tri g g e r trace. This change C makes the program run about 9% faster than TRIG40EX on the '386. C c **** . H 0 W THIS PROGRAM DIFFERS FROM VERSIONS OF TRIG*EX.FOR'OLDER C THAN TRIG40EX.FOR: C (i) The e l e c t r i c a l traces are extracted from the large input data C array E(*,*) AFTER the correct tri g g e r point has been found and C are stacked within the main program (subroutine STACKE has been C deleted). This change makes the program run about 5% faster on C the '386 (with '387) computer. C ( i i ) The f i n a l seismic stack FSTACK(*) i s now calculated i n C exactly the same way as the f i n a l e l e c t r i c a l stacks (ie. C FSTACK(*) i s now calculated to be the average of a l l traces C stacked rather than the average of a l l the intermediate stacks.) C Subroutine LYTKB2 i s a modified version of subroutine LYT C C INPUT PARAMETERS: C C FILNAM - name of the binary f i l e containing the input data C RLEN - number of bytes i n one record i n the input data f i l e C NELEC - number of a u x i l i a r y (usually e l e c t r i c a l ) channels to be C stacked i n addition to the t r i g g e r channel C DATLEN - number of records to be read from the binary input data f i l e C each time a swath of data i s read from disk C NLOOP - number of times the input data i n memory i s replaced by C reading a new swath of data from disk C LEVEL - the amplitude that the geophone signal must r i s e above i n C order to r e g i s t e r as an i n i t i a l t r i g g e r point C FSTART - The record number that i s to mark the s t a r t of data to be C stacked i n the input data f i l e . The f i r s t trigger signal C encountered i n the data following FSTART becomes the C reference trigger s i g n a l . C PRE - the number of points to precede the t r i g g e r point i n each C t r i g g e r signal C MAXLAG - The maximum.number of samples that an i n i t i a l t r i g g e r point C may be s h i f t e d . If the envelope cross-correlation i s found C to peak at a lag greater than MAXLAG or less than -MAXLAG C the current tri g g e r trace i s discarded ( i t i s not added to C the stack). C SKIP - the number of samples following the current trigger point C that are skipped over by the algorithm when i t begins C searching for the next t r i g g e r point C EAHP - amplification factor by which the stacked e l e c t r i c a l channels C are to be multiplied p r i o r to being output to the ASCII output C f i l e OUTDATl (no a m p l i f i c a t i o n i s applied p r i o r to output to C the binary output f i l e OUTDAT8) C NFT - The t o t a l number of samples i n each trace extracted from the C input data for stacking. The maximum value for NFT i s 512. C ** NOTE: NFT MUST BE A POWER OF TWO IN ORDER TO SATISFY THE C ASSUMPTIONS OF THE FFT SUBROUTINE USED BY THIS PROGRAM. C Q ********************* ********************************** ************* c C OUTPUT FILES: C C OUTDATl - An ASCII f i l e l i s t i n g the f i n a l stacks of a l l channels C OUTDAT2 - An ASCII f i l e l i s t i n g the reference geophone signal and C the f i n a l stack on the geophone channel C OUTDAT3 - Comments (text) d e t a i l i n g the progress of stacking C OUTDAT7 - A binary f i l e l i s t i n g the intermediate stacks for a l l C channels. Each record contains one REAL*4 sample from each C channel. The record length i s given by 4*(no. of channels). C OUTDAT8 - A binary f i l e l i s t i n g the f i n a l stacks for a l l C channels. Each record contains one REAL*4 sample from each C channel. The record length i s given by 4*(no. of channels). C OUTDAT9 - A text f i l e that l i s t s the input parameters and gives a C b r i e f summary of the program output. C c ******************************************************************** C NOTES: C - F i n a l stacks are output for the t r i g g e r (seismic) C channel and NELEC e l e c t r i c a l channels. (The number C of e l e c t r i c a l channels NELEC i s spec i f i ed by the C user.) C - intermediate t r i gge r channel s tacks are wr i t t en to C a n output data f i l e C *** - intermediate e l e c t r i c a l Channel stacks are C wr i t t en to a binary output data f i l e C - the DC of fse t on the seismic s i g n a l i s estimated C (based on the DC offse t of the f i r s t 2*NFT seismic C channel samples) and removed before envelopes are C ca lcu la ted C - the DC of fse t for each e l e c t r i c a l t race i s C ca lcu la ted and removed before s tacking C - stacked t races are of constant f o l d along t h e i r C en t i r e length Q *************************************************************** INTERFACE TO SUBROUTINE LYTKB2[NEAR] (X,NT,NFT,ENVX) INTEGER NT,NFT REAL* 8 X, ENVX END INTERFACE TO SUBROUTINE STACK [NEAR] (X,DATA, NFT, TOTLAG) INTEGER NFT,TOTLAG REAL*8 X, DATA END COMPLEX*16 ANS(S12) REAL*8 ENVX(512),REF(512),REFENV(512) , STKDAT(512),FSTACK(512) REAL*8 X<912),ANSR(1024),CCMAX REAL*4 EAMP INTEGER*2 E [ALLOCATABLE, HUGE] ( : , : ) , S [ALLOCATABLE, HUGE](:) INTEGER*2 ED(3,512),LEVEL,RLEN, RLEN2 INTEGER EDS(3,512),FEDS(3,512), KCC(100) , TRGPT, RLOOP,TEMP, FSTART INTEGER A,B,H,I,J,K,L,NFT,NT,LAG,MAXLAG,DATLEN,AVSAMP,PRE,POST INTEGER NFTA,NFTB,NFTC,THIN,SKIP,NSAMP,NLOOP,LASTRG,START,NELEC INTEGER TOTLAG, RECNUM, RECN, SUM, MEAN, I ERR, DCOFF, OFFSET CHARACTER*40 FILNAM EQUIVALENCE (ANS,ANSR) ALLOCATE ( S(23010), E(3,23010), STAT-IERR) IF (IERR .NE. 0) THEN WRITE(6,*)'Not enough storage for data; execution aborted. ' STOP END IF OPEN (UNIT-2, F I L E - ' TRIG42EX.PRM', STATUS-' OLD' ) . READ (2 , ' (A39) ' ) FILNAM READ (2,*) RLEN OPEN (UNIT-1, FILE-FILNAM, STATUS-' OLD' , ACCESS-' DIRECT' , + FORM-'UNFORMATTED' ,RECL-RLEN) OPEN (UNIT-3,FILE-'0UTDAT1') OPEN (UNIT-4, FILE-'OUTDAT2') OPEN (UNIT-5,FILE-'OUTDAT3') READ(2,*) NELEC RLEN2-4+NELEC*4 OPEN (UNIT-7, F I L E - ' OUTDAT7' , STATUS-'UNKNOWN' , ACCESS-'DIRECT', + FORM—'UNFORMATTED',RECL-RLEN2) OPEN (UNIT-8.FILE-'OUTDAT8',STATUS-'UNKNOWN', ACCESS-'DIRECT', + FORM-'UNFORMATTED' ,RECL-RLEN2) OPEN (UNIT-9,FILE-'OUTDAT9') READ(2,*) DATLEN READ(2,*) NLOOP READ(2,*) LEVEL READ(2,*) FSTART READ(2,*) PRE READ(2,*) MAXLAG READ(2,*) SKIP READ(2,*) EAMP READ(2,*) NFT THIN-1 AVSAMP-1 NT—NFT-MAXLAG POST-NFT-PRE-1 NSAMP—INT (DATLEN/THIN) -POST-MAXLAG NFTA-NFT/2 NFTB—NFT/2+1 NFTC—NFT+1 FLAG-0 WRITE(9,600) 600 FORMAT!//' SUMMARY OF INPUT AND OUTPUT FOR PROGRAM TRIG42EX:' ) WRITE(5,610) FILNAM WRITE(6,610) FILNAM WRITE(9,610) FILNAM 610 FORMAT!//' Input data was read from f i l e ' ,A39 , /1 WRITE(9,*)'Input Parameters:' WRITE(9,620) RLEN, NELEC WRITE(9,630) DATLEN, NLOOP WRITE(9,640) LEVEL, FSTART WRITE(9,650) PRE, MAXLAG WRITE(9,660) SKIP, EAMP WRITE(9,670) NFT 620 FORMAT!' RLEN: ' , 17 ,12X, 'NELEC: ' ,17) 630 FORMAT(' DATLEN:',17,12X,'NLOOP: ' ,17) 640 FORMAT!' LEVEL: ' , 17 ,12X, 'FSTART: ' , 17 ) 650 FORMAT(' PRE: ' ,17 ,12X, 'MAXLAG:' ,17) 660 FORMAT(' SKIP: ' ,17 ,12X, 'EAMP: ' ,F9 .3 ) 670 FORMAT(' NFT: ' , 1 7 , / / ) WRITE(9,*)'Output Summary:' WRITE(9,680) 680 FORMAT!/' Number of Traces i n the Intermediate S tacks : ' ) WRITE(9,*)' Stack No. of Traces' DO 5 I - l . N F T DO 3 J - l ,NELEC 3 EDS(J , I ) -0 FSTACK(I)-0.0D0 5 CONTINUE C Find the f i r s t s igna l and use i t as a reference t race . J-0 DCOFF-0 DO 10 I—FSTART,(FSTART+NFT*2*THIN),THIN J-J+l READ(1,REC-I) S(J) S(J)-S(J)-2048 DCOFF-DCOFF+S(J) 10 CONTINUE DCOFF-INT(DCOFF/J) OFFSET-DCOFF+2 048 DO 15 1-1, (NFT*2*THIN) , THIN S(I)-S(I)-DCOFF CONTINUE I-PRE+MAXLAG+1 IF (S(I) .GT. LEVEL) THEN A-0 B-0 DO 22 J-(I-AVSAMP),(1-1) A-A+S(J) DO 24 J-(I+1),(I+AVSAMP) B-B+S(J) A-INT(A/AVSAMP) B-INT(B/AVSAMP) IF ((A.LT.LEVEL) .AND. (B.GT.LEVEL)) THEN L-0 DO 26 J-(I-PRE),(I+POST) L-L+l REF(L)-DBLE(S(J)) CONTINUE WRITE(6,*) ' REFERENCE TRACE FOUND' WRITE(5,*)'REFERENCE TRACE FOUND' I-WFT+PRE+MAXLAG FLAG-1 ELSE I-I+l END IF ELSE I-I+l END IF IF (I .LT. (NFT+PRE+MAXLAG) ) THEN GO TO 20 ELSE IF (FLAG .EQ. 0) THEN WRITE(6,*) 'PROBLEM!I NO TRIGGER FOUND WITHIN',NFT,' PTS.' WRITE(6,*)'OF START OF RECORD.' STOP END IF END IF Find the envelope of the reference trace assuming the l a s t MAXLAG samples i n the trace are zeros. CALL LYTKB2(REF,NT,NFT,REFENV) Zero the l a s t MAXLAG samples i n the reference envelope. DO 30 I-NT+1,NFT REFENV(I)-0.0D0 CONTINUE Each time the following outer loop i s performed the program READS a new swath of DATLEN/THIN samples from disk, finds the t r i g g e r points i n the swath, extracts traces from each channel at the t r i g g e r points, and adds the traces to t h e i r respective stacks. An intermediate stack i s created and output each time the outer loop i s executed. RECNUM-0 K-0 DO 60 RLOOP-l,NLOOP DO 32 I-1,NFT STKDAT(I)-0.0D0 DO 33 J-1,NELEC 33 EDS(J,I)-0 32 CONTINUE C Calculate the record number corresponding to the s t a r t of the C swath of data to be read from disk. IF (RLOOP .EQ. 1) THEN START-FSTART*THIN ELSE START-INT((LASTRG+SKIP-PRE-MAXLAG)*THIN) END IF C Read a swath of DATLEN/THIN samples from disk. WRITE(6,*)'READING FROM DISK: LOOP*', RLOOP WRITE (5,*)'READING FROM DISK; LOOP*',RLOOP J-0 DO 35 I-START,(START+DATLEN),THIN J-J+l READ(1,REC-I) S(J), (E(L,J),L-1,NELEC) S(J)-S(J) -OFFSET DO 34 L—1,NELEC 34 E(L,J)-E<L, J)-2048 35 CONTINUE C Search through the data and f i n d t r i g g e r points using C simple c r i t e r i a based on signal l e v e l and slope. KCC(RLOOP)—0 FLAG-0 I-PRE+MAXLAG+1 40 IF (S(I) .LE. LEVEL) THEN I-I+l ELSE TEMP-I+INT (START/THIN) RECN-I*THIN+START WRITE(6,*)'TEST 1 PASSED AT RECORD',RECN WRITE(5,*)'TEST 1 PASSED AT SAMPLE*',TEMP,' REC' ,RECN A-0 B-0 DO 42 J-(I-AVSAMP),(1-1) 42 A-A+S(J) DO 44 J—(1+1),(I+AVSAMP) 44 B-B+S(J) A-INT(A/AVSAMP) B-INT(B/AVSAMP) IF ((A.LT.LEVEL) .AND. (B.GT.LEVEL)) THEN K-K+l TRGPT—TEMP L-0 DO 46 J-(I-PRE-MAXLAG) , (I+POST+MAXLAG) L-L+l 1 X(L)-DBLE(S(J)) 46 CONTINUE C The envelope of the trace extracted above Is computed C and cross-correlated with the envelope of the C reference trace. CALL LYTKB2(X(MAXLAG+1),NT,NFT, ENVX) DO 56 J-NT+1,NFT 56 ENVX(J)-0.0D0 CALL CORREL (REFENV, ENVX,NFT, ANS,ANSR) C The lag (LAG) coresponding to the maximum correlation is C now found. CCMAX-ANSR(l) LAG-0 DO 57 J-2,NFTA IF (ANSR(J) .GT. CCMAX) THEN CCMAX-ANSR(J) LAG-J-1 END IF 57 CONTINUE DO 58 J-OTTB,NFT IF (ANSR(J) .GT. CCMAX) THEN CCMAX-ANSR(J) LAG-J-NFTC END IF 58 CONTINUE WRITE(6,700) K, LAG, CCMAX WRITE(5,700) K,LAG,CCMAX 700 FORMAT(' TRACE',16,' TRIGGERED, LAG-',15,', CCMAX-',F13.1) C Stack the trace if its maximum cross correlation occurs at C a lag of no more than MAXLAG samples. IF (ABS (LAG) .GT. MAXLAG) THEN WRITE(6, *)' LAG IS GREATER THAN THE MAX. ALLOWABLE LAG.' WRITE(6,*)'THIS TRACE WILL NOT BE STACKED.' WRITE(5,*)'LAG IS GREATER THAN THE MAX. ALLOWABLE LAG.' WRITE (5,*)'THIS TRACE WILL NOT BE STACKED.' K-K-l ELSE TOTLAG-MAXLAG-LAG CALL STACK(X.STKDAT,NFT,TOTLAG) C Extract traces from each electrical channel, calculate and C remove their DC offsets (mean values), and add the traces to C their respective stacks. DO 585 H—1,NELEC L-(I-PRE-LAG) SUM-0 DO 583 J-1,NFT ED (H, J) -E (H, L) SUM-SUM+INT(ED(H,J)) L-L+l 583 CONTINUE MEAN - INT (SUM/NFT) DO 584 J-1,NFT 584 EDS (H, J) —EDS (H, J) +INT (ED (H, J) )-MEAN 585 CONTINUE KCC(RLOOP)-KCC(RLOOP)+1 END IF I-I+SKIP ELSE I-I+l END IF END IF IF (I .LT. NSAMP) GO TO 40 IASTRG-TRGPT+LAG WRITE(6, *) 'NORMALIZING INTERMEDIATE STACK',RLOOP DO 59 I-1,NFT FSTACK(I)-FSTACK(I)+STKDAT(I) STKDAT(I)-STKDAT(I)/DBLE(KCC(RLOOP)) DO 586 J—1,NELEC 586 FEDS(J,I)-FEDS(J,I)+EDS(J,I) 59 CONTINUE C Write the current intermediate stacks to a binary output file. DO 595 I-1,NFT ' RECNUM—RECNUM+1 WRITE(7,REC-RECNUM) FLOAT (STKDAT (I) ), + (FLOAT (EDS (L, I) ) /FLOAT (KCC (RLOOP) ), L-l,NELEC) 595 CONTINUE WRITE(9,800) RLOOP,KCC(RLOOP) 60 CONTINUE 800 FORMAT(16,6X,17) C Divide the final stacks by the total number of traces stacked and C write the final stacks to various output files. WRITE(6,*)'CALCULATING FINAL STACKS' DO 65 I-1,NFT FSTACK(I)-FSTACK(I)/DBLE(K) 65 CONTINUE DO 80 I-1,NFT WRITE(4,900) REF(I),FSTACK(I) WRITE(3,900) FLOAT(FSTACK(I)), (FLOAT(FEDS(L,I))/FLOAT(K)*EAMP, + L-l,NELEC) WRITE(8,REC-I) FLOAT(FSTACK(I)), (FLOAT(FEDS(L,I))/FLOAT (K), + L-l,NELEC) 80 CONTINUE 900 FORMAT (5F14. 4) WRITE(6,*) 'FINISHED!' WRITE(5,*j 'FINISHED!' WRITE(9,910) NELEC WRITE(9,920) RLEN2 WRITE(9,930) K 910 FORMAT!//' The trigger channel and',13,' electrical channels were +stacked.') 920 FORMAT!' The record length in the binary output files is',14,'.') 930 FORMAT (' Each final stack is the average of,16,' traces.') CLOSE (UNIT-1) CLOSE (UNIT-2) CLOSE (UNIT-3) CLOSE (UNIT-4) CLOSE (UNIT-5) CLOSE (UNIT-7) CLOSE (UNIT-8) CLOSE (UNIT-9) STOP END Q ************* ************ ********************** ************** ******** C C SUBROUTINE LYTKB2 C This subroutine calculates the envelope of the scalar series x ( t ) . C F i r s t the Fourier transform of the an a l y t i c signal corresponding C to X i s found. (The an a l y t i c signal i s X(t)-i*H(X(t)) where i -C SQRT(-l) arid H(X(t)) i s the Hilbert transform of X. The Fourier C transform of the analytic signal i s obtained from the the Fourier C transform of X by se t t i n g the amplitude of the negative frequency C components equal to zero, doubling the amplitude of the po s i t i v e C frequency components and leaving the amplitude of the DC component C unchanged.) The envelope ENVX of X i s the magnitude of the C a n a l y t i c s i g n a l . ( *** NOTE: This i s a more e f f i c i e n t version of C subroutine LYTKB; the FFT of the (real) scalar series X(t) i s C c a l c u l a t e d using subroutine REALFT rather than FOUR1.) SUBROUTINE LYTKB2 [NEAR] (X, NT, NFT, ENVX) REAL*8 X(NT),ENVX(NFT) REAL*8 DATA(1024) COMPLEX*16 CDUMMY(1) N^JFT*2 DO 100 J-l,NT 100 DATA(J)-X(J) DO 105 J-(NT+1),NFT 105 DATA(J)-0.D0 CALL REALFT (CDUMMY,NFT/2, 1,DATA) DATA(2)-0.DO DO 200 J-3,NFT 200 DATA(J)-2.0D0*DATA(J) DO 300 J-(NFT+1),N 300 DATA(J)-0.D0 CALL FOUR1(DATA,NFT,-1) DO 400 J-1,N 400 DATA (J)-DATA (J)/DBLE (NFT) DO 500 J-l,NFT W * 2 - l 500 ENVX (J) -DSQRT (DATA(I) **2+DATA (1+1) **2) RETURN END c ********************************************************************* C Subroutine STACK - This subroutine s h i f t s and stacks traces extracted C from the trigger channel. SUBROUTINE STACK [NEAR] (X, DATA, NFT, TOTLAG) INTEGER I,J,NFT,TOTLAG REAL*8 DATA(*),X(*) DO 10 I-1,NFT J-I+TOTLAG DATA (I) OATA (I)+X (J) 10 CONTINUE RETURN END SUBROUTINE CORREL (DATA1, DATA2,N, ANS, ANSR) PARAMETER (NMAX-1024) REAL*8 DATA1 (N) ,DATA2 (N) ,RFFT (2*NMAX) ,ANSR(2*N) COMPLEX*16 FFT(NMAX) ,ANS(N) EQUIVALENCE (FFT,RFFT) CALL TWOFFT (D AT A l , DATA2, FFT, ANS, N, RFFT) N02-DBLE(N)12.0 DO 11 I-l,N/2+l ANS (I)—FFT(I)*DCONJG(ANS(I))/N02 11 CONTINUE ANS (l)-DCMPLX (DBLE (ANS (1) ) , DBLE (ANS (N/2+1) ) ) CALL REALFT (ANS,N/2,-1, ANSR) RETURN END SUBROUTINE TWOFFT(DATA1,DATA2,FFTl,FFT2,N,RFFT) REAL*8 DATA1(N),DATA2(N),RFFT(2*N) COMPLEX*16 FFTl(N), FFT2(N),H1,H2,C1,C2 C1-DCMPLX(0.5,0.0) C2-DCMPLX(0.0,-0.5) DO 11 J-1,N FFT1(J)-DCMPLX(DATA1(J) ,DATA2 (J) ) 11 CONTINUE CALL FOURl(RFFT,N,l) FFT2(1)-DCMPLX(DIMAG(FFTl(1)),0.0D0) FFTl(1)-DCMPLX(DBLE(FFTl(1)),0.0D0) N2-N+2 DO 12 J-2,N/2+1 Hl-Cl*(FFTl(J)+DCONJG(FFTl(N2-J))) H2-C2*(FFTl(J)-DCONJG(FFTl(N2-J))) FFTl(J)—HI FFTl(N2-J)-DCONJG(HI) FFT2(J)-H2 FFT2(N2-J)-DCONJG(H2) 12 CONTINUE RETURN END SUBROUTINE REALFT (CDUMMY, N, ISIGN, DATA) COMPLEX*16 CDUMMY(*) REAL*8 WR,WI, WPR, WPI, WTEMP, THETA REAL*8 DATA(*) THETA-6.28318530717959D0/2.0D0/DBLE(N) Cl-0.5 IF (ISIGN.EQ.l) THEN C2—0.5 CALL FOURl(DATA,N,+l) ELSE C2-0.5 THETA—THETA ENDIF WPR—2. 0D0*DSIN (0 . 5D0*THETA) **2 WPI-OSIN(THETA) WR-1.0D0+WPR I-PI N2P3-2*N+3 DO 11 1-2,N/2+1 11- 2*1-1 12- I1+1 13- N2P3-I2 14- I3+1 WRS-SNGL(WR) WIS-SNGL(WI) H1R-C1* (DATA(I1)+DATA(I3) ) H1I-C1* (DATA(I2)-DATA(I4) ) H2R—C2* (DATA (12) +DATA (14) ) H2I-C2* (DATA(I1)-DATA(I3) ) DATA(II)-H1R+WRS*H2R-WIS*H2I DATA(12)-H1I+WRS*H2I+WIS*H2R DATA(I3)-H1R-WRS*R2R+WIS*H2I DATA ( 1 4 ) — H1I+WRS*H2I+WIS*H2R WTEMP-WR WR-WR*WPR-WI*WPI+WR WI-WI*WPR+WTEMP *WPI+WI 11 CONTINUE IF (ISIGN.EQ.l) THEN HlR-DATA(l) DATA (1)-H1R+DATA (2) DATA(2)—H1R-DATA(2) ELSE HlR-DATA(l) DATA(1)-CI*(H1R+DATA(2)) DATA (2) "CI* (H1R-DATA (2) ) CALL FOURl(DATA,N,-l) ENDIF RETURN END SUBROUTINE FOUR1 (DATA,NN, ISIGN) REAL*8 WR,WI,WPR,WPI,WTEMP, THETA REAL*8 DATA(2*NN) N-2*NN J - l DO 11 I-1,N,2 IF(J.GT.I)THEN TEMPR—DATA(J) TEMPI-DATA(J+l) DATA(J)-DATA(I) DATA(J+1)-DATA(1+1) DATA(I)-TEMPR DATA(1+1)-TEMPI ENDIF M-N/2 1 IF ( (M.GE.2) .AND. (J.GT.M) ) THEN J-J-M M-M/2 GO TO 1 ENDIF J-J+M 11 CONTINUE MMAX-2 2 IF (N.GT.MMAX) THEN ISTEP-2*MMAX THETA—6.28318530717959D0/ (ISIGN*MMAX) WPR—2.D0*DSIN(0. 5D0*THETA) **2 WPI-t>SIN (THETA) WR-1.D0 WI-O.DO' DO 13 M-1,MMAX,2 DO 12 I-M,N,ISTEP J-I+MMAX TEMPR—SNGL(WR)*DATA(J)-SNGL(WI)*DATA(J+l) TEMPI—SNGL(WR)*DATA(J+l)+SNGL(WI)*DATA(J) DATA(J)-DATA(I)-TEMPR DATA(J+l)-DATA(1+1)-TEMPI DATA(I) -OATA(I)+TEMPR DATA(1+1)-DATA(1+1)+TEMPI CONTINUE WTEMP-WR WR-WR*WPR-WI*WPI+WR WI-WI*WPR+WTEMP*WPI+WI CONTINUE MMAX-ISTEP GO TO 2 ENDIF RETURN trig42ex.prm (September 20,1991 3:48am) d:\data\b37.dat 8 (RLEN - record length i n input data f i l e ) 3 (NELEC - no. of e l e c t r i c a l channels) 970 (DATLEN) 6 (NLOOP) 24 (LEVEL) 1 (FSTART) 225 (PRE) 70 (MAXLAG) 200 (SKIP) 1.0 (EAMP) 512 ( NFT - MUST BE A POWER OF 2 !!!!!!!) to 1 Q *********************************************************** 124 2 C Program PZM.FOR - w r i t t e n by K.Butler, l a s t r e v i s e d OCT 9/91 3 C • 4 C T h i s program c a l c u l a t e s the frequency and time 5 C domain r e p r e s e n t a t i o n s of the e l e c t r i c p o t e n t i a l at p o i n t s i n 6 C a co n d u c t i v e whole space due t o the presence of a d i s t r i b u t i o n 7 C of p i e z o e l e c t r i c spheres which are p o l a r i z e d by a s p h e r i c a l l y 8 C spreading e l a s t i c wave. C a l c u l a t i o n s are performed i n the 9 C frequency domain and time s e r i e s are obtained using the FFT 10' C subroutine DFO0R2 from the MTSG l i b r a r y at UBC. 11 C I t i s assumed t h a t the e l a s t i c wave i s compressional and i t s 12 C p a r t i c l e d i s p l a c e m e n t i s given by the f u n c t i o n 13 C 14 C u ( r , t ) = p(t) [*] H(T)*(G/r)*exp(-i*wn/sqrt(2)*T)*sin(wn*T) . 15 C 16 C where [*] i n d i c a t e s c o n v o l u t i o n 17 C p ( t ) i s an a r b i t r a r y pressure f u n c t i o n (p(t)=0,t<0) 18 C H(t) = 0, t<0 i s the Heaviside step f u n c t i o n 19 C 1, t>0 20 C r i s the r a d i a l d i s t a n c e from the se i s m i c source 21 C T = t - ( r - a ) / v 22 C wn = 2 s q r t ( 2 ) * v / 3 a i s a s c a l a r constant 23 C G = (a**2)/(2sqrt(2)*mu) i s a s c a l a r constant 2 4 C mu i s the shear modulus 25 C i i s the square root of (-1) 26 C t i s time 27 C v i s the v e l o c i t y of the s e i s m i c wave . 28 C a i s the radius of the c a v i t y i n which the 29 C pre s s u r e i s a p p l i e d . 30 C The above ex p r e s s i o n i s the f a r - f i e l d approx. t o the p a r t i c l e 31 C displacement caused by a the a p p l i c a t i o n of an a r b i t r a r y 32 C p r e s s u r e f u n c t i o n t o the i n t e r i o r of a s p e r i c a l c a v i t y w i t h i n 33 C a s o l i d whole space. The expression was d e r i v e d by Sharpe 34 C ( J . A. Sharpe, Geophysics 7, 1942) . 35 C T h i s program assumes t h a t the pressure f u n c t i o n i s give n by 36 C 37 C p ( t ) = H(t) *pmax* < ( (Aexp(l)) **2)/4) * (t**2) * (exp(-At) ) 38 C where A = 2*wn i s a s c a l a r c onstant. 39 C A l l of the spheres are assumed to have the.same volume and 40 C p i e z o e l e c t r i c t e n s o r d ( i j k ) . These values and other parameters 41 C are s p e c i f i e d by the user i n a parameter f i l e c a l l e d PZM.PRM. 42 C 43 C **** This program gives the user the option 44 C of s p e c i f y i n g the co o r d i n a t e s of each sphere i n d i v i d u a l l y 45 C or of p r o v i d i n g i n f o r m a t i o n t h a t w i l l i n s t r u c t s u b r o u t i n e 46 C GEOMPZ to set up a p l a n a r p i e z o e l e c t r i c body made up of a 4 7 C v a r i a b l e number of spheres 48 C 4 9 Q **************************************************************** 50 C 51 INTEGER I,J,NODMAX,NNODES,NNODEX,NNODEY,NS,SPHMAX,NT,NF 52 INTEGER NNZERO,NFMAX,NODIM,RECNUM,ISPEC1,ISPEC2 53. PARAMETER (NODMAX=150,SPHMAX=2000,NTMAX=256,NFMAX=NTMAX/2+l) 54 INTEGER ITM(NODMAX),IFM(NODMAX) 55 REAL*8 X(NODMAX),Y(NODMAX),Z(NODMAX),MAXT(NODMAX),MAXF(NODMAX ) 56 REAL*8 XP(SPHMAX),YP(SPHMAX),ZP(SPHMAX),W(NFMAX),F(NFMAX) 57 REAL*8 PRT (NTMAX) , U (NTMAX) , VP (NTMAX) , POT (NTMAX) , D (3, 6) 58 REAL*8 PI,SPHRAD,A,V,R,WN,WNR2,DT,DW,MU,G,TN,DIST,TEMP 59 REAL*8 XS,YS,ZS,C,RP,PERMIT,RES,VOL,XD,YD,ZD,RESPER,PMAX 60 COMPLEX*16 CP0T(N0DMAX,NFMAX),CSIGMA(6,NFMAX),P(3,NFMAX) 61 COMPLEX*16 UF(NFMAX),VPF(NFMAX),POTF(NFMAX),PRW(NFMAX) 62 COMPLEX*16 CWPRW(NFMAX),CPP(NFMAX) 63 COMPLEX*16 CI,C2,CGWNR,CE,CPDR,CTEMP,CSCALE 64 CHARACTER*16 FNAME1,FNAME2,FNAME3,FNAME4,FNAME7 65 EQUIVALENCE (U,UF),(VP,VPF),(POT,POTF) 66 FNAME1='PZM.PRM" 67 FNAME2='SPHERES.PRM' 68 FNAME3='GRID.PRM' 69 FNAME4=1PZM.TXT' 70 FNAME7=1 MAG.OUT' 71 72 OPEN(UNIT=1,FILE=FNAME1) 73 OPEN(UNIT=2,FILE=FNAME2) 74 OPEN(UNIT=3,FILE=FNAME3) 75 OPEN(UNIT=4,FILE=FNAME4) 76 OPEN(UNIT=7,FILE=FNAME7) 77 OPEN(UNIT=8,ACCESS='DIRECT',FORM='UNFORMATTED', 78 + STATUS= 1 UNKNOWN',RECL=4) 79 C 80 C I n i t i a l i z e the p i e z o e l e c t r i c moduli by s e t t i n g them a l l 81 'c equal t o zero. 82 C 83 DO 4 1=1,3 84 DO 3 J=l, 6 85 D(I,J)=0.D0 86 3 CONTINUE 87 P fl 4. c CONTINUE j O O 89 c Read some u s e r - s p e c i f i e d model parameters from an input f i l e 90 c 91 READ(1,*) NT 92 READ(1,*) DT 93 READ(1,*) V 94 READ(1,*) MU 95 READ(1,*) PMAX 96 READ(1,*) A 97 READ(1,*) RES 98 READ(1,*) DIST 99 READ(1,*) SPHRAD 100 READ(1,*) NNZERO 101 DO 5 1=1,NNZERO 102 5 READ(1,*) J,K,D(J,K) 103 READ(1,*) DMAX 104 READ(1,*) ISPEC1 105 READ(1,*) ISPEC2 106 C 107 C Before proceeding any f u r t h e r check t h a t the a r r a y 108 C dimension NFMAX i s at l e a s t as l a r g e as the number of 109 C f r e q u e n c i e s f o r which r e s u l t s are t o be generated. 110 c 111 NF = NT/2 + 1 112 IF (NFMAX .LT. NF) THEN 113 WRITE(6,*)'Array dimension NFMAX i s t o o 'small.' 114 WRITE(6,*)'Increase NFMAX to at l e a s t ' , N F , ' . • 115 STOP 116 END IF 117 c 118 c S p e c i f y or c a l c u l a t e some constants t h a t are t o be 119 c used f r e q u e n t l y i n the modelling c a l c u l a t i o n s . 120 c 121 PI = 3.1415927D0 122 PERMIT = 8.85D-12 123 VOL = (4.D0/3.D0)*PI*SPHRAD**3 124 WN = (2.D0*DSQRT(2.D0)*V)/(3.D0*A) 125 DW = (2.D0*PI)/(NT*DT) 126 G = (A**2)/(2.D0*DSQRT(2.DO)*MU) 127 WNR2 = WN/DSQRT(2.D0) 128 CI = DCMPLX(WN**2,0.DO) 129 c 130 c Write the p i e z o e l e c t r i c matrix t o an output f i l e . 131 c 132 WRITE(4,850) 133 850 FORMAT(' The matrix of p i e z o e l e c t r i c moduli f o l l o w s : ' ) 134 DO 10 1=1,3 135 10 WRITE(4, 860) (D (I, J) , J=l, 6) 136 860 FORMAT(6G11.3) 137 C 138 C I f the f l a g ISPEC1 was set to zero by the user, c a l l a sub-139 C r o u t i n e t o c a l c u l a t e the co o r d i n a t e s of the p i e z o e l e c t r i c 140 C spheres. Otherwise ( i f ISPEC=1), the c o o r d i n a t e s of each 141 C sphere w i l l be read from an input f i l e l a t e r i n the program. 142 C 143 IF (ISPEC1 .EQ. 0) THEN 144 CALL GEOMPZ(XP,YP,ZP,SPHRAD,NS) 145 ELSE 146 READ(1,*) NS 147 DO 20 1=1,NS 148 20 READ(1,*) XP(I) , YP(I),ZP(I) 149 END IF ' 150 WRITE(6,*) ' 151 WRITE(6,870) NS 152 WRITE(4,870) NS 153 870 FORMAT(/' There are',16,' spheres i n the p i e z o e l e c t r i c body. 154 ) IF (NS .LT. 20) THEN 155 WRITE(4,875) 156 875 FORMAT(' The c o o r d i n a t e s of the spheres are l i s t e d below:' 157 W R I T E ( 4 , * ) ' X Y Z' 158 DO 25 1=1,NS 159 25 WRITE(4,880) XP(I),YP(I),ZP(I) 160 880 FORMAT(5F10.3) 161 WRITE(4,*)' 162 END IF 163 C 164 C I f the f l a g ISPEC2 was set t o zero by the user c a l l 165 C a su b r o u t i n e t o c a l c u l a t e the c o o r d i n a t e s of the p o i n t s 166 C at which the e l e c t r i c p o t e n t i a l i s t o be c a l c u l a t e d . 167 C Otherwise, the c o o r d i n a t e s of the p o i n t s are read from 168 C the i n p u t parameter f i l e . 169 c 170 IF (ISPEC2 .EQ. 0) THEN 171 c. CALL GE0M0B(X,Y,Z,NNODES,NNODEX,NNODEY) 172 WRITE(6,*)'ERROR-The c o o r d i n a t e s of each o b s e r v a t i o n ' 173 WRITE(6,*)'must be s p e c i f i e d i n the parameter f i l e . ' 174 WRITE(6,*)'Make sure you have set ISPEC=1.' 175 STOP 176 ELSE 177 READ (1,*) NNODES 178 DO 27 1=1,NNODES 179 27 READ(1,*) X ( I ) , Y ( I ) , Z ( I ) 180 WRITE(4,890) NNODES 181 WRITE(6,890) NNODES 182 END IF 183 890 FORMAT(/' The p o t e n t i a l w i l l be c a l c u l a t e d at",14,' p o i n t s . 184 C 185 c I n i t i a l i z e the elements i n the e l e c t r i c p o t e n t i a l a r r a y by 186 c s e t t i n g them a l l equal t o zero. 187 c 188 DO 30 1=1,NNODES 189 DO 2 9 J=1,NF 190 29 CPOT(I,J)=DCMPLX(0.DO, 0.DO) 191 30 CONTINUE 192 193 C 194 C C o n s t r u c t two a r r a y s of the temporal and ra d i a n 195 C f r e q u e n c i e s at which c a l c u l a t i o n s w i l l be performed 196 C i n the frequency domain. A l s o make a r r a y an a r r a y CPP 197 C f o r a term t h a t i s a f u n c t i o n of frequency and i s used 198 C . r e p e a t e d l y i n the c a l c u l a t i o n s t h a t f o l l o w . 199 C 200 RESPER = 2.DO/(3.D0*RES*PERMIT) 201 CTEMP = DCMPLX(VOL/(4.D0*PI*PERMIT),0.D0) 202 DO 32 1=1,NF 203 W(I) = DBLE(1-1)*DW 204 F(I) = W(I)/<2.D0*PI) 205 CPP(I) = CTEMP*DCMPLX(0.D0,W(I))/DCMPLX(RESPER, W(I)) 206 32 CONTINUE 207 C 208 c C a l l a su b r o u t i n e t o c a l c u l a t e the pre s s u r e f u n c t i o n 209 c and i t s F o u r i e r t r a n s f o r m (to be s t o r e d i n arrays PRT(i) 210 c and PRW(i) r e s p e c t i v e l y . T h e n c o n s t r u c t array CWPRW tha t 211 c i s the product of the pr e s s u r e ' s F o u r i e r transform and iw 212 c where w. i s r a d i a n frequency and i i s the sq. root of -1. 213 c T h i s a r r a y c o n t a i n s the t r a n s f e r f u n c t i o n by which the 214 c F o u r i e r t r a n s f o r m of the p a r t i c l e displacement step 215 c response must be m u l t i p l i e d i n order t o ob t a i n the F o u r i e r 216 c t r a n s f o r m of the a c t u a l p a r t i c l e displacement. 217 c 218 CALL PRESSUR(PRT,PRW,W,PMAX,WN, NT,NF, DT) 219 DO 33 1=1,NF 220 33 CWPRW(I) = DCMPLX(0.DO,W(I))*PRW(I) 221 C 222 C C a l c u l a t e the p a r t i c l e displacement and p a r t i c l e at a 223 c d i s t a n c e of DIST metres from the source. The c a l c u l a t i o n s 224 c are done i n the frequency domain and an i n v e r s e FFT i s 225 c used t o o b t a i n the time domain r e p r e s e n t a t i o n s . 226 c 227 R=DIST 228 TN=(R-A)/V 229 CGWNR = DCMPLX(G*WN/R,0.DO) 230 WRITE<7,*)' FREQUENCY DOMAIN DATA:' 231 DO 34 J=1,NF 232 CE = CDEXP(DCMPLX(0.DO,-W(J)*TN)) 233 C2 = DCMPLX(WNR2,W(J)) 234 UF(J) = CWPRW(J)*CGWNR*CE/(C1+C2*C2) 235 VPF(J) = DCMPLX(0.D0,W(J))*UF(J) 236 WRITE(7,915)F(J),CDABS(PRW(J)),CDABS(UF(J)),CDABS(VPF(J)) 237 34 CONTINUE 238 C 239 C The continuous ( a n a l y t i c ) F o u r i e r transforms UF and VF 240 C are m u l t i p l i e d by s c a l i n g f a c t o r s 1/DT and 1/(NT*NT) so 241 C t h a t the time s e r i e s returned by the IDFT w i l l have the 242 c c o r r e c t magnitudes. 243 c 244 NODIM = 1 245 CSCALE = DCMPLX(1.D0/(DT*NT*NT) , 0.DO) 246 DO 3 6 J=1,NF 247 UF(J) = UF(J)*CSCALE 248 VPF(J) = VPF(J)*CSCALE 249 36 CONTINUE 250 CALL DFOUR2(U,NT,NODIM,1,-1) 251 CALL DFOUR2(VP,NT,NODIM, 1,-1) 252 WRITE(7,*) ' TIME DOMAIN DATA:' 253 DO 38 1=1,NT 254 38 WRITE(7,908) 1000.D0*DT*DBLE(1-1),PRT(I),U(I),VP(I) 255 908 FORMAT(F9.4,3G11.3) 256 C 257 C At each of the NS p i e z o e l e c t r i c spheres present i n the whole 258 C space, c a l c u l a t e the F o u r i e r t r a n s f o r m of the p a r t i c l e 259 C displacement and c a l c u l a t e the 260 c p o t e n t i a l at each p o i n t i n the o b s e r v a t i o n g r i d due t o the 261 c sphere. 262 c 263 NSOUT=0 264 DO 70 1=1,NS 2 65 XS=XP(I) 266 YS=YP(I) 267 ZS=ZP (I) 268 R = DSQRT(XS**2 + YS**2 + ZS**2) 269 c 270 c I f the shot t o sphere d i s t a n c e R i s g r e a t e r than a maximum 271 c d i s t a n c e s p e c i f i e d by the user then do not use the sphere 272 c i n the c a l c u l a t i o n of e l e c t r i c p o t e n t i a l . 273 c 274 IF (R -GT. DMAX) THEN 275 NSOUT=NSOUT+l 276 ELSE 277 TN = (R-A) IV 278 CGWNR = DCMPLX(G*WN/R,0.DO) 279 DO 4 0 J=1,NF 280 CE = CDEXP(DCMPLX(0.DO,-W(J)*TN)) 281 C2 = DCMPLX(WNR2,W(J)) 282 . . UF(J) = CWPRW(J)*CGWNR*CE/(C1+C2*C2) 283 DO 3 9 K=l,3 284 39 P(K,J) = 0.D0 285 40 CONTINUE 286 C 287 C C a l l a subroutine t o c a l c u l a t e the frequency domain 288 C r e p r e s e n t a t i o n of the s t r e s s t e n s o r . 289 C 290 CALL STRESS(CSIGMA,UF,W,XS,YS,ZS,R,MU,V,NF) 291 C 292 C C a l c u l a t e the frequency domain p o l a r i z a t i o n v e c t o r 293 C f o r the sphere. 294 c 295 DO 50 J=1,NF 296 DO 46 K=l,3 297 DO 45 L=l, 6 298 P(K,J) = P(K,J) + DCMPLX (D (K, L) , 0 .DO) *CSIGMA (L, J) 299 45 CONTINUE 300 46 CONTINUE 301 50 CONTINUE 302 915 FORMAT(F9.1,6G11.3) 303 C 304 IF (NS .LE. 20) THEN 305 WRITE(4,910) I,R 306 WRITE(6,910) I,R 307 910 FORMAT(' Distance to sphere', 14,' :',F10.3) 308 END IF 309 C I 310 C C a l c u l a t e the e l e c t r i c p o t e n t i a l (as a f u n c t i o n of 311 C frequency) due to sphere I at each po i n t i n the 312 C o b s e r v a t i o n g r i d . The c o n t r i b u t i o n from each sphere 313 C i s summed t o c a l c u l a t e the t o t a l e l e c t i c p o t e n t i a l 314 c at each p o i n t . 315 c 316 DO 60 J=l,NNODES 317 XD = X(J)-XS 318 YD = Y(J)-YS 319 ZD = Z(J)-ZS 320 RP=DSQRT(XD**2 + YD**2 + ZD**2) 321 DO 55 K=1,NF 322 CPDR = P(1,K)*DCMPLX(XD,0.D0) + P(2,K)*DCMPLX(YD,0.D01 323 + + P(3,K)*DCMPLX(ZD,0.D0) 324 CPOT(J,K)=CPOT(J,K)+CPDR*DCMPLX(1.DO/(RP**3),0.D0) 325 55 CONTINUE 326 60 CONTINUE 327 END IF 328 IF (MOD(I+100,100) .EQ. 0) THEN 329 WRITE(6,*) I , ' SPHERES PROCESSED.' 330 END IF 331 70 CONTINUE 332 C 333 WRITE(4,930) NSOUT 334 WRITE(4,931) DMAX 335 WRITE(4,932) 336 930 FORMAT(/,14,1 spheres were l o c a t e d at d i s t a n c e s greater') 337 931 FORMAT(' t h a n ' , F l l . 4 , ' from the shopoint and were not') 338 932 FORMAT(' used i n the c a l c u l a t i o n of e l e c t r i c p o t e n t i a l . ' ) 339 C 340 DO 80 J=l,NNODES 341 DO 75 K=1,NF 342 75 CPOT(J,K) = CPOT(J,K)*CPP(K) 343 80 CONTINUE 344 C 345 C Write a header at the s t a r t of the b i n a r y output f i l e . 346 C 347 WRITE(8,REC=1) NT 348 WRITE(8,REC=2) NNODES 349 WRITE(8,REC=3) SNGL(DT) 350 C 351 C At each p o i n t i n the o b s e r v a t i o n g r i d perform an i n v e r s e 352 c FFT t o o b t a i n the e l e c t r i c p o e n t i a l as a f u n c t i o n of time 353 c The frequency domain data i s m u l t i p l i e d by the s c a l i n g 354 c f a c t o r 1/ (NT*DT*DT) so t h a t the time s e r i e s r eturned 355 c by the IFFT has the c o r r e c t magnitude. 356 c The c o o r d i n a t e s of each p o i n t , and the amplitude spectrum 357 c and time s e r i e s of e l e c t r i c p o t e n t i a l at each p o i n t are 358 c w r i t t e n t o a b i n a r y output f i l e . 359 c 360 RECNUM=3 361 DO 90 J=l,NNODES 362 MAXT(J)=0.D0 3 63 MAXF(J)=0.D0 364 ITM(J)=1 365 IFM(J)=1 366 RECNUM=RECNUM+1 367 WRITE(8,REC=RECNUM) SNGL(X(J)) 368 RECNUM=RECNUM+1 369 WRITE(8,REC=RECNUM) SNGL(Y(J)) 370 RECNUM=RECNUM+1 371 WRITE(8,REC=RECNUM) SNGL(Z(J)) 372 DO 85 K=1,NF 373 85 POTF(K)=CSCALE*CPOT(J, K) 374 CALL DFOUR2(POT,NT,NODIM,1,-1) 375 DO 87 K=1,NF 376 TEMP=CDABS(CPOT(J,K) ) 377 IF (TEMP .GT. MAXF(J)) THEN 378 MAXF(J)=TEMP 379 IFM(J)=K 380 END IF 381 RECNUM=RECNUM+1 382 WRITE(8,REC=RECNUM) SNGL(TEMP) 383 87 CONTINUE 384 DO 88 K=l,NT 385 IF (DABS(POT(K)) .GT. DABS(MAXT(J))) THEN 386 MAXT(J)=POT(K) 387 ITM(J)=K 388 END IF 389 RECNUM=RECNUM+1 390 WRITE(8,REC=RECNUM) SNGL(POT(K)) 128 391 88 CONTINUE 392 90 CONTINUE 393 C 394 WRITE(4,940) 395 WRITE(4,941) 396 WRITE(4,942) 397 940 FORMAT(/' The c o o r d i n a t e s and maximum e l e c . p o t e n t i a l ' ) 398 941 FORMAT(' i n both the time and frequency domains are') 399 942 FORMAT(' l i s t e d below:') 400 WRITE(4,945) 401 945 FORMAT(1 PT. ',2X, •X',6X, 'Y',6X, 1Z',5X, 'MAX TDOM. •, 402 + 4X,'TIME',3X,'MAX FDOM.',5X,'FREQ.') 4 03 DO 100 J=l,NNODES 404 100 WRITE(4,950)J,X(J),Y(J),Z(J),MAXT(J) , 405 + DBLE(ITM(J)-1)*DT*1000.DO,MAXF(J),F(IFM(J)) 406 950 FORMAT(I3,3F7.2,1X,D11.4,F9.3,D11.4,F10.2) 407 C 408 C C a l l a s u b r o u t i n e t o w r i t e out i n f o r m a t i o n d e s c r i b i n g the 409 C o b s e r v a t i o n g r i d and the magnitude and phase of the 410 C e l e c t r i c p o t e n t i a l at every p o i n t i n the g r i d . The output 411 C f i l e i s w r i t t e n i n a format that can be read by the 412 C c o n t o u r i n g program KCONTR. 413 C 414 C CALL CONOUT(POT,X,Y,Z,NNODES,NNODEX,NNODEY,NS,XP,YP,ZP,SPHRAD ) " 415 C 416 CLOSE(UNIT=8) 417 WRITE(6,980) 418 980 FORMAT(• S i m u l a t i o n complete'//) 419 STOP 420 END 421 C 422 C 423 c ********************************************************** 424 C Subroutine PRESSUR - T h i s subroutine c a l c u l a t e s the time and 425 C frequency domain r e p r e s e n t a t i o n s of the p r e s s u r e f u n c t i o n 42 6 C t h a t i s a p p l i e d t o the i n t e r i o r of a s p h e r i c a l c a v i t y . The 427 C time domain e x p r e s s i o n i s 428 C p ( t) = b * H ( t ) * ( t " 2 ) * e x p ( - a * t ) 429 C 430 C where H(t) = 1, t>=0 431 C = 0, t< 0 432 C a = 2*WN i s the s c a l a r decay constant 433 C b = PMAX*((a*exp(1))"2)/4 i s a s c a l a r n o r m a l i z a t i o n 434 C f a c t o r t h a t ensures the maximum pressure i s PMAX. 435 C 4 36 C The F o u r i e r t r a n s f o r m of the pressure f u n c t i o n i s gi v e n by 437 C P(w) = 2*b/(a+iw)"3 438 C where w i s r a d i a n frequency 439 C i i s the square root of (-1). 440 C 4 41 C ************************************************************** 442 SUBROUTINE PRESSUR(PRT,PRW,W,PMAX,WN,NT,NF,DT) 443 INTEGER I,NT,NF 444 REAL* 8 PRT (NT) , W (NF) , WN, A, B, DT, T, PMAX 445 COMPLEX*16 PRW(NF),CDENOM,CTEMP 446 A = 2.D0*WN 447 B = (PMAX/4.DO)*(A*DEXP(1.0D0))**2 448 DO 10 1=1,NT 449 T=DBLE(1-1)*DT 450 PRT(I) = B*(T**2)*DEXP(-A*T) 451 10 CONTINUE 452 DO 20 1=1,NF 453 CTEMP = DCMPLX(A,W(I)) 454 CDENOM = CTEMP*CTEMP*CTEMP 455 PRW(I) = DCMPLX(2.D0*B,0.DO)/CDENOM 456 20 CONTINUE 4 57 RETURN 458 END 459 C 460 C ************************************************************** 461 C Subroutine STRESS - C a l c u l a t e s the s t r e s s t ensor at p o i n t 462 C ( X ( l ) , X ( 2 ) , X ( 3 ) ) as a f u n c t i o n of r a d i a n frequency W. 4 g3 Q ************************************************************** 4 64 SUBROUTINE STRESS(CSIGMA,UF,W,XI,X2,X3,R,MU,V,NF) 465 INTEGER I,J,NF 466 REAL*8 W(NF),X(3),XI,X2,X3,R,MU,V,TI,T2 467 COMPLEX*16 CSIGMA(6,NF),UF(NF),C4,C5,C6,COM 468 X(1)=X1 469 X(2)=X2 470 X(3)=X3 471 DO 20 1 = 1, 3 472 TI = (3.DO - 4.DO*(X(I)/R)**2)/R 473 T2= -(1.D0 + 2.DO*(X(I)/R)**2)/V 474 DO 10 J=1,NF 475 CSIGMA(I,J) = DCMPLX(MO,0.DO)*DF(J)*DCMPLX<T1, W (J)*T2) 476 10 CONTINUE 477 20 CONTINUE 478 TI = 2.D0/R 479 C4 = DCMPLX(-T1*MU*X(2)*X(3)/R,0.D0) 480 C5 = DCMPLX(-TI*MU*X(1)*X(3)/R,0.D0) 481 C6 = DCMPLX(-T1*MU*X(1)*X(2)/R,0.DO) 482 DO 30 J=1,NF 483 COM = UF(J)*DCMPLX(T1,W(J)/V). 484 CSIGMA(4,J) = C4*COM 485 CSIGMA(5,J) = C5*COM 486 CSIGMA(6,J) = C6*COM 487 30 CONTINUE 488 RETURN 489 . END 4 90 C 491 C 492. Q ************************************************* 4 93 C Subroutine GEOMPZ - T h i s subroutine c a l c u l a t e s the c o o r d i n a t e s 494 C of a l l spheres making up the p i e z o e l e c t r i c body. The body may 495 C be a s i n g l e sphere or a v e r t i c a l l a y e r of spheres 496 C l y i n g i n the plane X=XPLANE. The l a y e r i s r e s t r i c t e d t o a t h i c l 4 97 C ness of one sphere but i t s l a t e r a l dimensions are s p e c i f i e d by 498 C the user. 499 C 500 c ************************************************************** 501 c 502 SUBROUTINE GEOMPZ(XP,YP,ZP,SPHRAD,NS) 503 INTEGER I,J,K,NS,NSY,NSZ 504 REAL*8 ZP(*),YP(*),XP(*),ZMIN,ZMAX,YMIN,YMAX,SPHRAD,SPHDIA 505 REAL*8 ZTEMP,XPLANE 506 READ(2,*) XPLANE 507 READ(2,*) YMIN,NSY 508 READ(2,*) ZMIN,NSZ 509 SPHDIA=SPHRAD*2.DO 510 K=0 511 DO 30 1=1,NSZ 512 ZTEMP=ZMIN+DBLE(I-1)*SPHDIA 513 IF (MOD(I+2,2) .EQ. 0) THEN 514 DO 10 J=1,NSY 515 K=K+1 516 XP(K)=XPLANE 517 YP(K)=YMIN+SPHDIA*DBLE(J-l)+SPHRAD 518 ZP(K)=ZTEMP ' 519 10 CONTINUE 520 ELSE 521 DO 20 J=1,NSY 522 K=K+1 523 XP (K)=XPLANE 524 YP(K)=YMIN+SPHDIA*DBLE(J-l) 525 ZP(K)=ZTEMP 526 20 CONTINUE 527 END IF 528 30 CONTINUE 529 NS=K 530 ZMAX=ZMIN+SPHDIA*DBLE(NSZ-1) 531 YMAX=YMIN+SPHDIA*DBLE(NSY-1)+SPHRAD 532 IF (NS .EQ. 1) THEN 533 WRITE(4,600) 534 WRITE(4,601) XP (1) , YP (1) , ZP (1) 535 600 FORMAT(/' The p i e z o e l e c t r i c body c o n s i s t s of one sphere') 536 601 FORMAT(1 c e n t r e d at x=',F8.3, ' y=',F8.3,' z=',F8.3, ' . ') 537 ELSE 538 WRITE(4,610) NS 539 WRITE(4,620) XPLANE 540 WRITE(4,625) YMIN,YMAX 541 WRITE(4,630) ZMIN,ZMAX 542 610 FORMAT(/' Number of spheres i n the p i e z o e l e c t r i c l a y e r : ' , I 543 620 FORMAT(' The l a y e r l i e s i n the v e r t i c a l plane x=',F8.3) 544 625 FORMAT(' and extends from y=',F8.3,' to',F8.3) 545 630 FORMAT(' and from z=',F8.3,' to',F8.3,'.'/) 546 END IF 130 54 7 RETURN 548 END 549 C 131 55 0 Q **************************************************************** * * * 551 C Subroutine GEOMOB - T h i s subroutine c a l c u l a t e s the c o o r d i n a t e s o f 552 C the p o i n t s at which the e l e c t r i c f i e l d i s t o 553 C be c a l c u l a t e d . 554 C A l l p o i n t s are l o c a t e d at the i n t e r s e c t i o n s of l i n e s i n a 555 C r e c t a n g u l a r g r i d l y i n g i n a h o r i z o n t a l plane. The user s p e c i f i e s 556 C the g r i d by su p p l y i n g the f o l l o w i n g parameters i n an input f i l e : 557 C 558 C 559 C 560 £ **************************************************** *** 561 SUBROUTINE GEOMOB(X,Y,Z,NNODES,NNODEX,NNODEY) 562 INTEGER I,K,L,KK,LL,NRMAX,NNODEX,NNODEY,NNODES,FLAG 5 63 PARAMETER (NRMAX=15) 564 INTEGER NSUBX(NRMAX),NSUBY(NRMAX) 565 REAL*8 X(*),Y ( * ) , Z 566 REAL*8 DX(NRMAX),DY(NRMAX),RX(NRMAX),RY(NRMAX),XTEMP 567 C 568 NNODEX=l 569 NNODEY=l 570 READ(3,*) Z 571 READ(3,*) NRX,NRY 572 READ(3,*) (RX(I),1=1,(NRX+1)) 573 READ(3,*) (NSUBX(I),1=1,NRX) 574 DO 10 1=1,NRX 575 DX(I)=(RX(I+1)-RX(I))/DBLE(NSUBX(I)) 57 6 NNODEX=NNODEX+NSUBX(I) 577 10 CONTINUE 578 READ(3,*) (RY(I),1=1,(NRY+1)) 579 READ(3,*) (NSUBY(I),1=1,NRY) 580 DO 20 1=1,NRY 581 DY(I)=(RY(I+1)-RY(I))/DBLE(NSUBY(I)) 582 NNODEY=NNODEY+NSUBY(I) 583 20 CONTINUE 584 C 585 C The nodes w i l l be numbered s e q u e n t i a l l y i n the y - d i r e c t i o n . 586 C 587 1=0 588 FLAG=0 '589 DO 50 L=1,NRX 590 DO 40 LL=1,NSUBX(L) 591 XTEMP=RX(L)+DBLE(LL-1)*DX(L) 592 31 1=1+1 593 X(I)=XTEMP 594 Y(I)=RY(1) 595 DO 34 K=l,NRY 596 DO 32 KK=1,NSUBY(K) 597 1=1+1 598 X(I)=XTEMP 599 Y(I)=RY(K)+DBLE(KK)*DY(K) 600 32 CONTINUE 601 34 CONTINUE 602 IF (L.EQ.NRX .AND. LL.EQ.NSUBX(NRX) .AND. FLAG.EQ.0) THEN 603 XTEMP=RX(NRX+1) 604 FLAG=1 605 GO TO 31 606 ' END IF 607 40 CONTINUE 608 50 CONTINUE 609 NNODES=I 610 WRITE(6,800) NNODES 611 800 FORMAT(' The p o t e n t i a l w i l l be c a l c u l a t e d at',16,' p o i n t s . ' ) 612 WRITE(4,810) NNODES 613 WRITE(4,811) 614 WRITE(4,812) Z 615 810 FORMAT(/• The e l e c t r i c p o t e n t i a l i s t o be c a l c u l a t e d at',16) 616 811 FORMAT(' p o i n t s a l l of which l i e i n the h o r i z o n t a l ' ) 617 812 FORMAT(' plane z=',F10.3,'.•) 618 RETURN 619 END 620 C 621 C g22 C **************************************************************** 132 623 C Subroutine CONOUT -624 C 625 c *************************************************** * * * .62 6 SUBROUTINE CONOUT(POT,X,Y, Z,NNODES,NNODEX,NNODEY,NS,XP,YP,ZP, R) 627 INTEGER TEMP(100),I,J,NNODES,NNODEX,NNODEY,START, NS 628 REAL*8' X (*) , Y (*) , Z, XP (*) , YP (*) , ZP (*> , R 629 COMPLEX*16 POT(*) 630 WRITE(7,*) NNODEX,NNODEY 631 WRITE(8,*) NNODEX,NNODEY 632 WRITE(7,700)(X(I),1=1,NNODES,NNODEY) 633 WRITE(7,700)(Y(I),1=1,NNODEY) 634 700 FORMAT(50F11.3) 635 DO 10 J=l,NNODEY 636 START=1+(J-l)*NNODEY 637 WRITE(7,710)(CDABS(POT(I)),I=J,NNODES,NNODEY) 638 10 CONTINUE 639 710 FORMAT(50D11.3) 640 C 641 C Determine which spheres l i e i n the same plane as the obser-642 C v a t i o n g r i d and w r i t e out the c o o r d i n a t e s of the c e n t r e s of 643 C those spheres so t h a t they may be p l o t t e d on the contour map. 644 C 645 J=0 646 DO 20 1 = 1,NS 647 IF (DABS(Z-ZP(I)) .LT. 0.01D0) THEN 648 J=J+1 649 TEMP(J)=I 650 END IF 651 20 CONTINUE 652 WRITE(7,*) J , 0 653 DO30 1=1,J 654 WRITE(7,700) XP(TEMP(I)),YP(TEMP(I)),R 655 30 CONTINUE 656 RETURN 657 END PZM.PRM 1 256 NT: no. of p o i n t s d e s i r e d i n time s e r i e s 2 0.000078125 DT: sample i n t e r v a l d e s i r e d f o r time s e r i e s 3 3000.0" V: s e i s m i c wave v e l o c i t y (m/s) 4 2.0D+10 MU: shear modulus 5 250.0D+06 PMAX: maximum pressure w i t h i n c a v i t y 6 1.500527 A: r a d i u s of c a v i t y 7 1000.0 RES: r e s i s t i v i t y of the whole space (ohm*m) 8 8.0 DIST: r a d i u s f o r w r i t i n g out d i s p l . and v e l . f u n c t i o n s 9 0.5 SPHRAD: r a d i u s of the p i e z o e l e c t r i c spheres (m) 10 5 NNZERO: number of non-zero p i e z o e l e c t r i c moduli 11 1,6,2.3D-12 I , J , D ( I , J ) : i , j and the value of modulus d ( i j ) 12 3,4,2.3D-12 13 2,1,1.15D-12 14 2,3,1.15D-12 15 2,2,-2.3D-12 16 23.2 DMAX: maximum d i s t a n c e from shot t o a c t i v e p i e z o . sphere 17 0 ISPEC1=1 i f a l l sphere l o c a t i o n s are t o be given i n t h i s f i l e 18' 1 ISPEC2=1 i f a l l o b s e r v a t i o n p t s . are t o be given i n t h i s f i l e 19 76 NP 20 13,-10,0 21 13, -5,0 22 13, -2,0 23 13, -1,0 24 13, 0,0 25 13, 1,0 26 13, 2,0 27 13, 5,0 28 13, 10,0 29 13, 15,0 30 13, 20,0 31 18,-10,0 32 18, -5,0 33 18, -2,0 34 18, -1,0 35 18, 0,0 36 18, 1,0 37 18, 2,0 38 18, 5,0 39 18, 10,0 40 18, 15,0 41 18, '20, 0 42 23,-10,0 43 23, -5,0 44 23, -2,0 45 23, -1,0 46 23, 0,0 47 23, 1,0 48 23, 2,0 49 23, 5,0 50 23, 10,0 51 23, 15,0 52 23, 20,0 53 18, 0,-5 54 18, 0,-2 55 18, 0, 0 56 18, 0, 2 57 18, 0, 5 58 18, 0,10 59 18, 0,15 60 18, 0,20 61 18, 5,-5 62 18, 5,-2 63 18, 5, 0 64 18, 5, 2 65 18, 5, 5 66 18, 5,10 67 18, 5,15 68 18, 5,20 69 18, 10,-5 70 18, 10,-2 71 18, 10, 0 72 18, 10, 2 73 18, 10, 5 74 18, 10,10 75 18, 10,15 76 18, 10,20 77 13, 2,-5 78 13, 2, -2 79 13, 2, -1 80 13, 2, 0 81 13, 2, 1 82 13, •2, 2 83 13, 2, 5 84 13, 2, 10 85 13, 2, 15 86 13, 2, 20 87 13, 2, 0 88 18, 2, 0 89 23, 2, 0 90 28, 2, 0 91 38, 2, 0 92 48, 2, 0 93 58, 2, 0 94 78, 2, 0 95 108, 2, 0 # # # # # # #LIST SPHERES.PRM 1 8.0 XP, (REAL*8) 2 -21.75,44 YMIN,NSY, (REAL*8,INTEGER) 3 -21.50,44 ZMIN,NSZ, <REAL*8,INTEGER) 4 ************************** 5 EXPLANATION OF PARAMETERS: 6 XP: x - c o o r d i n a t e of the p i e z o e l e c t r i c body 7 YMIN,ZMIN: s m a l l e s t y and z c o o r d i n a t e s f o r any sphere i n the 8 p i e z o e l e c t r i c body 9 NSY,NSZ: number of spheres i n the y and z d i r e c t i o n s APPENDIX B Mathematical derivations (1) Derivation of the stress tensor associated with a spherical, compressional elastic wave (2) Derivation of the Fourier transform of the approximate particle step response given by equation (3.9). 136 D Derivation of the stress tensor associated with a spherical, compressional elastic wave Consider a compressional elastic wave expanding radially from the origin of a spherical coordinate system. Let the particle displacement of as a function of time and radius be T?(r,t) = u (r,t) r If a Cartesian coordinate system OxjX^ is centred at the same point as the spherical system V 2 2 2 Xj + X2 + X3 and the particle displacement at point r = (x1,x2,x3) may be expressed in terms of Cartesian unit vectors as follows u(Xj,x2,x3,t) r = Uj + u 2x 2 + u3 x3 (B1) u x, u x2 u x3 where u t = — , u 2 = — , u3 = — . In an isotropic elastic medium the stress tensor may be expressed in a Cartesian coordinate system as follows (Sheriff and Geldart, 1985, p.35) \ (B2) where X and p. are the elastic constants of the medium ei; = are normal strains " dx; 3u; du: Ejj = Ejj = ^  + ^  i^ j are the shearing strains and A = e u + + e33 is the dilatation. Differentiation of the components of particle displacement u,, u2, and u 3 gives the components of strain. For example, 137 Figure B-l . Spherical and Cartesian coordinate systems with a common origin. An elastic wave is assumed to expand radially away from the origin. 138 . du . . dr but r = (xt + x2 + x3) .-. | ^ = X ! (x,2 + x 2 2 + x 3 2 ) - m = X l r 1 1 9 l l j 2 3 =* en = ( x i r ) ^  + u (r" - x i r ) 1 xt a u x2 => e n = ( — + — — - - 7 ) u 1 1 v r r 8xj r3 ' The expressions for the other strain components are found in a similar way. The general expressions for these components are 1 ,1 3 u x ^ eH = - (1 + X: - - T ) 11 r v 1 3x; r 2 ' _\ du du 2 Xj Xj u  Eu ~ e Ji - r ( x ' dxi + X i dx/ ' r3 The eSj are combined according to equations (B2) to yield the following expressions for the components of seismic stress at point (x!,x2,x3) located a distance r from the origin (Xl,x2,x3,t)=-rjjj (Xl,x2,x3,t) = r 2 ^ + 2 p du X-21 r 2 . + (k + 2 p ) Xj | ^ + X du du x ^ X j + X kax k , i * j * k X i a X j + x j a X j du 2 X ; X : U ID Derivation of the Fourier transform of the approximate particle displacement step  response given by equation (3.91: a2 -0)o(t-to)A/2 2>/2UJ us(r,t) = H(t-g-7=— e sinco0(t-t0) (3.9) , 0 t<0 where H(t) = 1 l t > Q and a, 0)o and t. are constants. a2 ^ A = 2V2lir * Define the Fourier transform of a function g(t) as oo G(co) = J{g(t)} = fg(t)ej(0tdt , —oo and define the inverse transform as oo g(t)= FHG(G>)} =^  /gtV^dco . —oo Then, the Fourier transform of (3.9) may be determined as follows: J{us(r,t)} = A J{H(t-t0) e'^'1^2 sinco0(t-t0)} . -j<at0 Using the time shifting property J {g(t-t0)} =G(co)e we find -/cotn -contA/2 J{us(r,t)}= AeJ ° J{H(t)e 0 sinoy} . (B3) Express the right hand side of (B3) as a convolution in the frequency domain rather than multiplication in the time domain as follows: -/cot. 1 -co^ tA/2 5{us(r,t)}=Ae-/ °[ — J{U(x)e 0 } ® J{sincoct}] (B4) where <8> denotes convolution. The two Fourier transforms on the right hand side of (B4) are relatively easy to evaluate: 140 7 {sinco0t} = - ; TE [8 (co - co0) - 8 (co + coj] (B5) where 8 (co) is the Dirac delta function, and 7{ H(t) = / H ^ 1 ^ dt OO = ± - r P<<°J& + » t I (a>0A/2 + ja>)16 J 0 -co„tA/2 1 J { H ( t ) e 0 } = , , c; . > (B6) (co0/V2+;co) The convolution of two functions F(co) and G (co) is O O F(co) ® G(co) = / F (x) G(co - x) dx . — O O Therefore, substituting (B5) and (B6) into (B4) and expressing the convulution in its integral form yields , , -./Ae" ;mt° 7 8(x-03o) - 8(x + co0) J{u.(r,t)} = J —x I , / - , . , — dx . s 2 -> (o0A/2 + ;(co - x) — O O Multiplying both sides of this equation by j/j we obtain , / Ae" ; C O t° 7 8(x-co0) - 8(x+ (0o) ^ ( u s ( r , t ) } = ^ ; ; W V - 2 ) . ( r o . T ) dx 5{us(r,t)} = 2 [ . K A / - 2 } . ( ( 0 . - ; W V - 2 ) . ( c o + roo )] (B7) Expressing the right hand side of (B7) in terms of a common denominator and simplifying the resulting expression we find 141 J{us(r,t)} =Ae ox. (/co + co0A/2 ) 2 + co0 2 (B8) Let Us(r,co) = J{us(r,t)} . Recall that A = Substituting for A in equation (B8), the following expression is obtained for the Fourier transform of the approximate particle displacement step response: Us(r,co) = 2>/2 p.: 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0052923/manifest

Comment

Related Items