MULTIPLE DYNAMIC MATCHING AND MEDIAN A P P L I E D T O SONIC W E L L L O G S FILTERS by W. S C O T T P O W E L L LEANEY B.Se.(Geophysics), T h e University of Manitoba, 1983 A T H E S I S S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F THE REQUIREMENTS FOR T H E DEGREE OF M A S T E R OF SCIENCE in T H E F A C U L T Y O F G R A D U A T E STUDIES Department of Geophysics and Astronomy We accept this thesis as conforming to the required standard T H E UNIVERSITY O F BRITISH November COLUMBIA 1986 © W a l t e r Scott Powell Leaney, 1986 In presenting degree freely at this the available copying of department publication this or of thesis in partial fulfilment University of British Columbia, for and reference thesis by this for his thesis scholarly or for her Department V6T DE-6(3/81) 1Y3 purposes Columbia the requirements I agree I further that agree may be It is representatives. financial gain permission. The University of British 1956 Main Mall Vancouver, Canada study. of shall not that the permission granted an advanced Library shall by understood be for allowed the for make extensive head that without it of copying my my or written 11 ABSTRACT Nonlinear signal matching is a generalization of cross correlation in that a discrete lag between signals is replaced with a variable lag function or 'matching function', m(x). T w o methods are reviewed which attempt to solve for m(i), namely the dynamic programming approach and the inverse theory approach. Both methods suffer from pitfalls and require the input of prior constraints to ensure convergence to the correct solution. The goal of this work has been to develop a method that can handle simple or complex matching problems and can succeed without any prior knowledge constraints. T h e multiple dynamic matching method is the result. It uses a significance threshold to extract a set of ridge points from a similarity matrix and applies dynamic programming to obtain a set of matched sections. These significant matched sections or 'subpaths' are then combined into a set of complete matching functions and a 'mean local confidence' norm is evaluated to determine the best overall match. It is argued through a model of change between signals, that given that the correct matching function is known, the presence of large amplitude changes between signals can cause the correct matching function to appear suboptimal under similarity norms. Multiple dynamic matching, because it generates suboptimal solutions as well, does not overlook the correct matching function. Typically the top three interpretations as ranked by mean local confidence will contain the expert's choice for the correct matching function. The use of median filters to preprocess for matching has been investigated. the data and enhance well log features A new 'median decomposition' is discussed as well, and in the context of a scale - space point of view to filtering, it is argued that median scale space is the proper choice for blocky waveforms. Finally, the connection between multiple dynamic matching and pattern recognition is explored, and matching iii iteratively through scale is proposed as a means of further generalizing the multiple dynamic matching method, making efficient high resolution matching possible. iv TABLE OF CONTENTS ABSTRACT ii LIST O F T A B L E S vi LIST O F F I G U R E S vii ACKNOWLEDGMENTS ix CHAPTER I INTRODUCTION . 1 1.1 Correlation and Matching Functions 1 1.2 Waveform Change and the Matching Problem 3 1.3 Dynamic Waveform Matching 4 1.4 T h e Inverse A p p r o a c h to Signal Correlation 8 1.5 Multiple Dynamic M a t c h i n g 11 1.6 T h e Well L o g Matching Problem 13 C H A P T E R II P R O C E S S I N G SONIC W E L L L O G D A T A W I T H M E D I A N F I L T E R S . . . . 16 2.1 Introduction to M e d i a n Filtering 16 2.2 Median Filtering Applied to Well Logs 18 2.3 C o m p o u n d Median Filtering 19 2.4 Median Decomposition 23 2.5 T h e ' B l o e t r u m ' 24 2.6 Median Scale Space 27 2.7 Automatic Event Detection 29 2.8 Processing Sonic Well Logs For Matching 31 C H A P T E R III MULTIPLE DYNAMIC MATCHING 38 3.1 Introduction 38 3.2 T h e Similarity Matrix 38 3.2.1 The Semblance Functional 39 3.2.2 Smoothness and Sampling 41 3.3 Local Confidence 44 3.4 Ridge Points 47 3.5 Subpath Determination (Dynamic Matching) 48 3.6 T h e Function Network (Multiple Dynamic Matching) 52 3.7 Ranking the Functions 54 3.7.1 Mean Local Confidence 54 3.7.2 Simplicity 57 3.7.3 Sufficient Diversity 57 3.8 A p p l y i n g the Matching Function 60 V 3.9 Computing T i m e 62 3.10 Sources of Error 63 3.11 Multiple Dynamic M a t c h i n g and Pattern Recognition 65 3.12 Automatic Stratigraphic Correlation 66 C H A P T E R IV CONCLUSIONS 4.1 Improvements 69 69 4.2 Parameter Estimation 70 4.3 Speculations on the Future and Related Ideas 72 4.3.1 Extending Well L o g M a t c h i n g to 3D 72 4.3.2 Seismic Signal Matching 73 4.4 Summary 74 Bibliography 75 Appendix 78 vi L I S T O F T A B L E S 1. F i v e classes of geologic change and digital reponse change 12 2. Standard errors for different niters 20 vii LIST OF F I G U R E S 1.1 Global and local path constraints 5 1.2 M a t c h i n g of envelopes of local earthquake data 6 1.3 Inverse theory solutions to a matching problem 10 1.4 Examples of matching functions for four classes of structural change 2.1 Simple examples of median filtering 2.2 Sonic log convergence to root for an N — 12 median 2.3 Comparison of median, compound median and Walsh . . 14 16 filter filtering 19 21 2.4 Original log and its median components 24 2.5 Residual components and original reconstructed. 25 2.6 The Bloetrum 26 2.7 Paths of zero crossings through scale 28 2.8 Automatic event detection 30 2.9 Trend removal using frequency low-cut 32 2.10 Trend removal using median low-cut 35 2.11 Median low-pass and frequancy low-pass 35 2.12 Preprocessing steps 36 2.13 Five logs after preprocessing steps 37 3.1 Computation of the similarity matrix 39 3.2 Logs 1 and 2 after preprocessing (example m /i 41 3.3 2D power spectra of similarity matrices 42 3.4 Local confidence illustration 45 2 3.5 Ridge points in similarity matrix for m / i 48 3.6 Local path constraint 49 2 viii 3.7 Subpaths from ridge points 50 3.8 T h e function network from subpaths 52 3.9 T o p ten functions after M L C 54 3.10 Application of best three functions to test log 55 3.11 57 Top ten functions after S I M P 3.12 T o p three functions after M L C and S I M P 58 3.13 Top three functions applied after sufficient diversity 60 3.14 Local confidence vs. x 64 3.15 Automatic correlation of five logs 67 4.1 Areal distribution of the five logs 72 IX ACKNOWLEDGEMENTS T h e problem of nonlinear signal correlation and the dynamic programming approach to solve it was first communicated to me by Toshifumi Matsuoka and I am very thankful to him for that. Throughout this work K e n Whittall has made countless fruitful suggestions and has helped to keep my goal in perspective. Thanks K e n . M y supervisor D r . T a d Ulrych has been a great source of inspiration and has enabled me to laugh at myself when things seemed impossible. T h a n k s to him for his special brand of encouragement. I would also like to thank the entire department of Geophysics and Astronomy and all my friends there. Special thanks also to my parents for supporting my education. 1 Chapter 1 INTRODUCTION 1.1 C o r r e l a t i o n a n d M a t c h i n g Functions The general problem of correlating two or more signals occurs in many fields of study. Examples of correlation problems come from geology (stratigraphic correlation), geophysics (seismic signal correlation), and speech processing (word recognition). A n y time that sequences of data have been recorded in an experiment and a connection between signals is suspected, a correlation problem exists. Obtaining optimum signal alignment requires determining the samples of the measured response which correlate or match in two or more sequences of data. Conventionally this has been solved using a correlation or similarity functional where the signals or portions of the signals are lagged relative to one another. T h e unnormalized cross correlation or similarity between a reference signal r(x) ; x — 0,X S where and a test signal t(x) { = / Jo r(x)t{x + £)dz, ; x = 0, Y may be written as (1.1) is the global similarity at lag £, (and the limits of integration may be changed to a window (Xi,X ) 2 for a local similarity measure). Finding the optimum lag simply 2 involves computing over a range of £'s and determining the m a x i m u m . However, some data such as well log and seismic data may undergo variable moveout so that a continuously varying lag is necessary to recover optimum signal alignment. T h u s the lag £ is replaced with a continuous lag function or 'matching function', m(x). The generalized similarity functional may now be written as (1.2) o where m(x) determines the shift, stretches and compressions of the test signal. One may want to determine the 'matching function' for a variety of reasons. connection of equivalent events may be wanted for an automatic interpretation The (for example connecting equivalent time horizons in well logs), or one may want to actually apply the matching function to align the data before additional processing. Another application of matching functions is to assign a waveform or portion of it to one of several classes as defined by their expected waveforms. Here one is looking for the simplest matching function to determine, for example, whether a seismic signal is due to an earthquake or nuclear blast. The application explored in this project is that of connecting structurally equivalent portions of digital sonic well logs. A l l applications of matching functions require an accurate and reliable method of determining the function itself, and if objectivity is required so that no a priori solution, the problem is particularily difficult. knowledge is included in the 3 1.2 W a v e f o r m Change and the M a t c h i n g P r o b l e m A simple mathematical model of change between waveforms is helpful to further explain the problem of correlation or matching. T h e change between a reference signal r(x) and a test signal t(x) may be modelled as due to a structural change m(x) amplitude distortion a(x), and an so that: r(x) It is desired to solve for m(x) = t[m(x)\ + a(x). and a(x) (1.3) and thereby quantify completely the change between signals. The trouble is that a(x) can only be determined once m(x) is known, and as m(x) cannot be separated from the test signal, its determination is a nonlinear problem. Also, the matching function may be simple or complex, smooth or discon- tinuous, and there is another complication as. well. amplitude distortion a(x) For most matching problems the is small so that nearly identical matching is possible if only m(x) is recovered. However, the presence of |a(x)| >> 0 for any x can cause the result of applying the correct matching function to appear suboptimal, so that any method that solves for m(x) by maximizing the coherence between r(x) mislead, determining an incorrect structural change. and <[m(x)] would be In these severe cases it would be beneficial to generate suboptimal solutions as well. The determination of the matching function is thus a complicated issue. The problem is nonlinear, the solution may be arbitrarily nonlinear, and amplitude distortions may cause the correct solution to appear suboptimal. One method that has been successful in the solution of nonlinear problems is that of dynamic programming (Bellman and Dreyfus 1962). However, only recently has dynamic programming been applied to 4 matching problems in speech processing (Rabiner, Rosenberg, and Levinson 1978, M y ers and Rabiner 1981), stratigraphic correlation (Gordon and Reyment 1979, Smith and Waterman 1980, Clark 1985) and seismology (Liu and Fu 1982). Dynamic program- ming was reviewed as a tool for geophysical data processing by Anderson and G a b y (1983). They discuss the method of dynamic waveform matching and its geophysical applications, and discuss some of the difficulties of the method with certain matching problems such as well to well correlation. Another method which is totally distinct from dynamic matching is the inverse approach to signal correlation (Martinson, Menke and Stoffa 1982). Both methods of waveform matching and their pitfalls will be reviewed in the next two sections, followed by an introduction to multiple dynamic matching and the well log matching problem in particular. 1.3 D y n a m i c W a v e f o r m M a t c h i n g Nonlinear matching shifts and warps (stretches and compresses) one of the signal axes (time or depth) until the corresponding features of the two waveforms are aligned. T h e optimum shift and warping of the signal axis, the matching function, may be determined by dynamic programming. T h e following is a brief review of this approach to waveform matching based largely on a paper by Anderson and Gaby (1983). Since it is desired to find the shift and warping that will minimize the difference between a reference waveform difference or distance ; i — l,n and a test waveform tj ; j — l , m , a measure may be computed between pairs of points to form a dissimilarity matrix, d ^ : Ti TO i=l j = l 5 T h e n the path through aV, which minimizes n 0 = (1-5) t=l is defined as the optimum path and its coordinates i,p[i) determine the shift and warp- ing of the test waveform axis for optimum waveform alignment. T h e matching problem is therefore viewed as a pathfinding problem through the matrix dij. Dynamic program- ming can efficiently determine the optimum path whenever the dynamic programming principle applies: Whenever the path from a starting choice of paths for point S to an intermediate influence the optimum travelling from minimum distance from S to G is the sum of the minimum minimum distance from point I does not S to a goal point G, then the from S to I and the distance I to G. (Anderson and G a b y (1983)). Dynamic programming is efficient because it assumes that the optimum local path is part of the optimum global path. In other words, only a small portion of dij local to the present path point needs to be computed to determine the next path point. To minimize equation (1.5) with dynamic programming several things must be specified: a) endpoint constraints - where the path may begin and end. b) local continuity constraints - constraints on the path from one point to the next. c) global path constraints - limitations on where the path can wander in (i,j) A graphical representation of contraints a and c is shown in Figure 1.1a. space. In most applications of dynamic waveform matching, the starting and end points are equivalent to knowing and constraining the correlation fixed, of points near the beginning 6 j(k) - i(k)+R F i g u r e 1.1 (a) E n d p o i n t and global path constraints, (b) type Hd local path constraint, (from Anderson and G a b y (1983)). E n d points (1,1) and ( N , M ) correspond to already interpreted matches between signals, and the sloped lines beginning from them are the result of the local path constraint in (b) (pathfinding may proceed from right to left or from left to right). Potential path points are commonly weighted as in (b), further restricting the simplicity of the path. and end of the test and reference waveforms. Local constraints on the next possible path point can take the form of bounds on proximity, slope, or slope weighting, wherein a potential path point dij may be weighted proportional to the resulting deviation of the local slope from one. shown in Figure 1.1b. T h e commonly used type H d local path constraint is The path constraints are equivalent to restricting the overall simplicity of the function, and along with the dynamic programming principle, make dynamic waveform matching an efficient method for solving for the matching function. However, dynamic matching often fails to recover the correct match when waveforms 7 1.4 1.0 0.6 0.2 5 30 TIME (s) F i g u r e 1.2 T h e matching of envelopes of local earthquake data using equations (1.4) and (1.5), and the path constraints of Figure 1. Dashed lines indicate matches. (From Anderson and Gaby (1983)). Notice the poor correlation of major events. T h e local path constraints and measure of difference are too restrictive to recover the correct match. show high variability in moveout and signature. This is because sections of waveforms which have no events of significant magnitude will produce an unclear match and poorly defined m i n i m a in dij, and this can mislead a pathfinding algorithm that looks ahead only locally. Also, pronounced moveouts or disappearances of significant events can require the pathfinding algorithm to look farther ahead than the prescribed local path constraint. Figure 1.2 shows the result of applying dynamic waveform matching using equations (1.4) and (1-5), and the path constraints of Figures 1.1a and 1.1b to match the envelopes of local earthquake data. Notice the poor correlation of primary and secondary events. This is due to an inadequate measure of dissimilarity (or similarity) and local path constraints which are too restrictive. T h e problems with the dynamic matching approach have been overcome by using more than one metric to compare the 8 signals (Liu and F u 1982), and by the inclusion of user supplied constraints (Gordon and Reyment 1979, Clarke 1985), where the user may constrain the path to include points representing already interpreted matches (the correlation of 'marker' events for instance). Obviously any method that requires user interpretation is not objective or amenable to the processing of large volumes of data. What is needed is a method that is not susceptible to being mislead by insignificant or unclear matches, can handle large and variable moveout, and can succeed without any user supplied interpretation. 1.4 The Inverse Approach to Signal Correlation The inverse theory approach to signal correlation avoids some of the pitfalls of the dynamic programming approach and can yield stable solutions to relatively simple matching problems. Martinson, Menke and Stoffa (1982) suppose that a reference signal r(t) is related to its distorted version d(x) through the matching function x(t). Thus r(t) — d'[x(t)]. They then parameterize x(t) as the sum of a linear trend and a truncated Fourier series: n-1 x(t) — a t + ^ 0 ai sin(t7rt/r) (1-6) t=i where T is the length of r[t). It is desired to maximize the coherence between r(t) and d[x(t)] which is given as: c = Io' r(t)d[x(t))dt fiT*{t)dtfZ'd*\z(t))dt W h e n the coherence is maximized the gradient of C with respect to the parameterization coefficients ai is zero: — oai = V C = 0. (1.8) A complicated expression for the partial deivatives may be derived from (1.7). The numerical values for the coefficients are then found by starting w i t h an initial guess for the solution a° and iterating to improve the coherence by a small amount A O . Using Taylor's theorem A C is expanded as: AC - dC dC dC - — A a + — A a ! + - + ——Aa„ = V C • A a . 0 , (1.9) % T h i s equation and n — 1 other constraint equations which come from requiring that the increment A a always be in the direction given by V C give a system of equations which when solved yield the change to the coefficients A a . Iteration stops when O ceases to increase. The method is applied to O differential sedimentation rates. 1 8 profiles from two ocean drill cores to determine A s with all nonlinear inverse methods, convergence to the correct solution is strongly dependent on the initial guess or starting function. M a r t i n s o n et al. obtain different results given a hand picked starting function and a simple linear starting function (see Figure 1.3). T w o points should be made w i t h regard to the inverse method of signal correlation. F i r s t , since global coherency is m a x i m i z e d and perturbations exist over the entire length of the matching function, local weaknesses in similarity do not greatly affect the final result and the method is more stable than dynamic matching. Secondly however, in difficult problems where the correct matching function is quite nonlinear, the inverse approach has difficulty in converging to the correct solution without a user supplied F i g u r e 1.3 Matching function determination with the inverse approach to signal correlation applied to O 1 8 data from two ocean drill cores. Solutions using (a) a hand-picked starting function (triangles) and (b) a linear starting function. (From Martinson et al. (1982)). Notice the sensitivity of the final solution to the starting function. starting function. Therefore, as with dynamic waveform matching, the inverse approach when applied to the most general class of matching problems, lacks objectivity and is not suitable for automatic data matching. 1.5 M u l t i p l e D y n a m i c M a t c h i n g A n approach to overcoming the problems of the previous two methods just discussed may now be clarified. Difficult matching problems could be reliably solved with the inverse approach if a starting function that is simple or complex as necessary could be determined automatically. Similarity, if local low confidence problems of the dynamic matching method in the presence of insignificant or unclear matches could be overcome, and local path constraints relaxed to allow for more complex matching functions, it too would satisfy the criteria of objectivity and reliability. The presence of large amplitude 11 distortions (|a(x)| > > 0) between signals which may cause the correct matching function to appear suboptimal could also be handled by either method if more than one starting function was used in the inverse approach, and if additional suboptimal solutions were found with a reliable dynamic matching algorithm. T h e previous two methods differ essentially in their definitions of optimality. The inverse theory approach maximizes a coherency functional over the full length of the signals and so uses 'global' optimality. Dynamic waveform matching minimizes the immediate difference between points or windows of points in the signals, and so uses 'local' optimality. It can be said that local optimality gives efficiency while global optimality gives stability. superior results. Combining the benefits of both methods would likely give Experiments with both matching schemes (Leaney 1985, 1986) have indicated the following features for a more robust waveform matching method: 1. To avoid insignificant and unclear matches, associate a significance threshold with local optimality to extract sections of the matching function where confidence is high, and use linear approximations where confidence is low. 2. Quantify final optimality over the entire function with a 'global' norm. 3. To avoid missing the correct solution due to |a(x)| >> 0, find suboptimal solutions as well as the optimal solution. T h e multiple dynamic matching method contains these features, and when applied to the problem of digital well log correlation which is known to be a difficult matching problem (Robinson 1978, Gordon and Reyment 1979, Schwarzacher 1982, Anderson and G a b y 1983), has produced good results. T h e method itself will be explained in detail in Chapter 3, but prior to that, a discussion of some of the peculiarities of matching well logs is needed. 12 1.6 T h e Well L o g M a t c h i n g P r o b l e m The problem of matching digital well logs poses problems not encountered with other data sets. Well logs are generally separated by kilometers or tens of kilometers and structural geology can change appreciably in that distance, meaning that m(i) can be considerably nonlinear. As well, the signature within a given formation may change due to differences in porosity, density, matrix composition, fluid content etc., so that amplitude distortions a(x) are large. The difference between two well logs may be expected to show any one of, or combination of five classes of digital response change. The following table summarizes the geologic changes and the corresponding digital response changes. geologic change digital response change 1. dipping or faulting shift 2. thinning/thickening compression/stretching 3. pinchout disappearance 4. faulting/fold ing repetition/in verted repetition 5. compositional change signature change Table 1. Summary of geologic and digital response changes between two well logs. The matching function m(i) determines the first four geologic changes, that is, the structural change between wells, while compositional changes internal to a formation or amplitude distortions a(x), may only be estimated given the correct matching function. 13 Figure 1.4 shows examples of matching functions corresponding to the four classes of structural change. Conventional cross correlation can handle only class 1 structural change while nonlinear matching methods can handle classes 1 - 3 of structural change. Class 4 structural changes are much more difficult to recover, but it will be seen that a simple modification to multiple dynamic matching would allow for the recovery of this class of structural change. Another problem is associated with well logs. Since they are generally not smooth signals in that they cannot be accurately represented by a linear combination of continuous functions, conventional filtering techniques cannot be applied with an expectation of undistorted waveform features. Signal matching problems usually require some pre- processing of the data, and since well logs are commonly signals composed of blocks and spikes, an alternative filtering technique that can remove short impulses without smearing block edges is needed. filtering of sonic well logs. Median filtering has been.found appropriate for the The data used in this project are five offshore sonic well logs with a sampling interval of 3 meters with each value being the weighted average of 10 readings spaced .3 meters apart. The next chapter deals with the theory of median filtering and its application to sonic well logs. 14 F i g u r e 1.4 Examples of matching functions for four classes of structural change and their analogous geological situations summarized in Table 1. T h e reference signal is along the x axis and the test signal is along the y axis, (l) Shift (recoverable by conventional cross correlation), (2) compression or stretching, (3) appearance or disappearance, (4) inverted repetition and repetition. Nonlinear matching methods can recover classes 1 - 3 of structural change. Multiple dynamic matching can also handle class 4 structural change. 15 Chapter 2 P R O C E S S I N G SONIC W E L L L O G D A T A WITH M E D I A N FILTERS 2.1 Introduction to Median Filtering Median filters operate by replacing the value at the middle of a sequence of numbers with the middle value of the sequence after sorting. The median of a sequence of numbers may be defined mathematically as that value that is a minimal 'distance' from all other values (Claerbout and M u i r , 1973). T h a t is, the median value minimizes the sum of the absolute differences between it and the other points of the sequence and so operates as an ' L i ' estimate of the expected value of the sequence. This is known to be a better estimate than the mean when the sequence contains wildly erratic values. T h e median operator M is also data dependent and is a nonlinear operator as it does not conform to the definition of linearity: F ( a x + b) — o F ( x ) + - f i b ) . For example take x = and b = (2, 3,4),and a = 2. Then ( M ( a x + b) =• 8] / the median operator is nonlinear. (1,2,3) [ a M ( x ) + M ( b ) = 7j and so These are properties that define the median, but of particular interest is the median as a filter, moving or sliding along a sequence of numbers. 16 ' i ' i i i i i i i i i i i i I I I I I I I I I I I I - I I I I I I I I JL_I i ; : i i : 11 I I I I !• I I l Figure 2.1 (a) Three point median filter applied to a binary signal (from Stewart (1986)), (b) convergence of a binary signal to its N = 1 root (from F i t c h , Coyle and Gallagher (1984)). The operation of sort - replace - slide is presented in Figure 2.1a for a three point median filter length applied to a binary signal. Notice the spike rejection and 'block pass' result. Generally any impulse that is less than N+l points long where 2JV + 1 is the filter length, will be rejected. However, on complex waveforms where sequence values change at each point, the ideal transfer 'behavior' of passing blocks greater than N in length is not seen. In real data examples the same length filter must be applied successively and even then the ideal transfer behavior is not observed. Gallagher and Wise (1981) have performed a theoretical analysis of the properties of median filters. They derive necessary and sufficient conditions for a signal to be invariant under successive median filtering. These conditions state that the signal must be locally monotone through a median filter unchanged. to pass A signal that is invariant under further median filtering is called a root signal, and it is shown that a root signal may contain only 17 'constant neighborhoods' and 'edges'. These two possible structures of root signals are denned for a median filter of length 2N + 1 as: 1) A constant neighborhood is at least TV + 1 identically valued points such that the constant neighborhood and the edge together are monotone. 2) A n edge is a monotone region between two constant neighborhoods of different value and the connecting monotonic region cannot contain any constant neighborhood. T w o other important characteristics are noted here. For a root signal to contain both regions of increase and decrease, the points of increase and decrease must be separated by a constant neighborhood (at least TV + 1 points of constant value). of length T will converge to its root signal in at most (T — 2)/2 Also, a signal successive filterings. Figure 2.1b shows the convergence of a binary signal to its root for a median filter of length 3. 2.2 M e d i a n F i l t e r i n g A p p l i e d to W e l l L o g s The block and edge pass property of median filtering is a desired property for the filtering of sonic well logs. Well logs are composed of blocks and impulses, and it is sometimes desired to remove the impulses while preserving the blocks and edges. Bednar (1983) has investigated the editing of acoustic impedance data with a Markov algorithm, running means and median filters and concluded that median filters produced superior results. Figure 2.2 shows the convergence of a portion of a sonic well log to its root signal for a median filter length of 25. T h e root signal contains only constant neighborhoods or blocks of greater than 12 samples, and edges or monotonic regions between them. It is noted that an edge may be as little as one sample long. The degree to which median filtering does not preserve perfect steps in the original data, 18 F i g u r e 2.2 Convergence of a portion of a sonic well log for a median filter of length 25. 19 and produce only blocks, is dependent on the distribution of values near an edge. Evans (1982) points out that nearby spikes may be drawn into a step (making it a series of steps), or, if values are widely distributed near the step, the corner may be rounded. Median filtering when applied to sonic logs may not preserve block edges ideally due to the erratic nature of the data, but it is superior to other filtering techniques. 2.3 Compound Median Filtering Median filtering for the removal of short impulses and the preservation of blocks and edges can be seen in Figure 2.2. Notice however the creation of multiple steps and rounded corners on some blocks. A s mentioned previously, Evans (1982) attributes these artifacts to the distribution of values near the steps. It would seem that the removal of erratic values before the application of a median filter would help to reduce the multiple steps and rounded corners. A n innovation to median filtering not found in the literature is employed here to remove these artifacts and is called 'compound median filtering'. C o m p o u n d median filtering operates by removing all impulses of length < N before application of a 2N + 1 length median filter. This is accomplished by median filtering the input signal with a filter of length (2n + 1) to its n = 1 root signal, then filtering this result to its n = 2 root signal and so on until n — N. Gallagher and Wise (1981) state that the root of a signal depends on the median filter window length. For example, the root for a window of length 5 is always a root for a window of length 3, but the root for a window of length 3 is not necessarily a root for a length 5 window. C o m p o u n d median filtering uses this fact and would only apply a filter of length 5 to the root signal of a length 3 filter. The result is a different root signal for the same length filter, and since any impulse of length < J V has been removed, multiple 20 step and rounded corner artifacts should be reduced. Figure 2.3 shows the same original log portion as in Figure 2.2, the root signal for an N — 12 median filter and the root signal for the corresponding compound median filter. Also shown is the Walsh low-pass version of the original log. T h e use of Walsh transforms in the representation of well log d a t a has been explored by Lanning and Johnson (1983). Walsh transforms map the data onto an orthogonal set of square wave functions. Each function is described by its 'zero crossings per second' or 'zps' (also called 'sequency'), and reconstruction m i n i m u m block length is then possible by zeroing energy in sequencies to any corresponding to blocks below a desired length and inverse Walsh transforming. The Walsh low-pass version in Figure 2.3 has blocks of length N, the same as the median and compound median filtered logs. B y computing the root mean squared difference or standard error between each of the filtering results in Figure 2.3 and the original log, the best representation original log at the given block length may be determined. of the The table below shows the standard error (SE) for each result. Filter T a b l e 2. type SE median 191.1 compound median 174.8 Walsh 215.7 low-pass Standard error between the original log and median, compound median and Walsh low-pass versions of the original. 21 0 Figure i i i i i i i i t 2.3 i i 100 11 i i i 11 1 1 1 1 11 11 1 1 1 1 1 200 i i i 11 1 1 i i i {i i i 11 1 1 i i 300 111 1 1 1 11 i [ i i i i i i i 11 i (a) Original data (same as in Figure 2.2), 111 400 111 i 111 i 11 1 1 1 1 500 11 i 111 i i i 1111 (b) root signal for an N = 12 median filter, (c) the root signal for the corresponding compound median filter, and (d) the Walsh low-pass version of (a) with equivalent shortest block length as in (b) and (c). T h e compound median result of (c) best represents the signal at the given m i n i m u m block length. 22 This indicates that compound median filtering is the superior technique for extracting the blocky portion of a signal. The improved result of compound median filtering over conventional median filtering brings the behavior of a compound median filter close to the 'ideal transfer behavior' attributed to median filters and mentioned earlier. A practical application of compound median filtering will be discussed in section 2.7. 2.4 Median Decomposition The following is a new and original nonlinear signal decomposition which is founded on the innovation of compound median filtering discussed in section 2.3. The effect of the compounding process is to successively remove impulses of length n = 1 , 2 , . . . N. Impulses of length 2 for example are only removed after those of length 1 have been removed. Since compound filtering removes impulses of a fixed length, the difference between present and past root signals will be a signal containing impulses of only one length. T h e signal may thus be decomposed into its constituent 'residual component' signals, each containing impulses of one duration. If o{J))3 — \,m represents the original sonic log, and / „ ( • ) represents the operation v of a median filter of length 2n + 1, then the 'residual component' with only impulses of length n is: a (j) n = v -i[j) n - v (j) n where f (j) = /nK-i(j)). B (2.1,2.2) When this is done for n — 1,7V, the signal has been decomposed into a root signal and 23 a set of TV residual components a (j), n and the original signal may be reconstructed as: TV vo(j) = v {j) + J2 n{J)- (2-3) a N It is noted that this is not a linear decomposition since a (j) n is an arbitrary and unknown function. Figure 2.4 shows each compound median root or 'median component' for N — 43 and Figure 2.5 shows the residual components a (j) n and the original signal reconstructed from equation 2.3. It should be mentioned that the reconstruction is exact and does not suffer from round off error as no multiplications or divisions are inolved in the computation of the median or residual components. 2.5 The 'Bloctrum' Just as there is a frequency spectrum showing energy at different sine wave periods, so may the energy or amplitude at different impulse or block lengths be computed from a median decomposition. The amplitude A n of all impulses of length n may be computed as: m A„ = J ] |o (j)| n for n=l,2,...JV. (2.4) Since compound median filtering decomposes the signal into impulses or blocks, amplitude as a function of block length may be called a 'bloctrum'. T h e bloctrum of a trend-removed log is shown in Figure 2.6. T h e bloctrum provides quantitative insight into the block distribution within the signal and may be employed as a quality measure in matching. For example, after 24 25 Figure 2.5 T h e residual components a ( j ) corresponding to the median n in Figure 2.4, and the original signal reconstructed from equation (2.3). components 26 40 35 i 30 .. 25 20 .. 15 10 .. 5 o lllllllllllll.llllllif 11 0 F i g u r e 2.6 20 11111 40 ,i—1 60 ,—L 80 100 T h e bloctrum corresponding to the median decomposition of the original log in Figure 2.4 but with N = 98. A m p l i t u d e s have been normalized by the length of the signal. matching, the block distribution within the logs should be similar, so a coherency measure between bloctrums should be higher for a better match. Also, the bloctrum may be used to reconstruct the signal with only blocks of significant energy. A threshold may be applied to the bloctrum or to each residual component individually so that only those residual components with significant amplitude would be used in a reconstruction. T h e signal may also be reconstructed exactly from fewer than N median since some components have no contribution to the signal. components 27 2.6 Median Scale Space T h e operation of median filtering extracts signal from noise where the division between the two is at an exact impulse duration or block length. A longer median filter length passes longer blocks, and this may equivalently be seen as reducing resolution and producing a broader scale view of the log. T h e log may thus be viewed through scale by changing the length of the median filter. T h e idea of scale in the areas of computer vision and pattern recognition has received a fair amount of attention in recent years. T h e goal is to quantitatively describe waveforms in one and two dimensions through a set of useful descriptive primitives. The local extrema of the signal and its derivatives are useful descriptive primitives, but the problem of scale consistently emerges as a source of difficulty. The location of edges for example, will vary in number and position, depending on the value of the signal at two or more points and the extent of the neighborhood around the point of measurement, so that a different description is obtained for different scales. What scale is the best to describe any particular waveform? W i t k i n (1985) has come up with a multiscale approach to signal description. The idea is quite simple. He chooses as his parameter of scale the standard deviation (a) of a Gaussian convolved with the signal. A s o is decreased continuously, beginning at a coarser scale, existing extrema move along the signal axis and new extrema occasionally appear at singular points. When the zero crossings are found along the signal axis at a number of intervals of scale, lines may be drawn in the (x,a) plane or 'scale space' to follow this primitive descriptor through scale. A s scale is coarsened, zero crossings disappear in pairs as local extrema disappear. Figure 2.7 shows the paths of zero crossings in scale space for an image intensity profile. Notice that the traces of zero crossings form arched lines at broader scales, due to the smearing affect of the Gaussian 28 F i g u r e 2.7 T h e paths of zero crossings through scale space. The scale parameter <r increases upwards. Notice the arches in the lines at broader scales due to the smearing effect of the Gaussian operator, (from Witkin (1984)) Compare with the paths of block edges through 'median scale space' in Figure 2.4. operator. If this operator was convolved with a sonic log, the smearing would be more pronounced due to the abrupt j u m p s in the signal. A n edge in the original signal would be shown at an incorrect location at a broad scale in ( i , a) space. A better parameter of scale in the case of blocky waveforms would be the compound median filter length. T h e paths of edges in 'median scale space' should remain straight, owing to the step preservation property, and should disappear rapidly over scale as the median window length is increased. Figure 2.4 shows a trend - removed sonic log through median scale space. T h e edges disappear at the the same x coordinate as scale is coarsened, indicating that their location has been preserved through scale. T h e superiority of median scale space for the description of sonic logs is clear and adds support for the use of median filters on these signals. More will be said later on the connection between a scale space 29 point of view and pattern recognition, and in particular, matching through median scale space will be discussed. 2.7 A u t o m a t i c Event Detection The result of preserving sharp block edges provided by compound median filtering may be exploited in the form of an automatic event detector. Events may be defined as abrupt excursions from the preceeding trend, and since this also defines a block in a well log, block edges are events to be detected. The event detection algorithm proceeds as follows. After compound median filtering has been applied to the log, a differencing procedure may be applied to generate a psuedo-reflectivity sequence, with events indicated where magnitudes are large. Since some block edges may exist over more than one point, locations which show adjacent values of the same sign are merged to the location of the largest value. This merging procedure is done over the entire psuedo-reflectivity sequence so that all values become either zeros or singular spikes. A simple thresholding may then be applied to extract only the largest events. The decision of the length of compound median filter to use is somewhat arbitrary. It must be as short as possible while still eliminating thin bed 'noise'. N — 12 was found suitable for the logs used in this project. A filter with Figure 2.8 shows a trend- removed log before and after compound median filtering, its psuedo-reflectivity, and the location of major events after merging and applying a threshold of 200 m / s . Figure 2.8 (a) original, (b) N = 12 compound median root signal, (c) the psuedo- reflectivity and (d) the major events after merging and applying a threshold of 200.0. 31 2.8 Processing Sonic Logs For Matching The goal of preprocessing the sonic data in this project is to enhance features for matching. In the case of well logs the most obvious matchable features are the block edges or steps. To enhance features for matching, two operations must be performed on the data. First, the velocity or background gradient which varies from well to well must be removed as it can induce a d.c. bias to any similarity measure computed between portions of the logs. It was found that conventional high pass frequency domain filtering (actually the time domain analog) was unsatisfactory for trend removal. This is because narrow high amplitude blocks require low frequencies for their representation, and the removal of low frequencies induces sidelobes adjacent to these blocks (see Figure 2.9). T o truely estimate the background trend, the trend upon which blocks sit, narrower high amplitude blocks should not be included. Median filtering can be used to estimate the background trend without 'feeling' blocks of a specified length (< TV), so the removal of a sufficiently long median filter output solves the problem of sidelobes. To obtain an estimate of the velocity gradient or trend that does not include the high amplitude blocks, a 105 point median filter is applied to the logs. This is then removed from the logs to accomplish the high-pass operation. Figure 2.10 shows the filtering step and the removal of the background trend from a sonic log. The second preprocessing operation performed on the sonic logs is low-pass filtering. Block edges are commonly imbedded in thin bed 'noise' and this can produce a greater uncertainty in any coefficient of similarity that is computed between portions of the signals, hindering the matching of block edges. Thus, a low pass operation is needed. Figure 2.11 shows the application of a 29 point median filter and an equivalent length 32 0 200 400 600 11111111111111 F i g u r e 2.9 800 1000 1200 11111111111111111111111111111 (a) original sonic log, (b) 105 point tapered mean of (a), (c) trend-removed log by subtracting (b) from (a). Notice the negative sidelobes induced adjacent to the large amplitude block located at 1020 - 1030 samples. 33 running mean applied an equal number of times. B o t h results remove the 'thin bed noise' but notice the preservation of edges in the median filter result. T h e processing of the sonic logs before matching thus involves two steps. First the removal of a 105 point median filter output from the original data, and then the application of a 29 point median filter for the low pass operation. The result is data with blocks greater than or equal to approximately 15 samples or 45 metres in length. Figure 2.12 shows the preprocessing steps one sonic log and Figure 2.13 shows all five logs used for matching after pre-processing. Unfortunately, for the purposes of processing sonic logs for matching, compound median filtering was found unsuitable, perhaps because it is too discriminating. It seems that a smoother transition between blocks (constant neighborhoods) makes matching easier. If some smoothing after compound median filtering was performed the result might be suitable for matching, but this possibility was not explored in the interest of keeping the preprocessing as simple as possible. 34 200 0 t 111111111 400 111111111 600 800 111111111 i 11111111 11111 1000 1200 111111111 111111111 111 1111111 t F i g u r e 2.10 (a) original sonic log, (b) 105 point median filter of (a), (c) trend - removed log (a minus b). N o negative sidelobes are seen adjacent to high amplitude blocks. 35 200 0 Figure 2.11 (a) 400 600 800 1000 1200 The trend-removed sonic log of Figure 2.10c, (b) 29 point running mean filter applied the same number of times (3) as was the 29 point median filter in (c) for convergence to the root signal. median filter result. Notice the preservation of block edges in the 36 0 200 400 11 600 800 1000 1200 111111111111 A 11 F i g u r e 2.12 11 11 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 M 1 1 1 1 1 1 1 1 T h e preprocessing steps applied to one sonic log. A 105 point median filter output is subtracted from the original and then a 29 point median filter is applied to the trend removed log. T h e log now has blocks of length greater than 14 samples in duration. 37 400 0 ' I II III I I | II I I I I I[ 800 | IIII I I II I|I I 1200 I I I I I II I | I I I I I I I I I|I II I I I I I I | I I I I I I I I I 1600 | M I II I III I|I II I I | I I I I |II I I1 IIII I I I I | I II I I I I I I | I I I IIII jIIII)IIII| I I |III I II I I I|I F i g u r e 2 . 1 3 T h e five logs after trend removal and median low-pass preprocessing. T h e background trend has been removed with a 105 point median filter and low-pass median filtering has been applied with a filter of length 29. T h e logs now have blocks of length greater than 14 samples. 38 Chapter 3 MULTIPLE DYNAMIC MATCHING 3.1 Introduction The foundation of multiple dynamic matching is in dynamic waveform matching, reviewed in section 1.3. It has as its basis a matrix known as a dissimilarity or similarity matrix depending on whether the functional values that make up the matrix are to be locally m i n i m u m or maximum for optimality. matching is a similarity The matrix used in multiple dynamic matrix although this is not necessary for the method to succeed. A s with dynamic waveform matching, multiple dynamic matching or ' M D M ' views the correlation problem as a pathfinding problem through the similarity matrix, but rather than finding only one path, many possible paths are found. T h e connection to the inverse approach to signal correlation (section 1.4) is that path optimality will be decided using global norms. In the sections that follow the details of M D M will be explained, beginning with the fundamentals of the similarity matrix. 3.2 The Similarity Matrix Since it is desired to maximize the similarity between the reference signal r(x) \,X and the test signal t(x) ; x = l,Y ; x — everywhere along their extent, the signals may be partitioned and a local similarity measure evaluated between windows of each signal, 39 where the relative lag between windows is given an upper bound IL. This bound on the relative lag is a measure of the greatest structural difference expected between the two signals and corresponds to the bandwidth of the similarity matrix Sij, computed as: L Sij — £ L Sij = m —l Y2i ( » A i , t ( i + l ) 4 i > ^ ) s r n-l Y 2 ( U+i)^x,tjAx,K), s r 1=0 j=l where the sampling increment between functional values is A x , n — (X — K)/Ax m = (Y — K)/Ax+ (3.1) + 1 and 1 are the number of points in the x and y directions respectively, s is the particular similarity functional and K is the window length over which similarity is computed. A s indicated in equaton (3.1), functional values are computed for (x, y) pairs corresponding to the locations of the first point in each of the windows of the reference and test signals respectively, and computation proceeds up the diagonals of the matrix for each relative lag /. This is done to make computation of the matrix most efficient, as advantage can be taken of operations already done. Figure 3.1 illustrates the computation of the similarity matrix 3.2.1 The The Sij. T h e Semblance Functional similarity functional used for the computation of Sij in this project is semblance. semblance functional for two signals r(x) and t(x) at windows i and j and for a window length K is: 8 = SlUK+fc-i+'i+fc-i)' . 2 " Sfc=i(ifc-i + r + tj k-i) + ( 3 . 2 ) 40 F i g u r e 3.1 Computation of the similarity matrix. ilarity window length is K, The half bandwidth is L, the sim- the matrix sampling increment is A x , and computation proceeds up the diagonals to take advantage of operations already done. Semblance is a measure of the o u t p u t / i n p u t energy ratio and is related to energy normalized cross correlation and summation or stacking. Neidell and Taner (1971) conclude from tests on multichannel seismic data that semblance has greater discriminating power in the case of ambiguously coherent events. Semblance has the range (0, l) and not (—1, l) as does cross correlation, and so has no built in phase or anticorrelation discrimination ability. In fact this is desired since in the case of matching sonic logs there is no utility in considering anticorrelations. Semblance is also related to the F-statistic 41 (Douze and Laster 1979). Although their experiments were performed on synthetic examples, their work indicates that larger semblance values indicate greater confidence in there being coherent energy between signals. 3.2.2 S m o o t h n e s s a n d S a m p l i n g Since the matching problem will be solved by pathfinding through the similarity matrix, it must be sufficiently smooth to avoid errors. It is still desired to sample the matrix (every A x ) to reduce its size and the cost of the algorithm, but this must not be done at the expense of accuracy. The question then becomes what is the largest A x that can be used to provide greatest efficiency while still retaining acceptable accuracy. T h e answer depends on the low pass median filter length. Median filtering leaves blocks of greater than N in length (for a filter of length 2N + 1), so from sampling theory each block must be sampled at least once to avoid aliasing. This puts an upper bound on the matrix sampling increment A x , and means that A x < N. Now that A x can be fixed, the matrix must be computed such that it is sufficiently smooth. It turns out that the smoothness of the similarity matrix can be controlled by the window length K of the semblance functional. The exact transfer function of the semblance functional is difficult to quantify, but as it is known to behave similarily to a stack or mean, it can be said that longer windows produce smoother matrices. In an attempt to establish the best value of K given A x and N, a 2D Fourier analysis has been done for similarity matrices computed using different semblance window lengths. It may be conjectured that below some value of K, noise of significant energy might appear in the matrix, making it unreliable for pathfinding. The data from which the matrices were computed for the 2D analysis are logs 1 and 2 after preprocessing 42 0 I I I M I I 400 I I | I I I I I I I I I | I I I I 800 I I t I I | I I I I | I 1200 I I I I I I I I | I I I I I I I I I | I I I I I I I I I | I I I I I I I I I | I I I I I I I I I | I I I I I I I I I | I M ( I I I I I | I I I I I II I I | F i g u r e 3.2 Logs 1 and 2 after preprocessing. These two logs will be used in an example of multiple dynamic matching which will be referred to as m / i . 2 (see Figure 2.12). These logs will also be used as an example to illustrate the multi- ple dynamic matching method and will be referred to as example m / i , with log 2 as 2 the reference log. Figure 3.2 shows logs 1 and 2 after the preprocessing steps of trend removal and low-pass median filtering (N — 14). Figure 3.3 shows similarity matrices computed for example m / i with A x =10 and K — 20,60,100 and their corresponding 2 2D power spectra. Notice that significant energy remains below the Nyquist frequency for all examples so this cannot be used to determine the best semblance window length. T h e 'character' of the matrices is different however. It is known that shorter windows allow for the matching of finer scale features. Given the shortest block length in the sonic logs (N), one might want the semblance window to be of approximately the same length. However, too short a window leads to an unreliable estimate of similarity (Douze 43 20 40 60 80 - 0.25 0.00 0.25 F i g u r e S.3 The similarity matrices (left) and 2D power spectra (right) for (a) K = (b) K = 60, and (c) K — 100. 20, A l l matrices show energy dropping to a very low value well below the Nyquist frequency, so this can not be used as a criterion for choosing the best value for K. (The spatial Nyquist frequency is expressed in terms of cycles per sample and so equals 0.5) The similarity matrix in (b) shows the best character relative to a priori knowledge. 44 and Laster 1979), so longer window lengths are needed for reliable similarity estimates. It is also known that trends in similarity should appear at roughly 45 degrees in the matrix, and should be clearly defined. The matrices of Figure 3.3 may now be viewed in a different light. For K — 20, contours of similarity follow block edges closely, with higher values where blocks of similar amplitude coincide. T h i s matrix does not conform to the expected 'good character' matrix and would clearly be unreliable for the extraction of trends in similarity. T h e K — 60 matrix shows the good character of well defined trends in similarity at roughly 45 degrees. The K — 100 matrix also shows coherent trends, but since it is desired to keep K as small as possible while still preserving the coherence of similarity trends, K — 60 is a better choice. somewhat arbitrary. Given that A x w 2/3 The selection of K is thus N to satisfy the block sampling criterion, K « 6 A x has been found appropriate to ensure the smoothness and character of the similarity matrix. 3.3 L o c a l Confidence The semblance coefficient between windows of points in the logs will be large (close to 1) when portions of the logs have similar amplitude and shape. In the similarity matrix semblance coefficients are located side by side, defining a surface, and different structures on this surface correspond to different situations of matching. For example, if amplitudes are roughly constant over a portion of the logs and events are small, any coefficient of similarity will not change appreciably between any two windows of these 'quiet' portions of the logs. The result is a flat region or plateau in the similarity matrix and a poorly defined match. Low confidence should therefore be placed in matches between portions of the logs with events of insignificant magnitude. Conversely the concept of high local confidence may be explained by the following illustration. 45 Imagine windows of each log 'sitting' on blocks of large and similar amplitude. A s one window remains stationary and the other is slid forward over the other log, it passes over the block edge and the amplitude drops rapidly. Similarity values drop also due to the sudden large difference in amplitudes between the windows. This will occur if the stationary and sliding windows are interchanged, so that a peak will appear in the similarity matrix in both x and y directions at some coordinates (x, y) corresponding to the vicinity of the two block edges. If the windows are slid at constant increments A x along both logs (following a diagonal in the matrix), the peaks will flatten where a match becomes unclear, and will sharpen where a significant match can be seen. T h e oceurance of repeated significant matches appears in the similarity matrix as a well defined trend or ridge. Local confidence in a match is therefore related to the magnitude of local similarity S, at a point on the matching function ( i , m(x), that is to S (z)> Iim and to the quality of the ridge in the direction of m(x). To quantify the quality of a ridge, the second derivative may be computed perpendicular to the direction of m(x) as it measures the downward curvature. The local confidence (LC) norm is now defined as: L C = - x,m(x)S" 5 m ( l ) l (3.3) Figure 3.4 shows portions of the logs of example m / i and the corresponding similarity 2 matrix. Where event amplitudes are large coherent trends in similarity are seen, and where events are smaller trends are less well defined. E a c h well defined ridge represents a region where matching confidence is high, and may be viewed as a portion of the matching function, m(x). Notice that ridges can run roughly parallel to one another, indicating more than one possible significant match for any one event on a log. This is due to the periodicity in the blocks of the logs. 46 The best overall matching function will be the function for which local confidence remains high throughout its length. This is easily expressed in terms of the mean local confidence ( M L C ) norm: ! p P =o where P is the length of the matching function. This is a global norm as it assesses the optimality of the entire matching function. More will be said on the computation of M L C in section 3.5.1. 3.4 Ridge Points Since nothing is known a priori about the optimum m(x), scanning Sij along all possible paths to find the maximum M L C is impractical. To handle this problem, points are located for which the local confidence is known to be large. may be constrained The magnitude of Sij to be large by requiring that it exceeds a significance threshold. T h e negative of the second derivative may be constrained to be large by locating only points which are peaks perpendicular to m(x). local direction of m ( x ) , Nothing is known a priori about the but it can be assumed to be in one of four directions defined by the eight points surrounding a central point. A simple test is then implemented to ensure that the point appears as a peak in three out of four directions. These two criteria make the location of local confidence maxima or 'ridge points' efficient, as the matrix Sij is simply scanned to find the points for which the above criteria are satisfied. T h e following summarizes the test for ridge points. 47 F i g u r e 3.4 Portions of the preprocessed logs and the similarity matrix with contours of similarity labelled so that 1 is low and 3 is high. Notice that where event amplitudes are large, at the beginning of the traces, trends in the similarity matrix are well defined. Toward the end of the traces where events are smaller, trends in similarity are less coherent. 48 1. Sij > .65 (significance threshold). 2. 3 of the following are true (a) Sj > Sij^i () S i,3 > Sij > 5 j _ i j _ ! and Sij > Si+ij+i > Si_i > Si ij-i b (c) t (d) Sij s i-i,j and Sij a i J + n d s i,j > Sij + i > S i and Sij i+l,j + These criteria allow the ridge points to be located adjacent to one another along the top of coherent and significant trends in the similarity matrix. similarity matrix for example m /\ 2 Figure 3.5 shows the with one contour at 0.5 and ridge points as x's. 3.5 S u b p a t h D e t e r m i n a t i o n ( D y n a m i c M a t c h i n g ) Establishing a set of potential pieces of the correct matching function requires joining the ridge points into 'subfunctions' or 'subpaths'. This is done with a simple dynamic programming scheme. A ridge point is defined to be a member of a subpath if (l) it is within a prescribed distance D from the previous ridge point and (2) the slope of the segment that joins them is non-negative. The subpaths are thus constrained to be non-decreasing. Figure 3.6 shows the local path constraint used for subpath finding. The algorithm to determine the set of subpaths proceeds as follows. T h e ridge points are sorted by their distance from the origin, and subpath finding begins with the ridge point nearest the origin. If the distance d from the last ridge point in the subpath to the next nearest ridge point exceeds D, that is if d > D, then the subpath is terminated. W h e n the first subpath is terminated, all ridge points belonging 49 REFERENCE Figure 3.5 Location of ridge points as x's with one contour at 0.5. 50 F i g u r e 3.6 T h e local path constraint used in dynamic matching to connect ridge points into subpaths. Ridge points 1 and 2 are legal choices, ridge point 3 is not. to it are removed from the set of 'possible subpath points', and the next ridge point nearest the origin begins a new subpath. The procedure is repeated and all subpaths have been found when the set of possible subpath points is empty. Also, if a subpath is less than a prescribed length L it is deemed insignificant and is removed from the set of subpaths. T h i s puts a m i n i m u m significant length on the subpaths as prescribed by the definition of local confidence discussed in section 3.3. Figure 3.7 shows the ridge points of Figure 3.5 connected into subpaths with D — 3.5 and L = 4.2. These parameter values were established from trial and error and may be converted to well log samples for generality by multiplication by A x . More will be said on the determination of these parameters in Chapter 4. 51 REFERENCE F i g u r e 3.7 T h e connection of ridge points into subpaths with dynamic matching. Ridge points were connected using D = 3.5 and L = 4.2. 52 3.6 The Function Network (Multiple Dynamic Matching) The set of subpaths found by dynamic matching must now be combined into a set of complete matching functions or 'function network'. Since confidence in a match is low between subpaths, linear approximations are made to join the subpaths into complete functions. To determine which subpaths to join, only those which when joined preserve a non-decreasing function are considered. Joining a subpath with its nearest non-decreasing neighbor may be too restrictive, so the nearest two subpaths are included as possible. The algorithm to build complete functions from subpaths is nontrivial when two choices are allowed at the end of every subpath. The approach taken to build the function network is to view each subpath as a state and allow two possible transitions from each state. A transition network is established in this manner, where in the case of subpaths, the two next nearest subpaths are the transitions allowed from each subpath. The subpath nearest the origin starts the function building process. A new function is added for each branch from a subpath, and has attributed to it, all the subpaths of the function built up to that branch. Paths are terminated with a subpath that has no possible next subpaths. To ensure that all paths have been found, all subpaths are tested for their membership in a complete path, and any subpath that is not part of a complete path becomes the starting subpath for a new set of paths. A n estimate of the total number of paths is discussed in section 3.9. Figure 3.8 shows the function network for the set of subpaths shown in Figure 3.7. 53 REFERENCE F i g u r e 3.8 T h e function network found from the subpaths in Figure 3.7. Each subpath may branch to the next two nearest subpaths. There are a total of 180 functions. 54 3.7 R a n k i n g the Functions T h e set of possible matching functions found from multiple dynamic matching must now be ranked to determine the optimum matching function. 3.7.1 M e a n Local Confidence It was argued in section 3.3 that significant trends in correlation appear as well defined ridges in the similarity matrix. T h e concept of local confidence was developed to measure the quality of the local matches, and now that a set of complete functions has been found, the mean local confidence ( M L C ) norm may be computed to measure the quality of each entire matching function. Equation (3.4) 1 MLC = is restated here. P - — E S x,m(x)5" m ( a ; ) x. p= 0 It requires computing S x<Tn () x and — S' x m ( )± x along each function m(x) which must be stepped along at a fixed increment to avoid biasing. Since this will generally lead to coordinates not on a lattice point of the matrix, interpolation is necessary. Simple distance-weighted bilinear interpolation is employed to compute S x<m (y x and — S" ^ . m ± T h e M L C norm thus requires interpolation of values at the coordinates (x, ra(x)) , (x — 1 1/2 A x , m ( x — A x ) + A y ) , and (x + A x , m(x + A x ) — A y ) , where ( A x 1 + Ay ) = 1. A three point finite difference is then used to compute the second derivative part of local confidence as: — Sx m(x) t x ~ Sx-Ax,m(x-Az) + Ay ~ 2 - S ( -) + S XTn x x + Ax,m(x-|- A x ) - A y (^-5) 55 A l l matching functions are now ranked by their M L C ' s to determine the best functions. Figure 3.9 shows the top ten functions based on M L C and the top ranking function as a thickened line. Figure 3.10 shows the results of applying the top three matching functions to the test log. Due to the relatively high amplitude distortion a(x) between well logs, additional constraints are helpful to determine the interpreter's optimal m(x). 3.7.2 Simplicity Other properties beside the mean local confidence are important in determining the optimum matching function. solution. The goal in most problems is to find the In the case of matching well logs the simplest structural produce optimal matching is a reasonable simplest change that will requirement. A measure of simplicity which is easy to compute is the ratio of linear path length to total path length. T h e linear path length may be easily approximated by the distance from the first to last points in the matching function. The simplicity (SIMP) norm for the kth matching function with ridge points j = 1, M is: S I M P = K*M - si) + ( " ' t ( x L ) - ^ K ) ) ] ( 3 6 ) The top 10 functions as ranked by M L C are then ranked by S I M P . T h e ten matching functions after this ranking are shown Figure 3.11, and the top three are shown applied to the test log in Figure 3.12. 56 REFERENCE F i g u r e 3.9 T h e top ten functions based on M L C with the top ranking function thickened. Figure 3.10 The results of applying the top 3 matching r=reference, t=test, l,2,3=matching function ranking. functions to the test 58 REFERENCE F i g u r e 3.11 T h e top ten functions after ranking by M L C and S I M P . The top ranking function is shown thickened. Compare with Figure 3.9. 59 F i g u r e 3.12 T h e top three functions after M L C and S I M P ranking applied to the test log. Compare with figure 3.10. 60 3.7.3. Sufficient D i v e r s i t y Supplying alternative interpretations requires that they be sufficiently different from one another if they are to provide useful information. T w o of the matching functions in the top three are usually so similar that another more diverse interpretation is desired. To ensure the sufficient diversity of the optimal matching functions, some measure of difference between them should be exceeded. A measure of difference between the best ten matching functions as ranked by M L C and S I M P would be: x Dij = \ i{ ) ~ j{ )\ m x m x f o r * = !> 1 0 '. 3=* + 1,10- (3-7) Cluster analysis may be done using this dissimilarity matrix to separate similar matching functions from sufficiently diverse ones. The top ranking functions from each of the clusters may then be taken to form the set of best three sufficiently diverse functions. Instead of implementing cluster analysis, the half-matrix of equation 3.4 was computed for the top ten functions in example m /\, and new functions were brought into 2 the top three if any two showed Dij < 2.0. The result of implementing this method of sufficient diversity is shown in Figure 3.13. The three functions now provide the interpreter with more information. 3.8 A p p l y i n g the M a t c h i n g F u n c t i o n Once the best functions have been found, a visual check as to the quality of the matches is useful. Figures 3.10, 3.12 and 3.13 were produced by applying the top three matching functions (after different rankings) to the test log. The application of a matching function requires the stretching or compressing (warping) of the log axis as F i g u r e 3.13 T h e best three functions applied to the test log after they are constrain to be sufficiently diverse. 62 determined by a segment of the piecewise linear function. If the local slope exceeds 1.0, the test log is locally compressed, and if the local slope is less than 1.0, the test log is locally stretched. Then the log is resampled to maintain a constant sampling interval. T h e resampling of the log after warping the signal axis requires the inperpolation of log values. Interpolation has been done for examples in this project using a simple distance - weighted linear scheme. filtering It has been found quite satisfactory given the rather heavy of the sonic data, but for higher resolution matching some other (te. cubic spline) interpolation may be needed. Another method of display, other than applying the matching function, is to connect matched locations in each log with a dashed line (see Figure 1.2). Graphically this is equivalent to locating a point on the axis of one of the logs and projecting it on to the matching function. This point is then projected on to the other signal axis to determine the matched location. This mode of display is particularily useful when matching structurally equivalent portions of well logs and will be used in section 3.12. 3.9 Computing Time Factors that affect the computing time necessary for a solution are the similarity matrix bandwidth b, the matrix sampling increment ( A x ) , and the number of subpaths (n). The matrix bandwidth 6 is set by the user based on a priori information about the complexity of the geology over the separation bound is easily estimated. of the wells. Some practical upper The size of the matrix is related directly to the bandwidth 6 and inversely to the matrix sampling increment A x . Since higher resolution means a shorter low-pass median filter length and finer matrix sampling, and since matrix size is proportional to ( 6 / A x ) , the cost of the algorithm rises rapidly with resolution. 2 63 The effects of bandwidth and resolution on computing time are easy to quantify. Another cost factor is the total number of matching functions found. This is an i m - portant factor of cost since the M L C norm is computed for each function. T h e total number of functions is related to the number of subpaths. T h e number of subpaths in turn is a function of the local path constraints and the underlying similarity of the signals being matched. Predicting the number of subpaths is therefore a difficult task. Given that the number of subpaths is known however, the total number of functions may be roughly estimated. The function building algorithm progresses forward through the similarity matrix, combining subpaths with two others at each 'level' of the network. A s levels are passed, all subpaths behind the forward most points in the functions are eliminated from possible combination. T h e power to which 2 should be raised is thus reduced by the number of subpaths passed at each level. T h e number of subpaths passed at any level is related to the spatial period of the ridges and the bandwidth of the matrix. Once the number of ridges r across the matrix is estimated, the number of paths for each starting subpath is then estimated by 2 ^ " / r _ ^ e _ 1 ^ where n is the number of subpaths and e — 1 is subtracted from the power for the subpaths in the ending level. This number is then multiplied by the number of starting subpaths. In example m / i . w i t h n = 30, r = 4, and e = 3, the estimate of the total number of paths 2 is 128 and the actual number is 135. Memory requirements should also be mentioned. T h e similarity matrix Sij is not fully used since only the central band needs to be stored for pathfinding. A more efficient use of memory would eliminate the unused 'triangles' above and below the central band, making Sij rectangular, but this would require a rotation of the matching function before application. T h i s saving in memory has not been implemented in the interest of preserving clarity in similarity matrix figures. 64 3.10 Sources of E r r o r The error in any matching function may only be thought of relative to the resolution or scale of structures left in the logs. A function which is correct at one scale will probably not be correct at another scale. T h e matching of significant events - the block edges or steps, is the most important aspect of matching and will be considered first. The error in m(x) for any x is at a minimum when local confidence is at a maximum. This occurs over subpaths where ridge points have been located. A s explained in section 3.3, ridge points (and subpaths) are found where events of significant magnitude have been matched, so the matching of major events is done with minimal error over subpaths. The errors in matching over subpaths are due only to the matrix sampling increment. A n y significant block edge will be matched to within A x , but ± A x / 2 is a better estimate. The greatest errors occur where confidence is low. This happens over the linear approximations between subpaths. Over these portions of the function, error in matching is not easily quantified in terms of samples. However, a plot of the local confidence vs reference depth can indicate where errors are greatest. m / i ( x ) and L C vs 2 shows match x. Another source of error should be mentioned also. filtering Figure 3.14 does not accurately The degree to which median represent the underlying blockiness of the original signal can also induce errors in the match. Matching errors of this type and those due to low local confidence may be reduced by increasing resolution and matching finer scale structures. Using a matching function found from broad scale matching as a constraint on finer scale matching can reduce errors and make matching of small scale structures possible. 65 0.4 - 0.1 ] 0 F i g u r e 3.14 , 20 , 40 1 , 60 80 L o c a l confidence vs. REFERENCE AXIS , , 100 120 1 1 140 x for the top ranking matching function in figure 3.12. 3.11 Multiple Dynamic Matching and Pattern Recognition Exploring the analogy between multiple dynamic matching or ' M D M ' and pattern recognition affords insight into how M D M works. A t a given scale obvious matches of major events 'jump out' at the interpreter. He/she scans along the signals thinking something like 'this matches with this and that matches with that...' etc. This pattern matching process is represented in the similarity matrix by coherent ridge trends where each ridge represents a distinct local interpretation. When the interpreter's thought becomes 'or that could match with that', this is equivalent to skipping to a new ridge in the similarity matrix, and another interpretation is included as possible. M D M shows this doubt in interpretation as a branch in the 'interpretation tree' of most likely matching functions (see Figure 3.8). A l l the most probable interpretations are found 66 automatically by this process of connecting local interpretations together. Finally the MLC and S I M P norms mimic the confidence criteria of the expert and the best match is chosen. The analogy may be continued through scale as well. Just as the interpreter matches 'marker' events first and then looks for subtleties, so may M D M be able to 'fine tune' a broad scale match in subsequent processing. High resolution matching down to a few samples should be possible using a 'scale down' procedure where the m(x) at one scale is used as a constraint on the m(x) at the next scale down, effectively restricting the bandwidth of Sij and so also the possible change in ra(x). 3.12 Automatic Stratigraphic Correlation The major objective of matching well logs is to ascertain an estimate of the geologic change between wells. T h i s means estimating the local dip or structure through the correlation of equivalent time horizons or events, and estimating the change in the measured response. and The local dip or structure is given, as is discussed in section 1.5, by m(x), while the change in the measured response is given by a(x). 1.2 Since the matching function quantifies the correlation between signals, it may be used to connect structurally equivalent depths in two wells. To illustrate this connection, matching functions were found between the five sonic well logs in Figure 2.12. T h e method of automatic event detection discussed in section 2.5 was implemented to find the major events in sonic log 1. Figure 3.15 shows these major events followed between wells using the matching functions for logs 1 and 2, 2 and 3, 3 and 4, 4 and 5. Each of the matching functions used was top ranked based on the M L C and S I M P norms, and since event 67 detection was also done automatically, the resulting structural interpretation is totally objective. F i g u r e 3.15 The automatic correlation of five well logs from four matching functions. M a j o r events were found in log 1 and followed across the section. Equal spacing of the logs is arbitrary. 69 Chapter 4 CONCLUSIONS 4.1 Improvements T h e method of multiple dynamic matching has been shown to be successful in the problem of matching trend - removed and median filtered sonic well log data. Never- theless some improvements are needed to increase matching resolution and to further generalize the method. A s discussed in section 3.9, the error in the match is greatest where linear matches are assumed between subpaths. Reducing the errors in the matching function requires higher resolution matching, but as discussed in section 3.9, this causes a rapid increase in the cost of the algorithm. One possibility that should allow high resolution matching at reasonable efficiency is that of matching through scale. The idea is to obtain a broad scale matching function, apply the function to less severely filtered versions of the logs and proceed to determine the higher resolution matching function. Since the broad scale estimate of the function has been applied to the data, it has been mapped to the diagonal of the next higher resolution similarity matrix which is now square, and the matrix bandwidth may be reduced dramatically since large scale structural changes have been removed. Thus matching at higher resolution may be done more efficiently due to the 'smaller' matrix. The matching function at the next resolution may be viewed in terms of its deviations from the diagonal, and once the previous 70 warping of the test signal axis has been accounted for, the deviations may be added to the broadest scale function. Resolution in the matching function may be iteratively built up by this method, and with the bandwidth of the similarity matrix reduced at each iteration, the cost is kept at a m i n i m u m . There is a problem with this approach however. The broad scale match has been obtained from filtered versions of the logs that have flat regions (constant neighborhoods). Sections of severe stretching and compressing in the matching function when applied to a fiat region of the test log, do not affect the quality of the match over these regions. Unfortunately, severe slopes in a matching function when applied to a less filtered test log will produce pronounced and unreasonable distortions in the signal, making it unsuitable for further matching. T h u s the application of a broad scale matching function may not be appropriate before higher resolution matching. T w o possible alternatives may help this situation and allow for iteratively matching through scale. The first involves smoothing the broader scale function before appli- cation, reducing the severe slopes and so also any unreasonable distortions in the less filtered logs. The second alternative does not apply the broader scale mapping at all, but instead, a narrow band of a higher resolution similarity matrix would be computed around the already determined broad scale function. The next matching function is thus restricted to be close to the broad scale function. This approach complicates the computation of Sij but may prove superior to smoothing the matching function before its application. It is believed that higher resolution matching is possible with either of the approaches discussed above, but they must be tried and assessed based on the usual criteria of accuracy and efficiency. 71 4.2 Parameter Estimation The method of multiple dynamic matching requires the fixing of the parameters N (half the low-pass median filter length), K (the similarity window length), and D (the proximity threshold for subpath determination). A rough estimate of K has been established (given that A x « 2/ZN) as RS 4N or 6 A x , as this provides a sufficiently smooth similarity matrix with the proper character. The median filter length (2AT + 1) determines the resolution of the match and is selected by the user (as is the matrix bandwidth 6), but the proper choice for D is difficult to establish given the other parameters. The proximity threshold D is the maximum allowable separation of ridge points in any subpath, and although a value of 3.5 (or 35 log samples) has worked well for the matching examples in this project and tests have indicated that broader scale matching using N — 29, A x = 20 K '= 120 and D — 35 have produced similarily good results, the value of D in other matching problems will likely have to be altered to suit the particular application. The generality of the method thus breaks down if the local path constraint (D) cannot be determined automatically in every application. one subpath becomes two is a crucial question in multiple dynamic matching. Where The method requires enough subpaths to fully explore the space of possible matching functions, but not too many to make the solution prohibitively costly. Intuitively it seems that there should be a range of values for D for which the number of subpaths should remain constant. This kind of testing could be done to establish D for matching. The m i n i m u m subpath length L has been fixed at 42 log samples and this has also proven to be a good value at different scales. B u t the fixing of D, L and K will necessarily depend on the particular application and perhaps most strongly on the bandwidth and nature (eg. blocky or continuous) of the input signals. The automatic fixing of these 72 parameters and ultimate generalization of M D M requires experiments with other data sets and other matching problems. 4.3 Speculations on the Future a n d Related Ideas 4.3.1 E x t e n d i n g Well L o g M a t c h i n g to 3 D T h e structural interpetation of Figure 3.15 equal distance appart. shows the logs arbitrarily spaced an In reality they are distributed as shown in Figure 4.1. It is evident that only 4 matching functions have been found and that there are still 6 remaining to be determined since there are a total of 10 possible matches. Figure 4.1b shows the 'net' of matching combinations possible with this areal distribution of five wells. The intersection points indicate locations where more control is possible on the structural interpretation. For each major event a surface may be approximated allowing errors at the intersection points. In this way an objective three dimensional interpretation is possible using only well log information. More complex geologies would require more sophisticated fitting schemes, and due to nonuniqueness, more than one interpretation should be generated, but each interpretation may be ranked based on optimality criteria, thus giving the interpreter an objective basis with which to compare his own interpretation. Once a three dimensional structural interpretation has been settled upon, by whatever method, the spatial variation of the measured response, (ie. may be estimated. acoustic velocity), The determination of matching functions for a set of well logs can therefore automatically estimate the structural and response change in an area. 73 1 1 F i g u r e 4.1 (a) T h e areal distribution of wells and the matches corresponding to those in Figure 3.14. (b) T h e matching 'net' possible with the five wells. Intersection points indicate where more control is possible for the mapping of structural events over an area. 4.3.2 Seismic Signal M a t c h i n g In seismology, the importance of signal correlation was recognized early on, especially for the determination of velocity functions. Indeed, the interpretation of seismic sections depends on the spatial correlation of events at different times. The spatial - temporal correlation function implicitly determined by the interpreter approximates the geology, so any artificial intelligence system that attempts to automatically interpret a seismic section must first determine the correlation functions between traces. T h e matching of seismic traces is a considerably easier problem than matching well logs, owing to the closer proximity of adjacent measurements and the lower frequency content of the signals. These factors mean that the matching function is constrained to be quite simple. In the context of multiple dynamic matching this means a narrower similarity matrix and efficient high resolution matching. The multiple interpretation 74 advantage of M D M could be effectively used in seismic interpretation, since locations where a second or third interpretation is probable may be nagged and further investigation focussed on that area, perhaps allowing the computer to continue the interpretation with the second best choice. Nonlinear signal matching thus provides automatic picking of events across a section. A s a bonus, storing adjacent matching functions gives the spatial - temporal correlation function, making improved trace interpolation possible. Also, multichannel deconvolution, dip moveout filtering, and automatic velocity analysis can be approached from a new and different point of view. 4.4 Summary T h e nonlinear signal matching problem has been introduced by contrasting it with the simple idea of cross correlation, and the discrete lag in the cross correlation expression was replaced with a continuously varying lag or matching function, m(x). The matching problem was further explained by a model of change between waveforms that decomposes the problem into a structural change m(x), a(x). and an amplitude distortion It was argued that in the presence of large amplitude distortions, the correct matching function may appear suboptimal under similarity norms and as large amplitude distortions are likely to occur in any difficult matching problem, a method is needed that can provide suboptimal as well as optimal solutions. T w o distinct methods for the solution of the nonlinear signal matching problem were reviewed, namely the dynamic programming approach and the inverse theory approach. It was found that both methods have advantages and pitfalls, and that neither method is satisfactory when the goal is to objectively determine optimal and suboptimal solutions. 75 Multiple dynamic matching was introduced as a method that overcomes the pitfalls of both of the other methods. It is founded on a dynamic programming philosophy in that pathfinding through a similarity matrix is done, but with the innovation of a significance threshold to extract only the matches where confidence is high. Each of these 'subpaths' is then joined in combination with others to generate a network of possible complete interpretations. T h e advantage of the inverse approach is that optimality is based on global rather than local norms, and this feature has been incorporated into M D M in that the objective function to be maximized measures the quality of the entire matching function. Due to the large amplitude distortions between well logs the correct or interpreter's matching function may be ranked as suboptimal, so a simplicity constraint can be invoked to improve the ranking of the 'interpeter's choice'. The results of M D M are good for the example explored here, and given the difficulty of the well log matching problem, the method should perform well in other simpler matching problems. It was also conjectured that higher resolution matching with reduced error would be possible by matching through scale, where a broader scale match is used as a constraint on finer scale matching. T h e 'scale space' approach is receiving considerable interest in the artificial intelligence fields of computer vision and pattern recognition, and since signal matching is essentially a pattern recognition process, matching through scale is a promising prospect to improve the results of multiple dynamic matching. In fact the idea of scale space was discussed in chapter two, and it was argued that the median is a superior parameter of scale than that employed by Witkin (1985), especially when looking at blocky waveforms. Median scale space and the new signal decomposition presented for the first time in section 2.4, median decomposition may prove to be suitable for the description of other intricate waveforms as well. Although the subject of median filtering was secondary in the present work, it represents, as does multiple dynamic matching, the need for new nonlinear approaches to solve difficult problems in signal processing. 77 Biblography Anderson K . R . and G a b y J . E . , 1983, D y n a m i c waveform matching, Information Sci- ences,31, p221-242. Bednar T . B . , 1983, Applications of median filtering to deconvolution, pulse estimation, and statistical editing of seismic data, Geophysics 48, No. 12, pl598-1610. Claerbout J . F . and M u i r F . , 1973, Robust modelling with erratic data, Geophysics, 38, p826. Clarke, R . M . , 1985, A fortran program for constrained sequence slotting based on minimum combined path length, Computers and Geosciences, 11, N o . 3, p605-617. Bellman R., and Dreyfus S., 1962, Applied dynamic programming, Princeton U . P . , Princeton N . J . Douze E . J . and Laster S.J., 1979, Statistics of semblance, Geophysics, 44, No. 12, pl999-2003. G o r d o n A . D . and Reyment R . A . , 1979, Slotting of borehole sequences, Journal of M a t h ematical Geology, 11, N o . 3, p309-327. Evans J . R . , 1982, Running median filters and a general despiker, Bulletin of the Seis- mological Society of America, 72, N o . 1, p331-338. F i t c h J.P., Coyle E . J . and Gallgher N . C , 1984, Median filtering by threshold decom- position, I E E E transactions on Acoustics, Speech and Signal Processing, A S S P - 3 2 , No. 6, pll83,1188. 78 L a n n i n g E . N . and Johnson D . M . , 1983, Automated identification of rock boundaries: A n application of the Walsh transform to geophysical well log analysis, Geophysics, 48, N o . 2, pl97-205. L i u H . H . and F u K . S . , 1982, A syntactic pattern recognition approach to seismic discrimination, Geoexploration, 20, pl83-196. Martinson D . G . , Menke W . and Stoffa P., 1982, A n inverse approach to signal correlation, Journal of Geophysical Research, 87, N o . B6, p4807-4818. Myers C S . and Rabiner L . R . , 1981, A comparitive study of several dynamic timewarping algorithms for connected word recognition, Bell system technical journal, 60, N o . 7, pl389-1409. Neidell N . S . and Taner M . T . , 1971, Semblance and other coherency measures for multichannel data, Geophysics 36, p482-497. Rabiner L . R . , Rosin berg A . E . and Levinson S . E . , 1978, Considerations in dynamic time warping of discrete word recognition, I E E E transactions on Acoustics, Speech and Signal Processing, A S S P - 2 6 , N o . 6, p575-582. Robinson J . E . , 1978, Pitfalls in automatic lithostratigrap hie correlation, Computers and Geosciences, 4, p273-275. Schwarzaeher W . , 1980, Models for the study of stratigraphic correlation, Mathematical Geology, 12, No. 3, p213-234. Smith T . F . and Waterman M . S . , 1980, New stratigraphic correlation techniques, Journal of Geology, 88, p451-457. 79 Stewart R . R . , 1985, Median filtering: review and a new F / K analogue design, Journal of the Canadian Society of Exploration Geophysicists, 21, N o . 1, p54-63. W i t k i n A . P . , 1984, Scale space filtering: a new approach to multi-sscale description, in Image Understanding, Norwood N . J . , Ablex Publishing Corporation. 80 Appendix T h e pages that follow contain the similarity matrices and ten best matching functions after ranking by M L C . The ranking of the chosen function is given in the captions of the ' b ' figures for each pair of logs. If the top 10 functions are ranked also by the S I M P norm, the chosen interpretation is then ranked as number 1. The chosen functions have been used to generate the structural cross section of Figure 3.15. 81 20 40 60 80 100 120 140 160 REFERENCE F i g u r e A l . a Similarity matrix for reference = log 3 and test = log 2, and the best 10 functions after ranking by M L C . T h e top ranking function is shown as a thickened line. 82 F i g u r e A l . b T h e top three functions after M L C diversity applied to the test log. The number 2 ranking function was chosen. 83 REFERENCE Figure A 2 . a Similarity matrix for reference = log 3 and test = log 4, and the best 10 percent of functions after ranking by M L C . T h e top ranking function is shown as a thickened line. 84 F i g u r e A2.b T h e top three functions after M L C applied to the test log. The number 1 ranking function was chosen. 85 REFERENCE Figure A 3 . a Similarity matrix for reference = log 5 and test = log 4, and the best 10 percent of functions after ranking by M L C . T h e top ranking function is shown as a thickened line. T h e matrix bandwidth is much narrower due to the close proximity of the wells. 86 ) 400 800 1200 1600 111111111111111111111111111111111111111111111111111H1111111111111111111M11111111H111111111 Pi T A I I I I II I I I I j M I I II I I I I I I I I I I I I III I I H I I I I I I F i g u r e A3.b T h e top three functions after M L C applied to the test log. T h e number 2 ranking function was chosen.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Multiple dynamic matching and median filters applied...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Multiple dynamic matching and median filters applied to sonic well logs Leaney, W. Scott Powell 1986
pdf
Page Metadata
Item Metadata
Title | Multiple dynamic matching and median filters applied to sonic well logs |
Creator |
Leaney, W. Scott Powell |
Publisher | University of British Columbia |
Date Issued | 1986 |
Description | Nonlinear signal matching is a generalization of cross correlation in that a discrete lag between signals is replaced with a variable lag function or 'matching function', m(x). Two methods are reviewed which attempt to solve for m(x), namely the dynamic programming approach and the inverse theory approach. Both methods suffer from pitfalls and require the input of prior constraints to ensure convergence to the correct solution. The goal of this work has been to develop a method that can handle simple or complex matching problems and can succeed without any prior knowledge constraints. The multiple dynamic matching method is the result. It uses a significance threshold to extract a set of ridge points from a similarity matrix and applies dynamic programming to obtain a set of matched sections. These significant matched sections or 'subpaths' are then combined into a set of complete matching functions and a 'mean local confidence' norm is evaluated to determine the best overall match. It is argued through a model of change between signals, that given that the correct matching function is known, the presence of large amplitude changes between signals can cause the correct matching function to appear suboptimal under similarity norms. Multiple dynamic matching, because it generates suboptimal solutions as well, does not overlook the correct matching function. Typically the top three interpretations as ranked by mean local confidence will contain the expert's choice for the correct matching function. The use of median filters to preprocess the data and enhance well log features for matching has been investigated. A new 'median decomposition' is discussed as well, and in the context of a scale - space point of view to filtering, it is argued that median scale space is the proper choice for blocky waveforms. Finally, the connection between multiple dynamic matching and pattern recognition is explored, and matching iteratively through scale is proposed as a means of further generalizing the multiple dynamic matching method, making efficient high resolution matching possible. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-06-20 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0052980 |
URI | http://hdl.handle.net/2429/25900 |
Degree |
Master of Science - MSc |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1986_A6_7 L42.pdf [ 3.45MB ]
- Metadata
- JSON: 831-1.0052980.json
- JSON-LD: 831-1.0052980-ld.json
- RDF/XML (Pretty): 831-1.0052980-rdf.xml
- RDF/JSON: 831-1.0052980-rdf.json
- Turtle: 831-1.0052980-turtle.txt
- N-Triples: 831-1.0052980-rdf-ntriples.txt
- Original Record: 831-1.0052980-source.json
- Full Text
- 831-1.0052980-fulltext.txt
- Citation
- 831-1.0052980.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052980/manifest