- Email: [email protected]
Error measures for template-fit geolocation based on light
ARTICLE IN PRESS
Deep-Sea Research II 54 (2007) 392–403 www.elsevier.com/locate/dsr2
Error measures for template-fit geolocation based on light Phil Ekstrom,1 Lotek Wireless, 114 Cabot Drive, St. John’s, NL, Canada A1C1Z8 Accepted 1 December 2006 Available online 20 February 2007
Abstract Template-fit geolocation offers an opportunity to estimate error in each individual day’s position results. The estimate requires determining an empirical calibration constant for latitude and another for longitude based on a trial against realistic noise in irradiance field data, which typically display strong point-to-point correlations. Such a calibration is reported here, making the error estimate available for use. The calibration is based on a 458-day set of blue-light irradiance data taken by an archival tag (Lotek LTD750) at a known location on land (48.571N, 122.941W). For the data studied, median values of the estimated daily standard deviations were 0.421 in longitude, 0.951 in latitude. The error estimates include the effects of weather and other data noise, latitude, and season, so (1) they vary widely from day to day and with season, and (2) they can be cautiously expected to remain valid in situations other than the particular one tested. A change in the character of the data noise is the only factor expected to have a significant effect on the calibration; this work should be repeated to sample the expected range of data noise types. r 2007 Elsevier Ltd. All rights reserved. Keywords: Errors; Error estimates; Geolocation; Position fixing; Statistical analysis
1. Introduction Some measure of accuracy is needed along with any estimate of geographic position if the position estimate is to be interpreted correctly as part of an animal track. While considerations of track continuity can be used to estimate random variations, some comparison with known positions is needed to provide an estimate of systematic effects and absolute accuracy. Ideally one seeks an error estimate that mirrors the effect of each day’s weather on that day’s position estimate. Tel.: +1 360 468 3375; fax: +1 360 468 3844.
E-mail address: [email protected]. Contact information Northwest Marine Technology, Inc., P.O. Box 427, Shaw Island, WA 98286, USA. 1
0967-0645/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.dsr2.2006.12.002
Previous studies of the threshold method of geolocation based on light, such as those by Musyl et al. (2001), have sought to characterise overall accuracy and to display seasonal variations in the accuracy of latitude determinations, but offered no way of associating an individual error estimate with each position estimate. However, it is clear that particular kinds of weather can distort the irradiance record so as to cause particular kinds of errors. There clearly are ‘‘bad days’’ for light-based geolocation, and conversely also ‘‘good days’’. It would be valuable to have our geolocation method tell us which was which. The more recently introduced template-fit method (Ekstrom, 2004) offers two error measures arising in each day’s calculation that are relevant to the accuracy of that day’s parameter estimates: the fit
ARTICLE IN PRESS P. Ekstrom / Deep-Sea Research II 54 (2007) 392–403
residual and the diagonal elements of the inverse Hessian matrix. They may be combined to form a daily estimate of standard deviation. This paper presents an overview of the templatefit method, examines template-fit geolocation results obtained from a 458-day set of irradiance data taken on land at a known location, considers the separate behaviour and significance of the two separate error measures, then proposes and tests against the data the combination to be used as an estimate of standard deviation. It obtains the necessary calibration factors for the particular noise present in the data record examined, which noise is primarily due to mid-latitude weather. 1.1. Symbols used
T a c L f R N p C, CLL C ff
EL ¼ Ef ¼
p p
ðC LL RÞ ðC ff RÞ
DL , Df
kL, kf
s sL ¼ kL E L sf ¼ kf E f
template function sun elevation angle ‘‘cloudiness’’ parameter longitude latitude residual sum of squares number of data points included in a day’s fit number of model parameters, in this case 4 inverse Hessian matrix, diagonal element for longitude, for latitude un-calibrated error measure for longitude un-calibrated error measure for latitude deviation of a position estimate from the known position scaling factors to be determined experimentally standard deviation standard deviation estimate for longitude standard deviation estimate for latitude
2. Materials and methods This section contains an overview of the templatefit method, a brief account of the data used, and
393
details peculiar to the specific implementation of the template-fit method as used in this study. 2.1. A summary of the template-fit method The template-fit method for geolocation begins with a template function for the variation of light with time during twilight, and fits that function to the actual record obtained on a particular day. There are four free model parameters: latitude, longitude, and ‘‘cloudiness’’ morning and evening. The best-fit latitude and longitude values are taken as the position estimate. The template is based on a simplified geophysical model developed in Ekstrom (2002). Light transmission through a clear atmosphere to the ground is dominated by the direct beam at midday, by singly scattered light (the same process that makes a clear sky blue instead of black) during ‘‘primary twilight’’, and by multiply scattered light during ‘‘secondary twilight’’. For blue light, the primary twilight regime extends between sun elevation angles of +31 and 51. Within that range and for sufficiently strong scattering (blue light), the analysis predicts for the variation of ground level irradiance vs. sun elevation angle (a) a universal function that can be approximated in closed form as T ¼ c expðu2 Þ=ð1 þ erfðuÞÞ; 2
up0,
¼ c expðu Þ=erfcðuÞ; u40, p where u ¼ ðR=2sÞ sin a. Here s ( ¼ 6.9 km) is the exponential scale height of the relevant part of the stratosphere and R ( ¼ 6378 km) is the radius of the earth, making u ¼ 21:5 sin a. Processes in the other two regimes do not lead to such rigidity of functional form, so the template is limited to the stated range in sun elevation where single-scattering dominates. Since the template is curved, it will match the shape of a data record only when correctly aligned in a, a fact that is crucial for geolocation. Numerical integration of Eq. (6) in Ekstrom (2002) yields an equally rigid and less approximate (but less compactly represented) template that is used for comparison with data in that work and for the fits in this work. The template derivation assumes a clear sky, but most of the scattering occurs in a stratospheric layer that is high enough to remain clear whatever the tropospheric weather. Any static configuration of tropospheric clouds below this scattering layer, or
ARTICLE IN PRESS P. Ekstrom / Deep-Sea Research II 54 (2007) 392–403
500
Log Irradiance, 32/decade
liquid water below that in which an animal may be swimming at constant depth (or for data corrected to constant depth), only attenuates light from the broad blue-sky source without significantly affecting the shape of T. Accordingly, the method absorbs the effect of static clouds along with instrument sensitivity and the attenuation due to any static water column into the value of the multiplicative factor c above and then identifies it as the ‘‘cloudiness’’ parameter, which can take on different values morning and evening. A second paper (Ekstrom, 2004) uses a standard astronomical coordinate transformation to rewrite sin a in terms of hour angle h, sun declination d, latitude f, and longitude L as sin a ¼ sin d sin f+ cos d cos h cos f. The hour angle h can in turn be written as h ¼ UT=4 L2E, where UT is the universal time of the observation in minute-of-day, L the longitude in degrees, and E the astronomical ‘‘equation of time’’. Both E and d depend only on the season. The result is a template for irradiance vs. time-ofday that depends on the season and has adjustable parameters latitude, longitude, and ‘‘cloudiness’’. The template is expressed on the same logarithmic scale as is used to represent the data and is fit to the data using the Levenberg–Marquardt method of nonlinear least squares. Morning and evening halftemplates are fit with common values for latitude and longitude, but separate values for cloudiness, leading to a four-parameter fit. The template can be seen in Fig. 1 best-fit to irradiance data from ‘‘good’’ and ‘‘bad’’ days for geolocation. The template is the two separate curved sections of heavy and (where separately visible) dashed lines. The ends of its morning and evening sections are connected in the figure by heavy (and where separately visible dashed) straight lines to the centre and edge of the figure to mark the +31 and 51 boundaries, and to make obvious the difference in morning and evening cloudiness values. This figure is discussed in more detail in Section 3.4. Derivatives of this now four-parameter template with respect to each parameter yield functions representing the data patterns corresponding to changes in that parameter. These are evaluated at each data point, and their products in all possible pairs are summed over all points within the range of validityPof the template to form the Hessian matrix, H jk ¼ i ½ðqT=qyj ÞðqT=qyk Þi ; where yj represents the jth one of the four parameters. The inverse of this matrix, often called the covariance matrix
400
300
200 0
200
400 Sample of Day
600
0
200
400 Sample of Day
600
500
Log Irradiance, 32/decade
394
400
300
200
Fig. 1. Time series irradiance data (light line) and best-fit templates (heavier lines, solid and dotted) for two successive days, Julian days 128 (top) and 129 (bottom). The horizontal lines connecting ends of the morning and evening template sections serve to guide the eye.
(though that is not always justified) and denoted by C, has a role in the conventional error analysis for least-squares fitting, and will play a comparable role in the estimates developed here. The inverse of the same Hessian matrix, suitably modified to implement a mixed fit strategy, is central to the Levenberg–Marquardt algorithm used in the fit process. Ekstrom (2004) makes the conventional suggestion that the relevant diagonal element of the inverse Hessian matrix be used as a measure of each parameter’s sensitivity to noise. It notes that these elements do represent the geometric illconditioning that is inherent in geolocation based on light when near the equator and/or the equinoxes, including the effects of extracting cloudiness
ARTICLE IN PRESS P. Ekstrom / Deep-Sea Research II 54 (2007) 392–403
parameters from the same data used for estimating location. Thus they can be used to extrapolate accuracy estimates calibrated at one location and one season to other locations and seasons. However, the noise in field irradiance data available for geolocation is dominated by strongly correlated processes such as cloud movement and imperfectly compensated attenuation changes due to diving, so that the rest of the conventional error treatment, which assumes statistically independent data, is inapplicable. This work pursues an extension to the case of correlated data. 2.2. Data-taking Blue-light irradiance data were taken at a dark rural site on Shaw Island, WA, USA, 48.571N, 122.941W, using a now-obsolete archival tag, Lotek model LTD750. A successor tag model having the same kind of fluorescent-conversion light detector and temperature-compensated clock is the Lotek LTD2310. The spectral acceptance band of the detector peaks at 470 nm with a sharp red-side edge near 500 nm, imposing the necessary restriction to blue light. Irradiance data were recorded as integers on a logarithmic scale having 15 points per decade. Dynamic range of the measurement included full noonday sun and extended downward far enough (nominally seven to nine decades depending on temperature) to always easily include the useful part of the twilight transient. That is, irradiance values were always available for all sun elevation angles well past 51 even in the darkest of cloudy weather. Beginning on 15 February 2000, data for 458 days were recorded in approximately 40-day segments, sometimes separated by short gaps, spanning 475 total days on the calendar. The recording interval was 128 s, yielding 675 points per day. 2.3. Geolocation processing After recovery from the tag, data were converted to an irradiance scale having 32 points per decade for compatibility with later model tags. However, quantisation noise in the data is that associated with the initial coarser scale of 15 points per decade. Data were analysed using the template-fit method as outlined in Section 2.1. After a manual determination of the start and end of each approximately 40-day group, processing proceeded automatically without further human intervention. Days were not otherwise selected, so the results can be expected to
395
represent fairly the weather and seasonal geometric difficulties of about 1.3 years of the natural environment on land, including three equinoxes. The method was not ‘‘tuned’’ in any way to the particular location, which can be regarded as a typical example of the natural environment on land far from city lights. Initial parameter values needed to help start the Levenberg–Marquardt nonlinear least-squares fit were obtained from an initial threshold-method analysis of the day’s data. A threshold corresponding to a sun elevation angle of zero, rather than a negative-angle threshold chosen for optimum accuracy was used in order to guarantee a unique threshold-method solution in all seasons. Two template fits were performed for each day, leading to what are referred to below as the ‘‘northern-start’’ and ‘‘southern-start’’ solutions. Both used an initial value for longitude equal to the threshold-method longitude estimate, but two different initial latitude values were derived from extreme values of the threshold latitude estimate. These extreme values were computed from day lengths that were, respectively, artificially lengthened and artificially shortened by 14 min of time (3.51 of hour angle), an estimated maximum error due to weather noise. When, near the equinoxes, the resulting uncertainty range in latitude included both hemispheres, one fit was started in each hemisphere half way between the equator and the extreme value estimate. When both fell in the same hemisphere, fits were started at the 14 and 34 points of the range spanned. The results for each day depend only on data gathered that day, with no carryover from previous solutions. Nor was any information about the known location available to the method, which was (except for the initial piecing together of the data record to deal with the obsolete tag’s memory limits) implemented so as to be capable of functioning autonomously at any geographic location. When the day’s fits are entirely successful, the results of the northern-start and southern-start are expected to coincide. It is a property of any geolocation method based entirely on blue-light measurements taken below +31 solar elevation that its results are largely independent of horizon detail. It is a further property of the template-fit method that its results are also largely independent of cloud cover so long as the degree of sky obscuration does not change significantly during the sunrise or sunset transient.
ARTICLE IN PRESS P. Ekstrom / Deep-Sea Research II 54 (2007) 392–403
396
3. Fit results This section presents a first graphically oriented look at the geolocation estimates and the accompanying diagnostic information, together with a qualitative discussion of several effects and relationships that are directly visible in the figures. 3.1. Angular quantities
a 100
50
b
3.2. Inverse Hessian-matrix elements Fig. 3 shows the behaviour of the inverse Hessianmatrix diagonal elements corresponding to the latitude and longitude parameters. The square root of each element is plotted, since that is directly commensurate with an error size and is what will be used below in the estimate of standard deviation. The element for longitude, trace a, shows little variation with season, while that for latitude, trace b, becomes very large just on the winter side of the equinoxes, the season where most of the difficulties occur that are visible in Fig. 2. 3.3. Number of points in the fit, N Fig. 4 shows the number of points in the fit, again northern-start shown solid, southern-start shown dotted, and the number theoretically expected for the known latitude and season shown as a smooth solid line. The data points included are those that lie between 51 and +31 sun elevation angle as judged by the method’s template, which can vary between iterations as the method’s current estimated value of
c 0
-50 100
200
300 JD2000
400
500
Fig. 2. a: Longitude and b: latitude results. In each case the northern-start results are shown by a solid line, the southern-start by a dotted line, and the true value by a solid horizontal line. c: Solar declination. Ordinate is ‘‘Julian Day 2000’’, time in days from noon, January 1, 2000.
Covariance element root
Latitude, Longitude, or Declination, Deg
Fig. 2 shows a composite view of several angular quantities with the intent of exhibiting general trends and time correlations rather than fine detail. The three traces collectively marked a show the longitude results with the northern-start trace drawn as a solid line and the southern-start drawn dotted, as will be the convention in most graphs that show both northern- and southern-start results. The known true location is a solid horizontal line. These three traces are barely distinguishable on this scale, and longitude accuracy does not depend visibly on season. The traces marked b are the comparable display for latitude. Trace c is the solar declination; zeros of this trace are the equinoxes. Note that the accuracy in latitude degrades severely near (especially just on the winter side of) the equinoxes, and a spurious southern hemisphere solution appears, which the southern-start solution follows as it migrates south during the first equinox shown, then north during the second. Since both it and the northern-start result represent actual solutions in the nearly-symmetric geometrical situation near the equinox, we cannot expect our diagnostic informa-
tion to help us choose between them, and will have to look for other criteria such as track continuity to exclude the inappropriate one. In two cases even the northern-start solution converges to the wrong hemisphere, a result of geometric ill-conditioning near the equinox combined with weather noise. We can expect our error measures to warn us about times when this kind of behaviour is likely, but not to specify how to resolve it. In field results, the situation would have to be resolved using such considerations as track continuity. Threshold methods exhibit a similar spurious solution unless the irradiance threshold is chosen to correspond to zero solar elevation, which choice is disadvantageous for accuracy.
0.4 b 0.2 a 0 100
200
300 JD2000
400
500
Fig. 3. a: Longitude and b: latitude results for the square root of the diagonal element of the inverse Hessian matrix.
ARTICLE IN PRESS P. Ekstrom / Deep-Sea Research II 54 (2007) 392–403
Root of Fit Residual per Point
80
N
60 40 20
15 10 5 0 100
100
200
300 JD2000
400
500
Fig. 4. Number of data points in the fit. Solid line, northernstart; dotted line, southern-start; smooth solid line, theoretically expected number for the known latitude.
the latitude parameter converges toward its final value. The value shown is for the last iteration. Since the sun sets most rapidly near the equinoxes, there are fewer such points then, and more near the solstices when the sun sets more slowly. Note that near the equinoxes, the relatively large errors in the tag’s estimation of latitude result in errors in the assumed rate of sunset and consequently cause more or fewer than the correct number of points to be included. This can have an indirect further effect on accuracy by including data from outside the 51 to +31 ‘‘sweet spot’’ that occurs during twilight. 3.4. Residual sum of squares The residual sum of squares (Fig. 5) exhibits no pronounced variation with season, but rather varies widely from day to day, reflecting the effects of weather. Since the number of points in the fit varies with the season, and since for a given degree of noisiness in the data the size of the residual is expected to be proportional to the number of points, it is the residual per point corrected for the number of parameters estimated, R/(Np), that is displayed here. The graph shows the square root of this quantity, since that is commensurate with a standard deviation. We can see part of what this quantity is trying to tell us by examining a superposition of template and data for two successive days that have very different residuals. Fig. 1 (upper panel) shows Julian day 128 (JD2000, referenced to the millennium boundary), which has indistinguishable template traces (heavy, solid, and dotted) that fit the data well enough to lie atop each other and entirely obscure the data (lighter line) in the fit region between 51 and +31 sun elevation. The root residual per point is the same for both fits and is only 0.91 light units, very little larger than the lower limit of 0.62 set by the
397
200
300 JD2000
400
500
Fig. 5. Square root of the residual sum of squares per point. Solid trace, northern-start; dotted trace, southern-start. The solid horizontal line is the theoretically expected minimum residual of 0.62 LSB/point due to the original coarse 15/decade resolution of the recorded data.
resolution to which the irradiance was originally recorded. Fig. 1 (lower panel) shows the corresponding record for the following day. Here the data trace is significantly distorted from the smooth pattern typical of a clear day, presumably by clouds that erratically changed their effective fractional sky cover and/or opacity during the day. They also did so during the sunrise transient, which caused the difficulty we see. Beginning with the two different starts the two attempts converged to different local minima in the residual, neither representing a very satisfactory fit to the data record. During sunrise that record does in fact not look much like anything within the range of shapes reachable by stretching, translating, or raising the template in response to parameter changes. The root residual per point in one case was 12.4 units, in the other 11.4 units on a logarithmic scale where 32 units equal a factor of 10 irradiance change. 4. Error measures This section explores in more theoretical detail the properties of the diagnostic information available from the daily fits, proposes the daily standard deviation estimate, and explores its performance relative to the actual observed errors in position estimation. 4.1. The inverse Hessian-matrix diagonal element The error analysis usually associated with leastsquares fitting (Press et al., 1988) assumes that noise in the data points being fit is equivariant and statistically independent. When this hypothesis is satisfied, then one may obtain an estimate of the variance in a parameter estimate as the product of the diagonal element of the inverse Hessian matrix
ARTICLE IN PRESS 398
P. Ekstrom / Deep-Sea Research II 54 (2007) 392–403
(called in that reference the covariance matrix) associated with the parameter (either CLL or C ff ) times the common variance of all the (equivariant) data points. The analysis can be readily extended to the case where the data points have individual known variances, but an extension to the highly correlated errors in real field irradiance data is not available. Still, as there is reason to expect theoretically and also empirically as we have seen illustrated in the figures above, that the inverse Hessian matrix does capture the inherent error-sensitivity of the estimate that is being made, and in particular the dependence of error sensitivity of the latitude estimate on season. Theoretically it also should capture the effects of latitude itself on the error sensitivity of both estimates.
4.2. Comparison with threshold method error sensitivity In contrast to the template-fit method, a threshold method first estimates day length as the time between the sunrise and sunset crossings of an irradiance threshold. Then it chooses as its latitude estimate the value that gives the observed day length at the known season between the occurrences of the sun elevation angle that is taken to correspond to the adopted threshold. This involves solving for the latitude f the same transformation equation we saw in Section 2.1: sin a ¼ sin d sin f+cos d cos h cos f, where a is the assumed sun elevation angle, h the hour angle of the sun at sunset, and d the sun declination. In this context the quantity 2h is the day length expressed as an angle (11 of angle ¼ 4 min of time). The sun declination can be approximated as sin d ¼ sin(280.461+3601n/365. 25) sin(23.441), where n is JD2000, the Julian Day referred to the start of year 2000. Unless the sunrise–sunset sun elevation (a) is exactly zero (geometric sunset), no closed-form solution is available and this equation must be solved iteratively. However, one can obtain a closed form for the derivative of implied latitude with respect to day length that is valid for any value of a: df/dh ¼ cos d cos2 f sin h/(sin dsin a sin f). When this derivative becomes large, a simple linear error model, Df ¼ df=dh Dh; predicts large errors in latitude caused by small errors in day length. In particular, the derivative is infinite when sin d ¼ sin a sin f.
This measure reflects the same geometric situation that degrades the accuracy of the template-fit method: near the equinox, twilight details do not depend strongly on latitude, so little latitude information can be extracted from them. Accordingly it is not surprising that this measure behaves qualitatively much like the inverse Hessian root. The match between them is closest for a threshold set near 5.51. Fig. 6 plots the two functions with the threshold method sensitivity multiplied by an arbitrary constant that makes the two agree near the summer solstice. The light solid trace is the corresponding threshold method sensitivity for a threshold of +31, multiplied by the same arbitrary factor. The graph illustrates that the error sensitivities to angular errors depend on the threshold chosen. These sensitivities exhibit their largest values at different sun declinations, and therefore at slightly different times of the year, a fact first pointed out by Hill and Braun (2001). This suggests the possibility of using the threshold method with a high threshold to fill in during the seasons when the template-fit method is at its worst. Naively assuming that a threshold method is as good as a template fit during the summer (which it appears not to be; see rough indications in Ekstrom, 2004) and that the day length accuracy attainable with a +31 threshold is as good as with one at 5.51, all of which is suggested visually by traces a through c of Fig. 6, it would appear advantageous to use the threshold method with a 31 threshold not only during the equinoxes, but for the entire interval from Julian day 270 through 440. Unfortunately the rate of change of irradiance with time is about a factor of 5.4 less rapid at +31 sun elevation angle as it is at 5.51, so comparable errors in irradiance can be expected to lead to errors in day length that are larger by that factor of 5.4 at +31 than they are at 5.51. To illustrate the impact of this effect, the light dashed trace of Fig. 6 is a copy of the +31 angular sensitivity trace, but multiplied by 5.4 for fair comparison with the 5.51 trace assuming comparable errors in determining irradiance. Even assuming that the absolute errors in a 5.51 threshold method are no worse than those of the template-fit method (which is visually suggested by the choice of arbitrary multiplier, but not asserted as true), when the difference in irradiance slope is taken into account the expected latitude error of a +31 threshold method is always substantially greater than either alternate method except very near the equinox, where none of
ARTICLE IN PRESS
Latitude SensitivityMeasure
P. Ekstrom / Deep-Sea Research II 54 (2007) 392–403
399
0.4
0.3
d d
0.2 c 0.1
a,b
c
0 50
100
150
200
250
300 JD2000
350
400
450
500
1 Fig. 6. The inverse Hessian-matrix root for latitude (a, heavy solid trace, identical to Fig. 2 trace b), compared with 20 times the latitude sensitivity measure for 5.51 threshold method (b, heavy dashed) and +31 threshold method (c, light solid trace) computed for the known true latitude. Trace d is a copy of c multiplied by a factor of 5.4 to account for the different irradiance slope at +31.
the three are useful, and therefore offers no advantage. 4.3. Proposed estimates of template-fit geolocation error In order to make use of the inverse Hessian elements in estimating geolocation error, it is necessary to have a measure of data noisiness to multiply by. The situation does not permit the usual first choice, the independently determined variance of uncorrelated, equivariant data, because the available data is at minimum strongly correlated. The most obvious (and traditional) candidate for replacing that first choice is some multiple of the residual sum of squares. As we have seen, it does measure the deviation of the data record from the form predicted by the template, and thus is expected to be large on a ‘‘bad day’’ when the assumptions underlying that template, notably that absorption due to clouds or water absorption does not change during twilight, are not realised. Making that choicep leads first to the un-calibrated p error measures, E L ¼ ðC LL RÞ and E f ¼ ðC ff RÞ. It would be more conventional to divide the residual sum of squares by the number of points in the fit, less the number of parameters estimated. That would be an appropriate choice if the dominant noise source were uncorrelated. In the present case, the noise is ordinarily over-sampled in time and strongly correlated in the record. These uncalibrated estimates will be eventually multiplied by two empirically determined constants, kL and kf, respectively, to obtain
estimates of the standard deviation to associate with the day’s position estimate. That is, the eventual calibrated estimates will be p sL ¼ k L E L ðC LL RÞ and and sf ¼ k E , leading to s ¼ k f f L L p sf ¼ kf ðC ff RÞ. 4.4. Avoiding spurious solutions As noted above, the erroneous southern-hemisphere solutions are real solutions that arise from the geometric near symmetry present near the equinox when the sun lies above the equator. Thus the method provides no diagnostic information that can serve to automatically choose the correct one of two solutions that disagree. However, while the method may be in doubt as to which hemisphere the animal is in, the experimenter is rarely in doubt. Track continuity can ordinarily resolve the question in situations where the issue arises at all. In line with this, the development from here on will make use of our comparable knowledge that the test site is in the northern hemisphere, and will do something that would be available to an experimenter in a real field situation. It will concentrate on the northern-start solutions in the rest of this paper, and will also reject as invalid the two northern-start solutions that ended up in the southern hemisphere during the time of equinoctial ill-conditioning. 4.5. Bias in the position results Given an estimate for the individual standard deviations of each day’s position estimates, the
ARTICLE IN PRESS P. Ekstrom / Deep-Sea Research II 54 (2007) 392–403
optimum estimate for the mean is an average weighted by the inverses of those standard deviations. Thus the observed systematic biases in the position estimates can be obtained in terms of the deviations D from the known position as BL ¼ SðDL =sL Þ=Sð1=sL Þ and Bf ¼ SðDf =sf Þ=Sð1=sf Þ. These sums can be computed now, while the calibration factors are still unknown because those factors cancel in the final division. That is, only the un-calibrated product error measure is needed. The results for the northern-start data are BL ¼ 0:041 and Bf ¼ 0:151. The standard deviations in those p weighted means are sL ¼ ðS1Þ=Sð1=sL Þ and p sf ¼ ðS1Þ=Sð1=sf Þ, where (S1) is used to represent the number of points in each of the sums. However, these quantities cannot be evaluated without values for the calibration factors, which are computed next.
400
Rank order
400
200
0 -40
-20 0 20 Normalized longitude error
40
Fig. 7. Two superimposed and barely distinguishable plots regarding longitude results: rank order plot of southern-start values of ðDL BL Þ=E L (solid), and best-fit scaled Gaussian CDF (dashed), which has s ¼ kL ¼ 0:90.
4.6. Determining kL and kf
Rank order
400
Once the proposed estimates of data variability are correctly calibrated, the normalised, bias-corrected quantity ðDL BL Þ=sL will be a random variable with zero mean and unit standard deviation. Multiplying through by the constant kL one may conclude that the quantity ðDL BL Þ=E L will have mean zero and standard deviation kL. Similarly, the quantity ðDf Bf Þ=E f will have mean zero and standard deviation kf. Thus one can determine the values of kL and kf by examining those distributions. The solid line in Fig. 7 is a rank order plot, that is, an empirical un-normalised Cumulative Distribution Function (CDF) of the first of these quantities. The superimposed dashed line is the theoretical CDF of a normal distribution (multiplied by the total number of data points to match vertical scaling) that is a best fit to it. While detailed examination of the tails shows some divergence between them, the two distributions can be barely distinguished in the figure. The best-fit value of standard deviation is 0.90, which one may identify with kL. Fig. 8 is a similar plot of the normalised latitude errors, with two outliers that lie in the southern hemisphere removed from the set, as explained in Section 4.4. Again there is a superimposed Gaussian, best-fit with a standard deviation of kf ¼ 1:17. Again the cumulative distributions largely coincide, but in this case there is a much longer negative tail on the data CDF involving about 3 of the remaining 456 points.
200
0 -40
-20 0 20 Normalized latitude error
40
Fig. 8. Rank order plot of southern-start values of ðDf Bf Þ=E f (solid), and best-fit scaled Gaussian CDF (dashed), which has s ¼ kf ¼ 1:17.
Inserting the calibration values just determined, p the two noise p estimates become sL ¼ 0:90 ðC LL RÞ and sf ¼ 1:17 ðC ff RÞ. As will be argued shortly, the following information should accompany any such calibration: the data used were taken at a recording interval of 128 s on land at mid-latitudes (Shaw Island as described) in approximately 40-day continuous blocks providing usable records for 458 of 475 total days. The missing days were determined by the accidents of data collection and did not represent any deliberate selection. Now with the scaling constants in hand, it is possible to evaluate the expected standard deviation in the weighted mean used to calculate the bias values, and to state those values with their one-s
ARTICLE IN PRESS P. Ekstrom / Deep-Sea Research II 54 (2007) 392–403
401
Standard deviation estimate
100
10
1
0.1 50
100
150
200
250
300 JD2000
350
400
450
500
Fig. 9. The daily standard deviation estimates for longitude (squares, solid line is its 0.421 median) and latitude (crosses, dashed line is its 0.951 median). Two latitude standard deviation estimates of 1221 and 2421 occurring near Julian day 430 lie off the graph.
limits as BL ¼ 0:041 0:021 and Bf ¼ 0:151 0:041. The template-fit method appears empirically at least to be reasonably unbiased. The values the daily position error estimates take on for this data are plotted in Fig. 9. The error estimates for latitude vary widely from values similar to or even slightly better than those for longitude near the summer solstice to very large values near the equinoxes. For this data the median value of sL is 0.421, the median value of sf is 0.951. 5. Discussion The development above yields an estimate for each day’s uncertainty in template-fit geolocation results, based on data obtained at a particular midlatitude dark site on land in all seasons and all weather. The method used is applicable to the completely automated extraction of geographic position from field data on the time sequence of blue-light irradiance. It remains to consider how broadly applicable these error measures may be to locations and conditions other than those explicitly represented by the data examined. The method is based on a geophysical model of sunset that relates sun elevation angle to observable irradiance. For sufficiently strong scattering (blue light) the shape of the resulting template function depends only on the ratio of two independently measurable quantities, the stratospheric scale height and the radius of the Earth, and it has no available adjustments. There is no need to ‘‘tune’’ the template to any particular location, and no obvious way to do so. The template shape is very robustly determined so long as the scattering is strong
enough. For example, a significant increase in stratospheric aerosol content, perhaps due to a volcanic eruption, will move the scattering peak to higher altitude and change the value of the ‘‘cloudiness’’ parameter, but will not alter the shape of the template function. The next step in developing the template uses an exact astronomical coordinate transformation to rewrite the sun elevation in terms of latitude, longitude, season, and time of day. As a result the place where the data happened to be measured can be considered as an entirely typical mid-latitude site. Since the effect of a change in longitude on the data record is entirely equivalent to resetting the clock, extension to any other longitude is immediate. Both latitude and season affect the amount of information that there is in the data record regarding latitude or, from another viewpoint, the degree of geometric ill-conditioning that occurs. This has a major effect on the accuracy with which the template parameters can be determined from the data, but that error sensitivity is reflected in the inverse Hessian-matrix elements which form part of the proposed error estimators. The fact that these quantities successfully represent the variation in accuracy with season has been examined and verified in this work, and the variation with latitude is just another consequence of the same coordinate transformation. It is reasonable to expect, though this work has not separately verified it, that the error measures studied here will adequately account for the geometric effects of site latitude as well. The most heuristic aspect of this work is the use of the residual sum of squares as a measure of data noise, and that is also the aspect whose extrapolation
ARTICLE IN PRESS 402
P. Ekstrom / Deep-Sea Research II 54 (2007) 392–403
to other situations is the least secure. In an important sense this quantity cannot directly see the noise we care about. If a pattern of noise in the data is such that it can be approximated successfully by a parameter change, that parameter change will occur when the template is fit to the data, causing an error in the parameter value. That noise pattern will be absorbed, and will not appear in the residual. If the change is to the latitude or longitude parameter, it will cause a geolocation error. What is left in the residual is the effect of those noise patterns that cannot be absorbed by parameter changes, and that therefore have had no effect on the geolocation result. When using the residual sum of squares as a measure of noise for estimating parameter errors, one necessarily assumes that total noisiness in the data is on the average partitioned in some approximately fixed ratio between the harmless data patterns that do affect its value and the harmful ones that do not. That assumption is implied by using one as a measure of the other. If the character of the noise changes, the ratio of harmful to harmless noise components in each day’s data may change, and may affect the appropriate value of the calibration constants, kL and kf. As a result, weather with very different statistical behaviour, or imperfectly compensated diving noise as opposed to weather noise, may call for a different value of calibration constants even though both are examples of highly correlated noise. Therefore one should regard these present k values with greater suspicion, the less similar the natural noise processes in a new situation are to weather at the test site where the data used here were taken. We can already see a hint of this effect in the different sizes of the two k constants; the noise encountered here had a greater effect on the latitude estimate than on the longitude. Since all other effects would seem to be accounted for in other ways, the only obvious reason for the difference is an accident that the particular data patterns most harmful to latitude determination were more prevalent in the weather noise encountered than were patterns harmful to the determination in longitude. Therefore the calibration parameters should be re-determined at every opportunity where there is reason to suspect a difference in noise properties. They should also be re-determined for any new data interval, and any determination of k should always report the data interval along with anything that might distinguish the noise processes involved. Doing so in a variety of situations will explore the
range of variation of the calibration parameters, along with the nature of the situations in which each value is appropriate. Until this exploration has provided a broad understanding, the values reported here can be used with caution to begin providing a way of distinguishing good days from bad days, and of forming some idea of how good and how bad they are, in field experiments on animals. 6. Conclusion Error estimates for each daily position are available as part of the template-fit method of light-based geolocation. These estimates are sL ¼ p p 0:90 ðC LL RÞ and sf ¼ 1:17 ðC ff RÞ, where the symbols are as defined in Section 1.1. The numerical proportionality constants given have been determined for the particular case of observations made at a mid-latitude site (Shaw Island, WA, USA) on land every 128 s in all seasons. For those data, the estimates had median values of 0.421 in longitude, 0.951 in latitude. It is desirable to re-determine the calibration constants in any situation where the noise correlations are suspected to differ significantly from those of the data used here. Acknowledgements Thanks are due to Raymond Glaze for accumulating the long-term irradiance records used in this study, to Hannah Halliday for helping scan the data for interesting features, to Steven Teo for suggesting an extension of the original poster version, and to a very helpful anonymous reviewer.
References Ekstrom, P.A., 2002. Blue twilight in a simple atmosphere. In: Shaw, J.A. (Ed.), Proceedings of SPIE Annual Meeting 2002: Atmospheric Radiation Measurements and Applications in Climate, vol. 4815, Paper 14. SPIE, Bellingham, WA, USA, pp. 73–81. Available at /http://www.lotek.com/library.htmS [doi:10.1117/12.482305]. Ekstrom, P.A., 2004. An advance in geolocation by light. In: Yasuhiko, N. (Ed.), Proceedings of the International Symposium on Bio-logging Science. Memoirs of the National Institute of Polar Research, Special Issue, 58. National Institute of Polar Research, Tokyo, Japan, pp. 210–226. Available at /http://polaris.isc.nipr.ac.jp/ penguin/oogataHP/pdfarticles/22p210-226.pdfS. Hill, R.D., Braun, M.J., 2001. Geolocation by light levels—the next step: latitude. In: Sibert, J.R., Nielsen, J.L. (Eds.),
ARTICLE IN PRESS P. Ekstrom / Deep-Sea Research II 54 (2007) 392–403 Electronic Tagging and Tracking in Marine Fisheries. Kluwer, Dordrecht, The Netherlands, pp. 315–330. Musyl, M.K., Brill, R.W., Curran, D.S., Gunn, J.S., Hartog, J.R., Hill, R.D., Welch, D.W., Eveson, J.P., Boggs, C.H., Brainard, R.E., 2001. Ability of archival tags to provide estimates of geographical position based on light intensity. In:
403
Sibert, J.R., Nielsen, J.L. (Eds.), Electronic Tagging and Tracking in Marine Fisheries. Kluwer, Dordrecht, The Netherlands, pp. 343–368. Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T., 1988. Numerical Recipes in C. Cambridge University Press, Cambridge, UK.
Deep-Sea Research II 54 (2007) 392–403 www.elsevier.com/locate/dsr2
Error measures for template-fit geolocation based on light Phil Ekstrom,1 Lotek Wireless, 114 Cabot Drive, St. John’s, NL, Canada A1C1Z8 Accepted 1 December 2006 Available online 20 February 2007
Abstract Template-fit geolocation offers an opportunity to estimate error in each individual day’s position results. The estimate requires determining an empirical calibration constant for latitude and another for longitude based on a trial against realistic noise in irradiance field data, which typically display strong point-to-point correlations. Such a calibration is reported here, making the error estimate available for use. The calibration is based on a 458-day set of blue-light irradiance data taken by an archival tag (Lotek LTD750) at a known location on land (48.571N, 122.941W). For the data studied, median values of the estimated daily standard deviations were 0.421 in longitude, 0.951 in latitude. The error estimates include the effects of weather and other data noise, latitude, and season, so (1) they vary widely from day to day and with season, and (2) they can be cautiously expected to remain valid in situations other than the particular one tested. A change in the character of the data noise is the only factor expected to have a significant effect on the calibration; this work should be repeated to sample the expected range of data noise types. r 2007 Elsevier Ltd. All rights reserved. Keywords: Errors; Error estimates; Geolocation; Position fixing; Statistical analysis
1. Introduction Some measure of accuracy is needed along with any estimate of geographic position if the position estimate is to be interpreted correctly as part of an animal track. While considerations of track continuity can be used to estimate random variations, some comparison with known positions is needed to provide an estimate of systematic effects and absolute accuracy. Ideally one seeks an error estimate that mirrors the effect of each day’s weather on that day’s position estimate. Tel.: +1 360 468 3375; fax: +1 360 468 3844.
E-mail address: [email protected]. Contact information Northwest Marine Technology, Inc., P.O. Box 427, Shaw Island, WA 98286, USA. 1
0967-0645/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.dsr2.2006.12.002
Previous studies of the threshold method of geolocation based on light, such as those by Musyl et al. (2001), have sought to characterise overall accuracy and to display seasonal variations in the accuracy of latitude determinations, but offered no way of associating an individual error estimate with each position estimate. However, it is clear that particular kinds of weather can distort the irradiance record so as to cause particular kinds of errors. There clearly are ‘‘bad days’’ for light-based geolocation, and conversely also ‘‘good days’’. It would be valuable to have our geolocation method tell us which was which. The more recently introduced template-fit method (Ekstrom, 2004) offers two error measures arising in each day’s calculation that are relevant to the accuracy of that day’s parameter estimates: the fit
ARTICLE IN PRESS P. Ekstrom / Deep-Sea Research II 54 (2007) 392–403
residual and the diagonal elements of the inverse Hessian matrix. They may be combined to form a daily estimate of standard deviation. This paper presents an overview of the templatefit method, examines template-fit geolocation results obtained from a 458-day set of irradiance data taken on land at a known location, considers the separate behaviour and significance of the two separate error measures, then proposes and tests against the data the combination to be used as an estimate of standard deviation. It obtains the necessary calibration factors for the particular noise present in the data record examined, which noise is primarily due to mid-latitude weather. 1.1. Symbols used
T a c L f R N p C, CLL C ff
EL ¼ Ef ¼
p p
ðC LL RÞ ðC ff RÞ
DL , Df
kL, kf
s sL ¼ kL E L sf ¼ kf E f
template function sun elevation angle ‘‘cloudiness’’ parameter longitude latitude residual sum of squares number of data points included in a day’s fit number of model parameters, in this case 4 inverse Hessian matrix, diagonal element for longitude, for latitude un-calibrated error measure for longitude un-calibrated error measure for latitude deviation of a position estimate from the known position scaling factors to be determined experimentally standard deviation standard deviation estimate for longitude standard deviation estimate for latitude
2. Materials and methods This section contains an overview of the templatefit method, a brief account of the data used, and
393
details peculiar to the specific implementation of the template-fit method as used in this study. 2.1. A summary of the template-fit method The template-fit method for geolocation begins with a template function for the variation of light with time during twilight, and fits that function to the actual record obtained on a particular day. There are four free model parameters: latitude, longitude, and ‘‘cloudiness’’ morning and evening. The best-fit latitude and longitude values are taken as the position estimate. The template is based on a simplified geophysical model developed in Ekstrom (2002). Light transmission through a clear atmosphere to the ground is dominated by the direct beam at midday, by singly scattered light (the same process that makes a clear sky blue instead of black) during ‘‘primary twilight’’, and by multiply scattered light during ‘‘secondary twilight’’. For blue light, the primary twilight regime extends between sun elevation angles of +31 and 51. Within that range and for sufficiently strong scattering (blue light), the analysis predicts for the variation of ground level irradiance vs. sun elevation angle (a) a universal function that can be approximated in closed form as T ¼ c expðu2 Þ=ð1 þ erfðuÞÞ; 2
up0,
¼ c expðu Þ=erfcðuÞ; u40, p where u ¼ ðR=2sÞ sin a. Here s ( ¼ 6.9 km) is the exponential scale height of the relevant part of the stratosphere and R ( ¼ 6378 km) is the radius of the earth, making u ¼ 21:5 sin a. Processes in the other two regimes do not lead to such rigidity of functional form, so the template is limited to the stated range in sun elevation where single-scattering dominates. Since the template is curved, it will match the shape of a data record only when correctly aligned in a, a fact that is crucial for geolocation. Numerical integration of Eq. (6) in Ekstrom (2002) yields an equally rigid and less approximate (but less compactly represented) template that is used for comparison with data in that work and for the fits in this work. The template derivation assumes a clear sky, but most of the scattering occurs in a stratospheric layer that is high enough to remain clear whatever the tropospheric weather. Any static configuration of tropospheric clouds below this scattering layer, or
ARTICLE IN PRESS P. Ekstrom / Deep-Sea Research II 54 (2007) 392–403
500
Log Irradiance, 32/decade
liquid water below that in which an animal may be swimming at constant depth (or for data corrected to constant depth), only attenuates light from the broad blue-sky source without significantly affecting the shape of T. Accordingly, the method absorbs the effect of static clouds along with instrument sensitivity and the attenuation due to any static water column into the value of the multiplicative factor c above and then identifies it as the ‘‘cloudiness’’ parameter, which can take on different values morning and evening. A second paper (Ekstrom, 2004) uses a standard astronomical coordinate transformation to rewrite sin a in terms of hour angle h, sun declination d, latitude f, and longitude L as sin a ¼ sin d sin f+ cos d cos h cos f. The hour angle h can in turn be written as h ¼ UT=4 L2E, where UT is the universal time of the observation in minute-of-day, L the longitude in degrees, and E the astronomical ‘‘equation of time’’. Both E and d depend only on the season. The result is a template for irradiance vs. time-ofday that depends on the season and has adjustable parameters latitude, longitude, and ‘‘cloudiness’’. The template is expressed on the same logarithmic scale as is used to represent the data and is fit to the data using the Levenberg–Marquardt method of nonlinear least squares. Morning and evening halftemplates are fit with common values for latitude and longitude, but separate values for cloudiness, leading to a four-parameter fit. The template can be seen in Fig. 1 best-fit to irradiance data from ‘‘good’’ and ‘‘bad’’ days for geolocation. The template is the two separate curved sections of heavy and (where separately visible) dashed lines. The ends of its morning and evening sections are connected in the figure by heavy (and where separately visible dashed) straight lines to the centre and edge of the figure to mark the +31 and 51 boundaries, and to make obvious the difference in morning and evening cloudiness values. This figure is discussed in more detail in Section 3.4. Derivatives of this now four-parameter template with respect to each parameter yield functions representing the data patterns corresponding to changes in that parameter. These are evaluated at each data point, and their products in all possible pairs are summed over all points within the range of validityPof the template to form the Hessian matrix, H jk ¼ i ½ðqT=qyj ÞðqT=qyk Þi ; where yj represents the jth one of the four parameters. The inverse of this matrix, often called the covariance matrix
400
300
200 0
200
400 Sample of Day
600
0
200
400 Sample of Day
600
500
Log Irradiance, 32/decade
394
400
300
200
Fig. 1. Time series irradiance data (light line) and best-fit templates (heavier lines, solid and dotted) for two successive days, Julian days 128 (top) and 129 (bottom). The horizontal lines connecting ends of the morning and evening template sections serve to guide the eye.
(though that is not always justified) and denoted by C, has a role in the conventional error analysis for least-squares fitting, and will play a comparable role in the estimates developed here. The inverse of the same Hessian matrix, suitably modified to implement a mixed fit strategy, is central to the Levenberg–Marquardt algorithm used in the fit process. Ekstrom (2004) makes the conventional suggestion that the relevant diagonal element of the inverse Hessian matrix be used as a measure of each parameter’s sensitivity to noise. It notes that these elements do represent the geometric illconditioning that is inherent in geolocation based on light when near the equator and/or the equinoxes, including the effects of extracting cloudiness
ARTICLE IN PRESS P. Ekstrom / Deep-Sea Research II 54 (2007) 392–403
parameters from the same data used for estimating location. Thus they can be used to extrapolate accuracy estimates calibrated at one location and one season to other locations and seasons. However, the noise in field irradiance data available for geolocation is dominated by strongly correlated processes such as cloud movement and imperfectly compensated attenuation changes due to diving, so that the rest of the conventional error treatment, which assumes statistically independent data, is inapplicable. This work pursues an extension to the case of correlated data. 2.2. Data-taking Blue-light irradiance data were taken at a dark rural site on Shaw Island, WA, USA, 48.571N, 122.941W, using a now-obsolete archival tag, Lotek model LTD750. A successor tag model having the same kind of fluorescent-conversion light detector and temperature-compensated clock is the Lotek LTD2310. The spectral acceptance band of the detector peaks at 470 nm with a sharp red-side edge near 500 nm, imposing the necessary restriction to blue light. Irradiance data were recorded as integers on a logarithmic scale having 15 points per decade. Dynamic range of the measurement included full noonday sun and extended downward far enough (nominally seven to nine decades depending on temperature) to always easily include the useful part of the twilight transient. That is, irradiance values were always available for all sun elevation angles well past 51 even in the darkest of cloudy weather. Beginning on 15 February 2000, data for 458 days were recorded in approximately 40-day segments, sometimes separated by short gaps, spanning 475 total days on the calendar. The recording interval was 128 s, yielding 675 points per day. 2.3. Geolocation processing After recovery from the tag, data were converted to an irradiance scale having 32 points per decade for compatibility with later model tags. However, quantisation noise in the data is that associated with the initial coarser scale of 15 points per decade. Data were analysed using the template-fit method as outlined in Section 2.1. After a manual determination of the start and end of each approximately 40-day group, processing proceeded automatically without further human intervention. Days were not otherwise selected, so the results can be expected to
395
represent fairly the weather and seasonal geometric difficulties of about 1.3 years of the natural environment on land, including three equinoxes. The method was not ‘‘tuned’’ in any way to the particular location, which can be regarded as a typical example of the natural environment on land far from city lights. Initial parameter values needed to help start the Levenberg–Marquardt nonlinear least-squares fit were obtained from an initial threshold-method analysis of the day’s data. A threshold corresponding to a sun elevation angle of zero, rather than a negative-angle threshold chosen for optimum accuracy was used in order to guarantee a unique threshold-method solution in all seasons. Two template fits were performed for each day, leading to what are referred to below as the ‘‘northern-start’’ and ‘‘southern-start’’ solutions. Both used an initial value for longitude equal to the threshold-method longitude estimate, but two different initial latitude values were derived from extreme values of the threshold latitude estimate. These extreme values were computed from day lengths that were, respectively, artificially lengthened and artificially shortened by 14 min of time (3.51 of hour angle), an estimated maximum error due to weather noise. When, near the equinoxes, the resulting uncertainty range in latitude included both hemispheres, one fit was started in each hemisphere half way between the equator and the extreme value estimate. When both fell in the same hemisphere, fits were started at the 14 and 34 points of the range spanned. The results for each day depend only on data gathered that day, with no carryover from previous solutions. Nor was any information about the known location available to the method, which was (except for the initial piecing together of the data record to deal with the obsolete tag’s memory limits) implemented so as to be capable of functioning autonomously at any geographic location. When the day’s fits are entirely successful, the results of the northern-start and southern-start are expected to coincide. It is a property of any geolocation method based entirely on blue-light measurements taken below +31 solar elevation that its results are largely independent of horizon detail. It is a further property of the template-fit method that its results are also largely independent of cloud cover so long as the degree of sky obscuration does not change significantly during the sunrise or sunset transient.
ARTICLE IN PRESS P. Ekstrom / Deep-Sea Research II 54 (2007) 392–403
396
3. Fit results This section presents a first graphically oriented look at the geolocation estimates and the accompanying diagnostic information, together with a qualitative discussion of several effects and relationships that are directly visible in the figures. 3.1. Angular quantities
a 100
50
b
3.2. Inverse Hessian-matrix elements Fig. 3 shows the behaviour of the inverse Hessianmatrix diagonal elements corresponding to the latitude and longitude parameters. The square root of each element is plotted, since that is directly commensurate with an error size and is what will be used below in the estimate of standard deviation. The element for longitude, trace a, shows little variation with season, while that for latitude, trace b, becomes very large just on the winter side of the equinoxes, the season where most of the difficulties occur that are visible in Fig. 2. 3.3. Number of points in the fit, N Fig. 4 shows the number of points in the fit, again northern-start shown solid, southern-start shown dotted, and the number theoretically expected for the known latitude and season shown as a smooth solid line. The data points included are those that lie between 51 and +31 sun elevation angle as judged by the method’s template, which can vary between iterations as the method’s current estimated value of
c 0
-50 100
200
300 JD2000
400
500
Fig. 2. a: Longitude and b: latitude results. In each case the northern-start results are shown by a solid line, the southern-start by a dotted line, and the true value by a solid horizontal line. c: Solar declination. Ordinate is ‘‘Julian Day 2000’’, time in days from noon, January 1, 2000.
Covariance element root
Latitude, Longitude, or Declination, Deg
Fig. 2 shows a composite view of several angular quantities with the intent of exhibiting general trends and time correlations rather than fine detail. The three traces collectively marked a show the longitude results with the northern-start trace drawn as a solid line and the southern-start drawn dotted, as will be the convention in most graphs that show both northern- and southern-start results. The known true location is a solid horizontal line. These three traces are barely distinguishable on this scale, and longitude accuracy does not depend visibly on season. The traces marked b are the comparable display for latitude. Trace c is the solar declination; zeros of this trace are the equinoxes. Note that the accuracy in latitude degrades severely near (especially just on the winter side of) the equinoxes, and a spurious southern hemisphere solution appears, which the southern-start solution follows as it migrates south during the first equinox shown, then north during the second. Since both it and the northern-start result represent actual solutions in the nearly-symmetric geometrical situation near the equinox, we cannot expect our diagnostic informa-
tion to help us choose between them, and will have to look for other criteria such as track continuity to exclude the inappropriate one. In two cases even the northern-start solution converges to the wrong hemisphere, a result of geometric ill-conditioning near the equinox combined with weather noise. We can expect our error measures to warn us about times when this kind of behaviour is likely, but not to specify how to resolve it. In field results, the situation would have to be resolved using such considerations as track continuity. Threshold methods exhibit a similar spurious solution unless the irradiance threshold is chosen to correspond to zero solar elevation, which choice is disadvantageous for accuracy.
0.4 b 0.2 a 0 100
200
300 JD2000
400
500
Fig. 3. a: Longitude and b: latitude results for the square root of the diagonal element of the inverse Hessian matrix.
ARTICLE IN PRESS P. Ekstrom / Deep-Sea Research II 54 (2007) 392–403
Root of Fit Residual per Point
80
N
60 40 20
15 10 5 0 100
100
200
300 JD2000
400
500
Fig. 4. Number of data points in the fit. Solid line, northernstart; dotted line, southern-start; smooth solid line, theoretically expected number for the known latitude.
the latitude parameter converges toward its final value. The value shown is for the last iteration. Since the sun sets most rapidly near the equinoxes, there are fewer such points then, and more near the solstices when the sun sets more slowly. Note that near the equinoxes, the relatively large errors in the tag’s estimation of latitude result in errors in the assumed rate of sunset and consequently cause more or fewer than the correct number of points to be included. This can have an indirect further effect on accuracy by including data from outside the 51 to +31 ‘‘sweet spot’’ that occurs during twilight. 3.4. Residual sum of squares The residual sum of squares (Fig. 5) exhibits no pronounced variation with season, but rather varies widely from day to day, reflecting the effects of weather. Since the number of points in the fit varies with the season, and since for a given degree of noisiness in the data the size of the residual is expected to be proportional to the number of points, it is the residual per point corrected for the number of parameters estimated, R/(Np), that is displayed here. The graph shows the square root of this quantity, since that is commensurate with a standard deviation. We can see part of what this quantity is trying to tell us by examining a superposition of template and data for two successive days that have very different residuals. Fig. 1 (upper panel) shows Julian day 128 (JD2000, referenced to the millennium boundary), which has indistinguishable template traces (heavy, solid, and dotted) that fit the data well enough to lie atop each other and entirely obscure the data (lighter line) in the fit region between 51 and +31 sun elevation. The root residual per point is the same for both fits and is only 0.91 light units, very little larger than the lower limit of 0.62 set by the
397
200
300 JD2000
400
500
Fig. 5. Square root of the residual sum of squares per point. Solid trace, northern-start; dotted trace, southern-start. The solid horizontal line is the theoretically expected minimum residual of 0.62 LSB/point due to the original coarse 15/decade resolution of the recorded data.
resolution to which the irradiance was originally recorded. Fig. 1 (lower panel) shows the corresponding record for the following day. Here the data trace is significantly distorted from the smooth pattern typical of a clear day, presumably by clouds that erratically changed their effective fractional sky cover and/or opacity during the day. They also did so during the sunrise transient, which caused the difficulty we see. Beginning with the two different starts the two attempts converged to different local minima in the residual, neither representing a very satisfactory fit to the data record. During sunrise that record does in fact not look much like anything within the range of shapes reachable by stretching, translating, or raising the template in response to parameter changes. The root residual per point in one case was 12.4 units, in the other 11.4 units on a logarithmic scale where 32 units equal a factor of 10 irradiance change. 4. Error measures This section explores in more theoretical detail the properties of the diagnostic information available from the daily fits, proposes the daily standard deviation estimate, and explores its performance relative to the actual observed errors in position estimation. 4.1. The inverse Hessian-matrix diagonal element The error analysis usually associated with leastsquares fitting (Press et al., 1988) assumes that noise in the data points being fit is equivariant and statistically independent. When this hypothesis is satisfied, then one may obtain an estimate of the variance in a parameter estimate as the product of the diagonal element of the inverse Hessian matrix
ARTICLE IN PRESS 398
P. Ekstrom / Deep-Sea Research II 54 (2007) 392–403
(called in that reference the covariance matrix) associated with the parameter (either CLL or C ff ) times the common variance of all the (equivariant) data points. The analysis can be readily extended to the case where the data points have individual known variances, but an extension to the highly correlated errors in real field irradiance data is not available. Still, as there is reason to expect theoretically and also empirically as we have seen illustrated in the figures above, that the inverse Hessian matrix does capture the inherent error-sensitivity of the estimate that is being made, and in particular the dependence of error sensitivity of the latitude estimate on season. Theoretically it also should capture the effects of latitude itself on the error sensitivity of both estimates.
4.2. Comparison with threshold method error sensitivity In contrast to the template-fit method, a threshold method first estimates day length as the time between the sunrise and sunset crossings of an irradiance threshold. Then it chooses as its latitude estimate the value that gives the observed day length at the known season between the occurrences of the sun elevation angle that is taken to correspond to the adopted threshold. This involves solving for the latitude f the same transformation equation we saw in Section 2.1: sin a ¼ sin d sin f+cos d cos h cos f, where a is the assumed sun elevation angle, h the hour angle of the sun at sunset, and d the sun declination. In this context the quantity 2h is the day length expressed as an angle (11 of angle ¼ 4 min of time). The sun declination can be approximated as sin d ¼ sin(280.461+3601n/365. 25) sin(23.441), where n is JD2000, the Julian Day referred to the start of year 2000. Unless the sunrise–sunset sun elevation (a) is exactly zero (geometric sunset), no closed-form solution is available and this equation must be solved iteratively. However, one can obtain a closed form for the derivative of implied latitude with respect to day length that is valid for any value of a: df/dh ¼ cos d cos2 f sin h/(sin dsin a sin f). When this derivative becomes large, a simple linear error model, Df ¼ df=dh Dh; predicts large errors in latitude caused by small errors in day length. In particular, the derivative is infinite when sin d ¼ sin a sin f.
This measure reflects the same geometric situation that degrades the accuracy of the template-fit method: near the equinox, twilight details do not depend strongly on latitude, so little latitude information can be extracted from them. Accordingly it is not surprising that this measure behaves qualitatively much like the inverse Hessian root. The match between them is closest for a threshold set near 5.51. Fig. 6 plots the two functions with the threshold method sensitivity multiplied by an arbitrary constant that makes the two agree near the summer solstice. The light solid trace is the corresponding threshold method sensitivity for a threshold of +31, multiplied by the same arbitrary factor. The graph illustrates that the error sensitivities to angular errors depend on the threshold chosen. These sensitivities exhibit their largest values at different sun declinations, and therefore at slightly different times of the year, a fact first pointed out by Hill and Braun (2001). This suggests the possibility of using the threshold method with a high threshold to fill in during the seasons when the template-fit method is at its worst. Naively assuming that a threshold method is as good as a template fit during the summer (which it appears not to be; see rough indications in Ekstrom, 2004) and that the day length accuracy attainable with a +31 threshold is as good as with one at 5.51, all of which is suggested visually by traces a through c of Fig. 6, it would appear advantageous to use the threshold method with a 31 threshold not only during the equinoxes, but for the entire interval from Julian day 270 through 440. Unfortunately the rate of change of irradiance with time is about a factor of 5.4 less rapid at +31 sun elevation angle as it is at 5.51, so comparable errors in irradiance can be expected to lead to errors in day length that are larger by that factor of 5.4 at +31 than they are at 5.51. To illustrate the impact of this effect, the light dashed trace of Fig. 6 is a copy of the +31 angular sensitivity trace, but multiplied by 5.4 for fair comparison with the 5.51 trace assuming comparable errors in determining irradiance. Even assuming that the absolute errors in a 5.51 threshold method are no worse than those of the template-fit method (which is visually suggested by the choice of arbitrary multiplier, but not asserted as true), when the difference in irradiance slope is taken into account the expected latitude error of a +31 threshold method is always substantially greater than either alternate method except very near the equinox, where none of
ARTICLE IN PRESS
Latitude SensitivityMeasure
P. Ekstrom / Deep-Sea Research II 54 (2007) 392–403
399
0.4
0.3
d d
0.2 c 0.1
a,b
c
0 50
100
150
200
250
300 JD2000
350
400
450
500
1 Fig. 6. The inverse Hessian-matrix root for latitude (a, heavy solid trace, identical to Fig. 2 trace b), compared with 20 times the latitude sensitivity measure for 5.51 threshold method (b, heavy dashed) and +31 threshold method (c, light solid trace) computed for the known true latitude. Trace d is a copy of c multiplied by a factor of 5.4 to account for the different irradiance slope at +31.
the three are useful, and therefore offers no advantage. 4.3. Proposed estimates of template-fit geolocation error In order to make use of the inverse Hessian elements in estimating geolocation error, it is necessary to have a measure of data noisiness to multiply by. The situation does not permit the usual first choice, the independently determined variance of uncorrelated, equivariant data, because the available data is at minimum strongly correlated. The most obvious (and traditional) candidate for replacing that first choice is some multiple of the residual sum of squares. As we have seen, it does measure the deviation of the data record from the form predicted by the template, and thus is expected to be large on a ‘‘bad day’’ when the assumptions underlying that template, notably that absorption due to clouds or water absorption does not change during twilight, are not realised. Making that choicep leads first to the un-calibrated p error measures, E L ¼ ðC LL RÞ and E f ¼ ðC ff RÞ. It would be more conventional to divide the residual sum of squares by the number of points in the fit, less the number of parameters estimated. That would be an appropriate choice if the dominant noise source were uncorrelated. In the present case, the noise is ordinarily over-sampled in time and strongly correlated in the record. These uncalibrated estimates will be eventually multiplied by two empirically determined constants, kL and kf, respectively, to obtain
estimates of the standard deviation to associate with the day’s position estimate. That is, the eventual calibrated estimates will be p sL ¼ k L E L ðC LL RÞ and and sf ¼ k E , leading to s ¼ k f f L L p sf ¼ kf ðC ff RÞ. 4.4. Avoiding spurious solutions As noted above, the erroneous southern-hemisphere solutions are real solutions that arise from the geometric near symmetry present near the equinox when the sun lies above the equator. Thus the method provides no diagnostic information that can serve to automatically choose the correct one of two solutions that disagree. However, while the method may be in doubt as to which hemisphere the animal is in, the experimenter is rarely in doubt. Track continuity can ordinarily resolve the question in situations where the issue arises at all. In line with this, the development from here on will make use of our comparable knowledge that the test site is in the northern hemisphere, and will do something that would be available to an experimenter in a real field situation. It will concentrate on the northern-start solutions in the rest of this paper, and will also reject as invalid the two northern-start solutions that ended up in the southern hemisphere during the time of equinoctial ill-conditioning. 4.5. Bias in the position results Given an estimate for the individual standard deviations of each day’s position estimates, the
ARTICLE IN PRESS P. Ekstrom / Deep-Sea Research II 54 (2007) 392–403
optimum estimate for the mean is an average weighted by the inverses of those standard deviations. Thus the observed systematic biases in the position estimates can be obtained in terms of the deviations D from the known position as BL ¼ SðDL =sL Þ=Sð1=sL Þ and Bf ¼ SðDf =sf Þ=Sð1=sf Þ. These sums can be computed now, while the calibration factors are still unknown because those factors cancel in the final division. That is, only the un-calibrated product error measure is needed. The results for the northern-start data are BL ¼ 0:041 and Bf ¼ 0:151. The standard deviations in those p weighted means are sL ¼ ðS1Þ=Sð1=sL Þ and p sf ¼ ðS1Þ=Sð1=sf Þ, where (S1) is used to represent the number of points in each of the sums. However, these quantities cannot be evaluated without values for the calibration factors, which are computed next.
400
Rank order
400
200
0 -40
-20 0 20 Normalized longitude error
40
Fig. 7. Two superimposed and barely distinguishable plots regarding longitude results: rank order plot of southern-start values of ðDL BL Þ=E L (solid), and best-fit scaled Gaussian CDF (dashed), which has s ¼ kL ¼ 0:90.
4.6. Determining kL and kf
Rank order
400
Once the proposed estimates of data variability are correctly calibrated, the normalised, bias-corrected quantity ðDL BL Þ=sL will be a random variable with zero mean and unit standard deviation. Multiplying through by the constant kL one may conclude that the quantity ðDL BL Þ=E L will have mean zero and standard deviation kL. Similarly, the quantity ðDf Bf Þ=E f will have mean zero and standard deviation kf. Thus one can determine the values of kL and kf by examining those distributions. The solid line in Fig. 7 is a rank order plot, that is, an empirical un-normalised Cumulative Distribution Function (CDF) of the first of these quantities. The superimposed dashed line is the theoretical CDF of a normal distribution (multiplied by the total number of data points to match vertical scaling) that is a best fit to it. While detailed examination of the tails shows some divergence between them, the two distributions can be barely distinguished in the figure. The best-fit value of standard deviation is 0.90, which one may identify with kL. Fig. 8 is a similar plot of the normalised latitude errors, with two outliers that lie in the southern hemisphere removed from the set, as explained in Section 4.4. Again there is a superimposed Gaussian, best-fit with a standard deviation of kf ¼ 1:17. Again the cumulative distributions largely coincide, but in this case there is a much longer negative tail on the data CDF involving about 3 of the remaining 456 points.
200
0 -40
-20 0 20 Normalized latitude error
40
Fig. 8. Rank order plot of southern-start values of ðDf Bf Þ=E f (solid), and best-fit scaled Gaussian CDF (dashed), which has s ¼ kf ¼ 1:17.
Inserting the calibration values just determined, p the two noise p estimates become sL ¼ 0:90 ðC LL RÞ and sf ¼ 1:17 ðC ff RÞ. As will be argued shortly, the following information should accompany any such calibration: the data used were taken at a recording interval of 128 s on land at mid-latitudes (Shaw Island as described) in approximately 40-day continuous blocks providing usable records for 458 of 475 total days. The missing days were determined by the accidents of data collection and did not represent any deliberate selection. Now with the scaling constants in hand, it is possible to evaluate the expected standard deviation in the weighted mean used to calculate the bias values, and to state those values with their one-s
ARTICLE IN PRESS P. Ekstrom / Deep-Sea Research II 54 (2007) 392–403
401
Standard deviation estimate
100
10
1
0.1 50
100
150
200
250
300 JD2000
350
400
450
500
Fig. 9. The daily standard deviation estimates for longitude (squares, solid line is its 0.421 median) and latitude (crosses, dashed line is its 0.951 median). Two latitude standard deviation estimates of 1221 and 2421 occurring near Julian day 430 lie off the graph.
limits as BL ¼ 0:041 0:021 and Bf ¼ 0:151 0:041. The template-fit method appears empirically at least to be reasonably unbiased. The values the daily position error estimates take on for this data are plotted in Fig. 9. The error estimates for latitude vary widely from values similar to or even slightly better than those for longitude near the summer solstice to very large values near the equinoxes. For this data the median value of sL is 0.421, the median value of sf is 0.951. 5. Discussion The development above yields an estimate for each day’s uncertainty in template-fit geolocation results, based on data obtained at a particular midlatitude dark site on land in all seasons and all weather. The method used is applicable to the completely automated extraction of geographic position from field data on the time sequence of blue-light irradiance. It remains to consider how broadly applicable these error measures may be to locations and conditions other than those explicitly represented by the data examined. The method is based on a geophysical model of sunset that relates sun elevation angle to observable irradiance. For sufficiently strong scattering (blue light) the shape of the resulting template function depends only on the ratio of two independently measurable quantities, the stratospheric scale height and the radius of the Earth, and it has no available adjustments. There is no need to ‘‘tune’’ the template to any particular location, and no obvious way to do so. The template shape is very robustly determined so long as the scattering is strong
enough. For example, a significant increase in stratospheric aerosol content, perhaps due to a volcanic eruption, will move the scattering peak to higher altitude and change the value of the ‘‘cloudiness’’ parameter, but will not alter the shape of the template function. The next step in developing the template uses an exact astronomical coordinate transformation to rewrite the sun elevation in terms of latitude, longitude, season, and time of day. As a result the place where the data happened to be measured can be considered as an entirely typical mid-latitude site. Since the effect of a change in longitude on the data record is entirely equivalent to resetting the clock, extension to any other longitude is immediate. Both latitude and season affect the amount of information that there is in the data record regarding latitude or, from another viewpoint, the degree of geometric ill-conditioning that occurs. This has a major effect on the accuracy with which the template parameters can be determined from the data, but that error sensitivity is reflected in the inverse Hessian-matrix elements which form part of the proposed error estimators. The fact that these quantities successfully represent the variation in accuracy with season has been examined and verified in this work, and the variation with latitude is just another consequence of the same coordinate transformation. It is reasonable to expect, though this work has not separately verified it, that the error measures studied here will adequately account for the geometric effects of site latitude as well. The most heuristic aspect of this work is the use of the residual sum of squares as a measure of data noise, and that is also the aspect whose extrapolation
ARTICLE IN PRESS 402
P. Ekstrom / Deep-Sea Research II 54 (2007) 392–403
to other situations is the least secure. In an important sense this quantity cannot directly see the noise we care about. If a pattern of noise in the data is such that it can be approximated successfully by a parameter change, that parameter change will occur when the template is fit to the data, causing an error in the parameter value. That noise pattern will be absorbed, and will not appear in the residual. If the change is to the latitude or longitude parameter, it will cause a geolocation error. What is left in the residual is the effect of those noise patterns that cannot be absorbed by parameter changes, and that therefore have had no effect on the geolocation result. When using the residual sum of squares as a measure of noise for estimating parameter errors, one necessarily assumes that total noisiness in the data is on the average partitioned in some approximately fixed ratio between the harmless data patterns that do affect its value and the harmful ones that do not. That assumption is implied by using one as a measure of the other. If the character of the noise changes, the ratio of harmful to harmless noise components in each day’s data may change, and may affect the appropriate value of the calibration constants, kL and kf. As a result, weather with very different statistical behaviour, or imperfectly compensated diving noise as opposed to weather noise, may call for a different value of calibration constants even though both are examples of highly correlated noise. Therefore one should regard these present k values with greater suspicion, the less similar the natural noise processes in a new situation are to weather at the test site where the data used here were taken. We can already see a hint of this effect in the different sizes of the two k constants; the noise encountered here had a greater effect on the latitude estimate than on the longitude. Since all other effects would seem to be accounted for in other ways, the only obvious reason for the difference is an accident that the particular data patterns most harmful to latitude determination were more prevalent in the weather noise encountered than were patterns harmful to the determination in longitude. Therefore the calibration parameters should be re-determined at every opportunity where there is reason to suspect a difference in noise properties. They should also be re-determined for any new data interval, and any determination of k should always report the data interval along with anything that might distinguish the noise processes involved. Doing so in a variety of situations will explore the
range of variation of the calibration parameters, along with the nature of the situations in which each value is appropriate. Until this exploration has provided a broad understanding, the values reported here can be used with caution to begin providing a way of distinguishing good days from bad days, and of forming some idea of how good and how bad they are, in field experiments on animals. 6. Conclusion Error estimates for each daily position are available as part of the template-fit method of light-based geolocation. These estimates are sL ¼ p p 0:90 ðC LL RÞ and sf ¼ 1:17 ðC ff RÞ, where the symbols are as defined in Section 1.1. The numerical proportionality constants given have been determined for the particular case of observations made at a mid-latitude site (Shaw Island, WA, USA) on land every 128 s in all seasons. For those data, the estimates had median values of 0.421 in longitude, 0.951 in latitude. It is desirable to re-determine the calibration constants in any situation where the noise correlations are suspected to differ significantly from those of the data used here. Acknowledgements Thanks are due to Raymond Glaze for accumulating the long-term irradiance records used in this study, to Hannah Halliday for helping scan the data for interesting features, to Steven Teo for suggesting an extension of the original poster version, and to a very helpful anonymous reviewer.
References Ekstrom, P.A., 2002. Blue twilight in a simple atmosphere. In: Shaw, J.A. (Ed.), Proceedings of SPIE Annual Meeting 2002: Atmospheric Radiation Measurements and Applications in Climate, vol. 4815, Paper 14. SPIE, Bellingham, WA, USA, pp. 73–81. Available at /http://www.lotek.com/library.htmS [doi:10.1117/12.482305]. Ekstrom, P.A., 2004. An advance in geolocation by light. In: Yasuhiko, N. (Ed.), Proceedings of the International Symposium on Bio-logging Science. Memoirs of the National Institute of Polar Research, Special Issue, 58. National Institute of Polar Research, Tokyo, Japan, pp. 210–226. Available at /http://polaris.isc.nipr.ac.jp/ penguin/oogataHP/pdfarticles/22p210-226.pdfS. Hill, R.D., Braun, M.J., 2001. Geolocation by light levels—the next step: latitude. In: Sibert, J.R., Nielsen, J.L. (Eds.),
ARTICLE IN PRESS P. Ekstrom / Deep-Sea Research II 54 (2007) 392–403 Electronic Tagging and Tracking in Marine Fisheries. Kluwer, Dordrecht, The Netherlands, pp. 315–330. Musyl, M.K., Brill, R.W., Curran, D.S., Gunn, J.S., Hartog, J.R., Hill, R.D., Welch, D.W., Eveson, J.P., Boggs, C.H., Brainard, R.E., 2001. Ability of archival tags to provide estimates of geographical position based on light intensity. In:
403
Sibert, J.R., Nielsen, J.L. (Eds.), Electronic Tagging and Tracking in Marine Fisheries. Kluwer, Dordrecht, The Netherlands, pp. 343–368. Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T., 1988. Numerical Recipes in C. Cambridge University Press, Cambridge, UK.