Limnol. Oceanogr., 26(2), 1981, 336-349 @ 1981, by the American Society of I,imnology
Immediate detection of heterogeneities in continuous multivariate, oceanographic recordings. Application to time series analysis of changes in the bay of Villefranche sur Mer
Station Zoologique, Abstract An index is described for immediate detection of discontinuities in continuously monitored, multivariate recordings. At every sampling point, the generalized distance is computed between the new multivariate observation and the centroid of the n previously recorded observations (D” to the centroid). These n points form a window, whose width determines the sensitivity of the index. By its mathematical formulation the distance takes into account the variations in range of the variables and their correlations. The validity of the index was first tested on simulated series into which different types of changes were arbitrarily introduced. Applying the index to a 4-year temperature and salinity series (weekly observations at the surface of a station in the bav of Villefranche sur Mer) leads to a coherent pattern of hydrological events. 06230 Villefranche sur Mer, France
The study of marine ecosystems rcquires simultaneous analysis of the changes in physical, chemical, and biological variables. Using continuous recorders (multiparameter probes, plotting in real time), we can observe the different stages of a dynamic system in time and space, each stage being defined by a set of chemical and biological measurements. These new techniques raise problems both for planning the sampling scheme and interpreting the complex data. On an oceanographic cruise, the planning is traditionally strict; it is developed a priori, and the coordinates of the stations are determined in advance. In classical design, both biological sampling and hydrological measurements yield only a few observations. Today, however, data collection is automated for many physical and chemical variables (temperature, salinity, pH, oxygen, nutrients, turbidity, light), and even for biological parameters like chlorophyll (in vivo fluorescence). Therefore, very large sets of observations are possible and data collection can respond to the oceanographic situation demonstrated by the data themselves. This means that cruise planning must be more H~exible in the short term, and it is necessary to develop a strategic
system to direct the sampling from moment to moment. As a means of dealing with this problem, I suggest a strategic approach based on rapid detection of discontinuities in the system of multivariate data that is being continuously recorded. For this purpose all the variables influencing sampling decision must be measured at comparable intervals. However, standardization of scales is far from being realized, since sampling interval is different for each variable: a few seconds (thermistors, salinometers, oxymeters, fluorometers), a few minutes (diffusiometers, particle counters, zooplankton pumping), or >1 h for vertical profiles with standard hydrological bottles. As soon as a heterogeneity is detected (hydrological front, high values of chlorophyll, patches of zooplankton, etc.) the sampling is modified: for instance the ship is stopped, and net tows or vertical hydrological profiles are made. Unfortunately, with the speed of accumulation of the records, it is not at all easy to recognize the exact point in time and space where hydrological changes occur. Moreover one can never be sure that high simultaneous fluctuations indicate the transition to another type of water mass or only an erratic phenomenon. I propose here an index that allows
ate stationarity is not usually fulfilled in reliable detection of the real discontinusitu. Cruzado (1971) used a technique ities in a continuous multivariate oceanbased on multiple correlations in succesographic data series. sive sets and his results are quite satisThe analysis for classical oceanographfactory, but his segmentation forbids a ic data applies well known statistical ‘precise determination of discontinuities. methods which assume the indepenInstead of taking successive sets of stadence of the observations. The methods tions, Webster (1973) used the generalused are variance, regression, and factor ized distance (D5 Mahalanobis 1936) beanalysis. With continuous monitoring, tween two groups of stations and then quantitative interpretation is quite differconsidered the distance between the two ent: the dependence between successive groups, separated from the preobservations in time and space can be following measured, and treatments in terms of pro- vious ones by only one unit of sampling. cesses are the most adequate. However, This method has been used (Ibanez spectral analysis, forecasting, and trans1976) to define the boundaries of homofer analysis always require a set of data geneous areas in zooplankton continuous with some constant properties, mainly series. Still, this technique, as will be stationarity. In spite of some data transseen further on, results in some artefacts formations (detrending, smoothing, dif- which are prejudicial to the interpretaferencing) it is necessary before a treat- tion. ment to find and distinguish the I present here an index that can be homogeneous area in the data to detect used in continuous sampling aboard ship, boundaries. We shall see how the index and its heuristicity is illustrated on a 4proposed can resolve that problem. year interpretative study of temperature A few workers have tried to formulate and salinity changes at one station in the the automatic identification of discontinbay of Villefranche sur Mer. uities in series analysis. Kelley (1971n) proposed the technique of principal com- Methods ponent analysis of successive groups of Identification of heterogeneities in stations: the distance between the first continuous accumulation of a multivariaxes of successive analysis is an indicaate data series can be achieved by using tion of the permanence of the change of one of four modes of comparisonenvironmental properties. However, the 1. Two groups of stations (two windows), case of two quite similar eigenvalues these two groups being successive. often leads to permutations among the 2. Two overlapping groups of stations, one higher order principal axes, and estimaunit of sampling apart, 3. Each new observation with a group of pretion of the distance may be strongly vious stations. biased. Also the correlations between the 4. The values predicted by a model with each variables and the principal components new observation. can change (not only in absolute value but also in sign) from one group to a folAn operator is defined: FkLo, where w lowing one and consequently alter the is the window width (number of stations), meaning of the distance. Kelley also sug- and k the position of the first station ingested that it is possible to use a multiple cluded in the window. For mode 1 the autoregressive model (Jones 1970). An windows FliZE and F,,..cwUJ are compared important difference between the model step-by-step, and for mode 2 the winexpectations and the data observed dows Fku’ and FI,.+,W are compared steppoints out a change in the environment. by-step. To obtain an index defining the This approach does not seem adequate in changes in a multivariate series, we must continuous analysis of oceanographic shift the windows along the series. If we data because the matrix computations are consider a shift equal to the window intricate, and the condition of multivariwidth (Kelley 1971b; Cruzado 1971), the
4 0 I Stations4’
30 1 Stat ions
of the four simulation
windows F,,“’ and F,,+,ll? and the windows F k blowand K- +zl,lu’ will be compared successively. When using mode 3 point i is compared to all the points of the window from i - x to i - 1. A sudden high value of the index from mode 1 means that a change has occurred from the group of stations included in the first window. For mode 2, this variation may result either from elimination of station k or for the k + 2w station input. Concerning mode 3, only one window moves along the series, so that an increase in the index results either from the w + 1 station input or from the removal of the first station of the preceding window. Simulation I will first compare the three previous models using simulated series. Forty observations were chosen at random for four simultaneously sampled variables (Fig. 1). The parameters of the four series are shown in Table 1 and we can see that the variables are Poisson-distributed, i.e. the variances about equal the means. To test the sensitivity of the index, I have placed very recognizable irregularities at particular spots in the multivariate series: an impulse and a local trend. An impulse is a sudden increase (or decrease) in the range of one or several variables which appears at one station,
Such a variation is larger than 2 SD. For example at station 20, I have simulated the very contrasted values, variable 1 = 1, variable 2 = 750, variable 3 = 500, variable 4 = 10. To test the influence of two impulses shifted by an interval smaller than the width of the window, I have simulated another impulse at station 25. Also I attempted to see how a gradient (increasing or decreasing) affecting one or several variables modifies the values of the index. In that case, four stations showing a local trend were placed between stations 19 and 20. The main results obtained with simulation models are presented according to the types of windows and the choice of different indices. The details of the simulations are reported only in the case of my proposed index. In order to compare two groups of observations, one can simply calculate the Euclidean distance between the centroids of the groups. This distance is quite independent of the dispersion of
Table 1. Distribution simulation variables. parameters of the four
0.66 1 0
10 Fig. 2. D’ to the ccntroid
20 on simulated series.
the variables and their correlations. With mode 1, the impulse was recorded by the index as soon as it appeared in the first window, and its effect was maintained as long it was part of either window. The local trend was expressed by a regular increase of the index corresponding to the progressive transition of the trend into the first window, then by a decrease when it gradually entered the second window, and by an increase again when it was located only in the last window. Considering the mode of two overlapped windows, I noticed a duplicated peak, corresponding to the impulse just as it goes into the second window. The local trend showed also duplication effects. I
conclude that this technique is not adequate to recognize the discontinuities in a multivariate series; the properties of the Euclidean distance are not precise enough and the use of two distinct windows produces artefacts. I have considered another type of index, to evaluate the changes in the correlation between the variables instead of the ranges. If we call the multiple correlation between the m variables for a group of stations R. 6 = 1,2, . . . , m) the total multiple correlation coefficient is R2r,,= + $ R2+
This index expresses the compactness
D’ to the ccntroid
D2 13.96 -
I 9.50~~-----------------~--~~----~~--------~~-------. --------.
A-.. 0 Fig. 4. 10 Impulse
effect at station 20 on the D2 to the ccntroid.
the group of variables. Rather than taking this value, I calculated an index corresponding to the difference of the successive R2T of two successive windows. With mode 1 of two windows shifted I could not identify the impulse, even by changing the width of the windows. If mode 2 of two overlapped windows was applied, I observed a duplication of the peaks. The definition of this index is not good in the sense that the changes of the signs of the Rj are not taken into account. Again the use of two windows is not adequate. To reconcile the effects of ranges and correlation variations, I have used the generalized distance of Mahalanobis between the centroids of two groups as proposed by Webster (1973):
D” = (x,
- x2)’ S-’ (x,
where X, and X2 are the vectors of means for each group of stations, S is the pooled variance-covariance matrix within groups. With mode 1, the impulse leads to a very weak increase corresponding to stations 21 and 22 instead of 20. Considering mode 2, I again find duplication of the peaks: the distance also increases when the impulse leaves the windows. None of the previous models gave precise, satisfying results, and consequently I have looked for an index and a type of window sliding, which avoids all these artefacts. It is based on the following principles: simultaneous comparisons of the ranges of the variables and of their
--10 Fig. 5. Effect of two impulses
>g 20 30 Stations 4o
at stations 20 and 25 on the D’ to the centroid.
of a local trend (from station 20 to 24) on the D” to the centroid.
correlations; sensitivity to the variations produced by the admission of a single new observation; convergence of the results using a relatively small window width. The first step is collecting a sequence of multivariate o’bservations, whose length is equal to the desired window length, Second, the centroid of that sequence is determined. At each subsequent step the distance is determined between a new observation and the centroid; heterogeneities are identified by large distance; the centroid is redetermined after deleting the earliest observation in the window and adding the new one; and so forth. The metric used is a particular form of the generalized distance, giving the distance between any observation and the centroid of an unique group (Dag&lie 1975): where x is the vector of the standardized coordinates (mean 0, variance 1) of the new point X, and R is the correlation matrix of the variables. Instead of using standardized variables, we can calculate an index with centered values and thus with the variance-covariance matrix. However, since the scales of variation differ from one parameter to another and we are mainly interested in profile variations, correlations only are considered. The basis of this index is simple: for an adequate window width, ccntroid posi-
tion will not significantly change as the window slides along the series. Heterogeneous or erratic values will not greatly alter the coordinates of the mean point, and progressive changes will gradually modify its position. On the other hand, a large distance from a given new point to the centroid can easily be recognized. The effect on centroid position caused by heterogeneous values leaving the window is suppressed by the large number of more typical values with which they are averaged. Our principal aim is to replace the multivariate series by a single series for visualization of the major heterogeneities. The profile of the values of D2 itself can be sufficient for detection, However, it is helpful to have an objective test of the significance of the values of D”. Provided that the distribution of the variables is multinormal, the distance can be tested by its equivalence to x2. As Cooley and Lohnes (1962) pointed out, the size of the hyperellipsoid of probabilities for a set of observations in a space of m dimensions is given by x2 = x’R--‘x with m df. It is obvious that this expression is identical to that for D”. If we again take the simulation example with 4 df, then if a value of D” is >9.48, the probability is <0.05 that the new observation X belongs to the group defined by the preceding stations. Although my data are very far from multinormality, the D’ of
series of surface temperature
at point B.
5% probability is plotted as threshold (dashed lines) on the later figures illustrating the properties of D”. This is particularly helpful for automatically scaled graphs such as these, where the ordinate scale varies from one series to another and can confuse the worker, especially at sea. Figure 2 shows the index obtained with the previous, simulated data using a window of 15 stations. The choice of the optimal window width is discussed later. There are variations in D2 along the series, but it remains below the criterion of 9.48. In particular there is a peak corresponding to station 33. The profile is the same with log-transformed data (Fig. 3). An impulse added at station 20 gives a significant peak, but the rest remain quite similar to the preceding examples (Fig. 4). Adding another impulse at station 25 shows the effect of two heterogeneities separated by an interval small-
er than the window width (Fig. 5). The succession of the two peaks occurs at the correct position. A local trend was simulated by including four supplementary stations between stations 19 and 20 over which the variables change progressively. The sequence of the D2 values (Fig. 6) indicates clearly the local increase starting from station 20. Time series The hydrological characteristics of station B in the bay of Villefranche sur Mer have been described by Bougis and Car& (1960), Bougis and Fenaux (1961), Bougis et al. (1965), and Bougis and Corre (1971). Here, I am only concerned with describing the hydrological changes given by the D2 to the centroid. I have considered a series of the hydrological 208 surface temperatures and salinities (weekly, representing about 4 years). The total depth of station B is about 80 m.
Figure 7 shows the temperature and salinity series from January 1968 to April 1972. I have computed first the index D” using a window of 52 observations. This width defines a centroid corresponding to yearly averages of temperature and salinity. Figure 8 shows the curve of the D2 obtained by sliding the window along the series. There are several peaks separated by rather homogeneous intervals of time. Representativity of the high D2 values-What do the peaks of the index mean with regard to the original data? First, a threshold has to be chosen in order to identify the most important values. Taking values of’ D” higher than 4, I consider the points located beyond the 10% probability ellipse (in the hypothetical case where binormality exists for the data, as discussed above). If these selected observations are transferred to a T/S diagram, we can verify that they effectively represent the remarkable values, those farthest from average temperature and salinity (Fig. 9). The D” peaks are effectively the marginal values of the distribution, and thus this index permits identification of the heterogeneities in a sequence. It is important to demonstrate that the calculated discontinuities are related to real hydrological changes. To do so, I have subdivided the high D2 values into four groups, each located in a quadrant; these quadrants have been defined by
the two lines passing through the mean point and are parallel to the axes of T/S. The months corresponding to the peripheral values of D2 in the four quadrants are recorded in Table 2. The significance of these peaks and the dates of their appearance are explained in a study of the hydrology of Villefranche sur Mer. Nival and Corre (1976) have identified three distinct periods of 4 months in the annual cycle. From April to July, there is great variability caused by the superposition of cyclical phenomena (caloric influxes or losses) with the events associated with geographical gradients (winds, currents, flow of the Roya and Var coastal rivers). From August to November, there is the beginning of homogeneity: this is the period when the thermocline and all vertical stability disappears. From December to March is the second period of homogeneity, characterized by cyclical movements on the vertical, caused by the cooling of superficial layers and also by abundant but erratic rainfalls. The peaks of D2 related to low temperature and high salinity correspond to the third period. However, since the vertical mixing processes are very important, we should expect close to mean values of salinity. To explain the high salinity I have considered the particular effects of westerly wind, which can raise subsuperficial waters (more saline). From meteorological recordings, it is clear that these events occur when there are winds
__ 153 38.06-. * . . . . ’ . . l . . . . : . . . . . l . : . . . . . . . l
‘. . . .
.. : . . . . . . . . . . l __--.-. . . . . .* : . . . . . . l
8 l .
a. + l
: . . . . . . . . . .
l l ’ .
. . . . .
. . ' '. . . 166 19 117 .
24 z . 23 .
. l . . . .
l * 158'
. 56 SC, . . 37.06 -197 204 I 15.41 Fig. 9. Position I T 63 160 162 64
20.41 25.41 To of high values of D” in T/S diagram,
~16 rn.sl, particularly two southern squalls between 17 and 21 m-s-’ in November 1970. The peaks of D’ related to high temperature and high salinity correspond to the end of the stage of preponderant heterogeneity, except for the values of October 1970 (however, weaker than for the other months). The peaks of D” here rcsuit from the temperature increase at the beginning of summer together with high values of salinity that are due to the lack of a feedback jet current in the bay (SaintGuily 1959), the thermocline being around 30 m deep. The peaks corresponding to warm water of low salinity are situated in the period of preponderant heterogeneity characterized by solar warming and abundant discharges of the coastal rivers
in spring. Rut we again find a particular month: September 1972. The months corresponding to low temperature and low salinity belong to the second period of homogeneity and the first part of heterogeneity. The low salinity can be explained by particularly
Table 2. Monthly vallles of D2, located diagram (see Fig. 9).
Low T/high S IIigh T/high
dates corresponding to high in the four quadrants of T/S
Nm70 Mar 71
Jun 68 Jul 69 Oct70
May Jun Jun Jul Sep
68 70 71 71 71
Feb Feb Mar Apr Jun Feb Mar
68 69 69 69 71 72 72
13.89 T 927 I
Weeks D2 12.82
of the D2 to the centroid,
heavy rainfalls in winter and in April and by snowmelt in May and June. Effect of the window width-The width of the window can be adjusted like the power of a microscope. To follow the fluctuations of my index, I have considered successive windows of 10, 20, 30, 40, 80, and 100 stations. Figure 10 shows the profiles of the index for 10, 30, and
100. The curve corresponding to a window of 10 observations is very unstable. The values are low and do not delimit particular zones. From 30 stations and on, the main peaks always occur at the same observation: 5, 6,23,24,60, 74, 117, 160, 170, 200, and 204. The zones between the peaks tend to become more and more monotonous as the width of the window
Table 3. Temperature and salinity variations for increases. So, starting from a relatively the centroid position of different window widths. small number of stations, the general profile of the curves of the index is not modMax T - Min T Max S - Min S ified, and the lesser discontinuities are 0.93 12.80 10 smoothed. The instability of D” is in9.44 0.60 20 versely related to the window width. Be6.13 0.44 30 cause this property results from the vari0.34 40 2.77 52 1.40 0.27 ations of the coordinates of the centroid 100 0.54 0.10 during the window translation, I have calculated the variations for six different widths of window (Table 3). The migration of the centroid decreases rapidly as then lagged successively one unit of samthe width increases. For a window of 52 pling forward, to the maximum T/2, where T is the total number of observastations, the ranged centroid position tions. For a lag 8, given represents 10% of the total temperature range and 16% of the total salinity range. Thus, there is a satisfying convergence of the results with increasing width of the window. where 8 = 0, 1, 2, . . . , T/2, and with the Selection of optimum window widthfirst term corresponding to the vector of Webster (1973) d iscussed the problem of the averages of a first group of observachoice a posteriori, after sampling, of the tions from t to T - 8 and the second term window width: if the window is too large, corresponding to the vector of the averdiscontinuities could be hidden because ages of a second group of observations of the presence of too many heterogenefrom 8+1 to T, we have ities in the window; if it is too small, too Auto D” (0) = 23’ S-l B. many discontinuities appear. It seems wise to consider a window equal or S-l is the inverse of the variance-covarislightly inferior to the mean distance be- ance within-groups matrix. tween the discontinuities of the multiTo eliminate the variance differences variate series. For this purpose, Webster between temperature and salinity, I first divided the data by the standard deviaadvised that the autocorrelation of each tion of each variable for each window. variable be examined and that the shift be defined corresponding, on average, to The matrix S becomes a correlation mathe values where the slopes of these cortrix. Figure 11 shows the Auto D” for the relations change. This shift, related to the surface temperature and salinity. There maximum drop of the autocorrelation, is an important change in the slope of the represents the mean distance between function after 8 = 57 weeks. The width the discontinuities. But if the variables of the window must be equal to this limit. are statistically very different, this techThus my preceding choice of 52 stations nique does not seem suitable, not to menwas suitable to show the main discontintion the time spent. In this cast, the minuities of the series. imum values of the autocorrelation will Although the Auto D’ may be helpful be very different and their average will in interpreting the data collected, it obnot be of interest. viously cannot be estimated on a ship. I have proposed (Ibancz 1976) a However, since very often the same variunique index, Auto Dz, which permits ables are sampled repeatedly, it can estimation of the scale of the changes of serve a priori to define suitable windows a multivariate series, The metric used is starting from previously recgrded data. the Mahalanobis D”. The calculation is Discussion and conclusion established as in the case of autocorrelation. The distance between the multiThe index that I propose seems to be variate series and itself is estimated and adequate for continuous sampling and for
of the f&ction
interpreting multivariate series. It remains to discuss some methodological or practical problems, The generalized distance describes at the same time the fluctuations of variable ranges and the changes in their correlations. In fact, Lefebvre (1976) showed by simulation that the D” is much more sensitive to the variations of the profiles than to the differences of absolute values. Moreover, there is a high similarity between the distance found by using the data and that found by using the log transforms. For continuous recording this aspect is quite attractive, because the variations in intensity for any variable are more influenced by merely random processes, which can be damped by transformations, than by the changes in the relationships between the variables. . As ascertained by simulation, a local trend or a succession of several impulses can be detected by the OZ. But we can ask whether a peak of D” indicates an erratic impulse in the original series or the transition to a new stable structure. By simulation (not presented here), I veri-
fied that when a boundary is like a hydrological front, with successive contrasted values, the D” value shows a lot of high individualized peaks. But if we imagine an instantaneous total change of the structure, the D’ shows at first a sudden increase, which is quickly weakened after two or three observations. Thus, the index reacts as in an impulse situation. Although this natural case of change is quite exceptional, the oceanographer has to refer to original records for the correct interpretation. On a ship the problem of the choice of window width is difficult because the best determination requires knowledge of all the data. Nevertheless, stabilization of the pattern of D” for a relatively small window has been demonstrated. Moreover, even if a very large width is chosen, the computation time remains short: using a small desk computer, with four variables and a window of a hundred stations, takes only 9 s for each D”. The interval of time between successive hydrological data sets is generally around 20 s, so that the computation can keep
of heterogeneities ->
pace even for a very substantial window width. I have not developed the possibilities of the D” for interpreting the data a posteriori, but I can propose many interesting treatments. For instance, to reduce the noise of the information, we can consider not the multivariate series itself but its first principal components. Also, since a series of D2 represents the time or space relationships of the variables, or both, it could be suitable for evaluating periodicities of the changes by spectral analysis. Going further, if we have two multivariate series, for example hydrological and meteorological, we can examine the instantaneous or the shifted correlations between both series of changes by using Q2. However, probably the main interest in this method concerns ecological problems. Classical marine ecological studies consist mainly of applying a statistical model: the estimation of the relationships between biological events and the intensity of one or several environmental variables (regression and multiregression analysis, canonical correlation, factor analysis, etc.). It seems attractive to be able to look also at the influence of changes in the environmental variables on biological processes, whether in real time or in shifted time. D2 is ideal for that, since it is a generalized measure of rate of change. References
BoUGIS, P., AND C. CARE&. 1960. Conditions
hydrologiques a Villefranchc sur mer pendant les annecs 1957 et 1958. Cah. Oceanogr. 12: 392408.
AND M. C. Conn~. 1971. Conditions hydrologiques k Villefranche sur mer pendant les an&es 1964, 1965, 1966 et 1967. Cah. Oceanogr. 23: 733-753. AND R. FENAUX. 1961. Conditions hydro-3 logiques a Villefranche sur mer pendant les ann&es 1959 et 1960. Cah. Oceanogr. 13: 627635. -, -, AND M. DEZILIERE. 1965. Conditions hydrologiques a Villefranchc sur mer pendant les an&es 1961, 1962 et 1963. Cah. Oceanogr. 17: 685-701. COOLEY, W. W., AND P. R. LOHNES. 1962. Multivariate procedures for the behavioral sciences. Wiley. CRUZADO, A. 1971. Sequential statistical analysis of continuous underway oceanographical data. Invest. Pesq. 35: 261-268. DAGN~LIE, P. 1975. Analyse statistique a plusieurs variables. Presses Agron. Gembloux. IRANEZ, F. 1976. Contribution k l’analyse mathdmatiquc des evcncments en ecologic planctonique. Bull. Inst. Oceanogr. Monaco 72, p. l-96. JONES, R. H. 1970. Change detection model for serially correlated multivariate data. Biometrics 26: 269-280. KEI,LEY, J. C. 1971a. A strategy for continuous multivariate analysis in oceanography. Invest. Pesq. 36: 175-179. -. 1971h. Factor analysis in real time. Invest. Pesq. 36: 179-183. LEFERVRE, J. 1976. Introduction aux analyses statistiques multidimensionnelles. Masson. MAHALANORIS, P. C. 1936. On the generalized distance in statistics Proc. Natl. Inst. Sci. India 12: 49-55. NIVAL, P., AND M. C. CO~RE. 1976. Variations annuelles des caracteristiques hydrologiques de surface dans la rade de Villefranche sur mer. Ann. Inst. Oceanogr. 52: 57-78. SAINT-GUILY, B. 1959. Measures de courant a l’ouvert de la baie de Villefranche sur mer. Cah. Oceanogr. 11: 602-604. WEBSTER, R. 1973. Automatic soil-boundary location from transect data J. Math. Geol. 5: 27-37.
Submitted: 11 December Accepted: 7 October