Introduction
A fascinating signal processing problem that is pertinent to many other applications often comes up in oil exploration. When looking for oil, the prospector often needs to know the slope of the underground formations. The oil, being lighter than underground water, filters upwards, to the top of underground domes. If you poke a hole near the top of such a dome, capped by salt or shale — you hit oil! If you miss the top of a dome, and have a dry well, you will want to know which direction will take you up the formation's slope. You move in this direction and drill your next well.
To find the slope of the formation, prospectors use a tool called a dipmeter. This tool simultaneously measures the resistivity of the surrounding formation at a few (usually four or six) places around the circumference of the bore hole. If the formation is dipping, measurements from two different probe sites look similar, but are shifted in relation to each other. From the displacements between the curves, the prospector can calculate the dip.
Astute readers will recognize in this example a more general signal processing problem: two time-varying signals from the same source each undergo separate, distorting processes. The processes are independent and each one distorts its respective signal by a time-varying amount, in both amplitude and delay. By comparing the two resulting signals, can you find corresponding points on each signal, and determine the relative delay between signals at each point?
This kind of problem can arise when trying to correlate signals that have passed through different media, or, say, when trying to register two graphic images. You may have encountered similar problems in your own applications. I would welcome hearing about your own approach to the problem.
The Dynamic Programming Approach
The dynamic programming approach involves three steps:
- Select "events" on
each signal.
- Find all possible matches
between events on each signal.
- Determine the best set of
non-contradictory matches.
The relative shift or delay between signals at each match point issues as a side-product of the calculations.
Detecting Events
But what do I mean by an "event"? An event is a feature on each curve that is likely to have been produced by the same physical event — the dipmeter passing through a formation, a cosmic-ray burst in signals from a distant galaxy, etc. There are at least three ways to pick corresponding events.
- a) Select fixed-length intervals on the first curve and match them as a whole to the second curve. This involves sliding a "window" along the second curve and computing a correlation coefficient for each possible displacement.
- b) Find inflection points on one curve and match them to inflection points on the other curve.
- c) Use segmentation and data trees [2] .
Each approach has its pros and cons:
- a) Selecting an interval does not require any special technique. Since the intervals are selected at regular steps, the interval selection does not depend upon the character of the signal. Thus, this technique produces correlations that are averaged and not visually discernable. This form of curve matching does not provide a good example for the sake of illustration.
- b) To find inflection points a
program can use the signal derivative. The program
selects points where the derivative reaches local
maxima and minima. It may be a good idea to average the
derivative somewhat to get rid of some noise. The
activity function described in [1]
provides a good means of noise filtering. From an
illustration standpoint the inflection points are easy
to see, so I have selected this method for my example
below.
- Note, however, that in general,
correlation on inflection points requires good quality
data, with low noise, and of sufficient frequency to
produce a statistically useful sample of points. Also,
you don't have to use inflection points; they just
happen to be significant to geologists. In your
application other features of the curves, say, peaks
and troughs, may be more significant.
- c) The segmentation trees
technique is a good method for finding events. However,
its discussion would take us away from the main goal of
this article. You can find more on this in [2]
.
Finding all possible matches
It is easy to find pairs of all possible matches. For every inflection point on the first signal we check all inflection points on the other signal within a certain distance. To be a possible match, the point on the second curve must satisfy the following conditions:
- It must be within a certain
distance from the first point, since in every practical
application there is a limit on how far along the curve
we should search. This limit is appropriately called
the search limit. If the two signals vary with time,
the search limit extends forward and back on the time
axis. In geology the signals vary with probe depth, and
the search limit is determined by the maximum dip
angle expected in this area.
- If the point on the first
signal is on the up slope, the point on the second
signal also must be on the up slope, and vice versa.
- The signal slopes should be similar. We satisfy this condition by requiring that the derivatives do not differ beyond a predetermined amount. (In some applications, it's appropriate to allow any slope to match any other slope, as long as the two are going in the same direction. In this case the derivative difference limit does not apply.)
Coding the checks above is relatively easy; hence, I do not give the code for this part. Figure 1 depicts the situation where all possible matches have been found.
Paring Down Matching Events
Now comes the heart of the algorithm. From all the possible matches in Figure 1, how do we select the best set?
The answer lies in the following simple model. Imagine that the events on the second curve are connected with springs (see Figure 2) . Before the matching starts, there is no tension in the springs; their unstressed lengths exactly fit the distances between points.
Now imagine that we start fitting the second curve to the first. If the curves are the same, we get a perfect fit. The points on the second curve match exactly to the points on the first one, and the springs do not stretch or compress.
However, the curves are not the same, so we must stretch or compress the second curve in some places. The distances between the points on the second curve start to change, the springs start moving, and a certain amount of tension appears in them. It seems natural to try to find the fit in which the total tension of all the springs would be minimized. This is exactly the solution provided by the algorithm below.
There is a minor complication, however. Not all points on the second curve may find a match. The algorithm must be able to handle a point with no matches. On the other hand, too few matches overall will produce an unsatifactory data set. Somehow the algorithm must be discouraged from rejecting too many less-than-ideal matches.
To get enough matches, I add a penalty to the total spring tension for leaving out any match. The amount of this penalty depends upon how good the rejected match is. The penalty makes the algorithm tend to keep as many good matches as possible.
(By the way, in dipmeter processing, if the charge is by the number of computed points, each match brings in about $1, so there is a good incentive to come up with many matches.)
Mathematical Description
In this section I present some formulas to explain elastic curve matching in more detail. Since understanding these formulas is not strictly necessary to using this technique, you can skip this section, as well as the next one, if you don't like math.
Let (Xi, Yi, Ci) denote a set of all possible matches, where i varies from 1 to N. For each i, Xi is the coordinate of the point on curve 1, Yi is the coordinate of the possibly matching point on curve 2, and Ci is the penalty for leaving this match out of the solution. (If you do not have a good way to determine the set of penalties Ci, you can set all of them to 1.)
Since a point on one curve may have a number of possible matches on the other curves, some of the Xs may be the same. However, all pairs (Xi, Yi ) are different for different i.
The tension of the spring can be roughly measured as the difference between its initial length and its final length. The tension between two adjacent points is thus equal to
|(Xj - Xj+1) - (Yj - Yj+1)|
The algorithm must select a set of indices j to minimize the total tension, while being penalized for each skipped correlation. In other words, the algorithm attempts to minimize the expression
where j denotes all indices from 1,N that are selected, k denotes all indices not selected, and E is the elasticity coefficient for the springs.
Dynamic Programming for Minimization
The problem of minimizing (1) can be solved using the mathematical technique known as dynamic programming. Dynamic programming was invented by the American mathematician Richard Bellman in 1957, and is used in operations research.
The dynamic programming algorithm consists of two steps. In the first step, the algorithm computes the best possible value of the optimization function at each intermediate point of the solution. In the second step, the algorithm goes back and decodes the path that led to this optimal solution.
As an illustration, consider a map of cities, connected by roads. Suppose that a driver must go from city A to city B in the shortest amount of time. The travel time between every two connected cities is known.
The dynamic programming algorithm will first find travel times to all cities immediately connected to A. Then it will check all two-leg journeys departing from A and select the quickest. For example, if there are two ways to get to city C, through D or through E, it will compare travel time AD + DC to AE + EC and select the best of the two. Then the algorithm checks all three-leg journeys and so on, until it has checked all possible journeys departing from A. Of course, when a journey ends at B, the algorithm does not check the remaining legs of the journey.
By now the algorithm knows the quickest time to get to B. To find the path, it will first see which was the previous city leading to B that allowed for this best time. Then the algorithm will find the city on the best path before this, and so on, until it gets back to A.
This is pretty standard, and you can find more detail in any book on the mathematical optimization or operations research [3] . The only part of our particular problem where we need some inventiveness is in finding the first and the last points in the solution.
In the example above, point A and points B were given. By contrast, for curve matching, we need the freedom to skip any number of points from the beginning of the interval, and any number of them from the end — as long as we find a good number of matches in-between.
To allow the algorithm to start matching in any place, I add two "virtual" points on both curves, one at the beginning of the curve and one at the end. (I have not seen this anywhere else in the literature). In the intermediate stages, both left and right virtual points are part of the solution, and the dynamic programming algorithm skips any number of points in between as it works. When I decode the path, however, I do not show these virtual points as part of the solution, so matching starts in any place.
The Algorithm in C++
The optimization function is called Select (Listing 1) and is straightforward. (Please forgive me for counting indices between 1 and N, and not between 0 and N-1; the function was originally coded in FORTRAN. Besides, since the two virtual points have to be prepared by the caller, they change the indexing anyway.)
The function expects all points to lie within the 1 through N positions of the D and X arrays. Points 1 and N should be virtual. Select's caller need not take any action to create the virtual points; it just provides the space for them. On the output, N is changed to the number of found matches, and arrays X and D contain the answers: X[i] matches to D[i] for i = 1,N.
The results of using Select are shown in Figure 3.
Generalizing the Optimization Function
With just a little effort, the function above can be generalized to allow for optimal selection of any objects, not just pairs of curve matches.
For example, you may need to select from a basket the best set of apples for an apple pie. Ideally you would pick all the apples, but they contradict each other in taste and color. If you can measure the amount of contradiction between every two apples, and the emotional cost of leaving out each apple, then your problem can be solved by the second version of the Select function.
All you have to do is derive your Apple class from the OptimizationObject class, override two functions in your implementation of Apple (Contradict and Tension), and set the cost of leaving out the apple in the Apple::grade public member.
Conclusion
Best-fit analyses often require some sort of mathematical optimization operation. When that optimization involves selection and rejection of paths, dynamic programming is a likely candidate for the job.
You can find more details, as well as some simple test cases in the references below. I will be glad to exchange ideas and to hear about other applications. I can be reached at 72623.3367@compuserve.com.
References
[1] Mark G. Kerzner. "Image Processing in Well Log Analysis," IHRDC, Boston, 1986.
[2] Mark G. Kerzner. "Obtaining Structural and Stratigraphic Dip Information from Dipmeter Logs Using Segmentation Trees and Hierarchical Multilevel Optimization," 1986 Annual SPE Conference, New Orleans.
[3] Richard E. Bellman and Robert Kalaba. Dynamic Programming and Modern Control (Academic Press, 1966). ISBN 012-084856-2.
Mark Kerzner learned programming in high school on the first Russian computer, Ural 1. He has been programming for the oil industry for the last fifteen years. He is currently attempting to create the ultimate dipmeter program (again).