- Email: [email protected]
Efficient radiative transfer model inversion for remote sensing applications
Remote Sensing of Environment 113 (2009) 2527–2532
Contents lists available at ScienceDirect
Remote Sensing of Environment j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / r s e
Efficient radiative transfer model inversion for remote sensing applications John Hedley a,⁎, Chris Roelfsema b, Stuart R. Phinn b a b
School of Biosciences, University of Exeter, EX4 4PS, UK Center for Remote Sensing and Spatial Information Science, School of Geography, Planning and Environmental Management, University of Queensland, Brisbane, Australia
a r t i c l e
i n f o
Article history: Received 13 May 2009 Received in revised form 16 July 2009 Accepted 18 July 2009 Keywords: Radiative transfer Inversion Look-up table (LUT) Successive approximation
a b s t r a c t A simple method for efficient inversion of arbitrary radiative transfer models for image analysis is presented. The method operates by representing the shape of the function that maps model parameters to spectral reflectance by an adaptive look-up tree (ALUT) that evenly distributes the discretization error of tabulated reflectances in spectral space. A post-processing step organizes the data into a binary space partitioning tree that facilitates an efficient inversion search algorithm. In an example shallow water remote sensing application, the method performs faster than an implementation of previously published methodology and has the same accuracy in bathymetric retrievals. The method has no user configuration parameters requiring expert knowledge and minimizes the number of forward model runs required, making it highly suitable for routine operational implementation of image analysis methods. For the research community, straightforward and robust inversion allows research to focus on improving the radiative transfer models themselves without the added complication of devising an inversion strategy. © 2009 Elsevier Inc. All rights reserved.
1. Introduction Inversion of physics-based radiative transfer models is an area of rapid development in remote sensing of both aquatic and terrestrial environments (Liang, 2007). In these approaches a parameterized forward model for reflectance such as Hydrolight (for aquatic applications) or PROSPECT, PROSAIL or DART (for terrestrial canopies) takes a series of parameters describing the optical properties of participating media or canopy structure and the model defines a mapping from parameter space to spectral space (Fig. 1a) (Darvishzadeh et al., 2008; Gastellu-Etchegorry et al., 2003; Mobley et al., 2005; Zhang et al., 2008). Image analysis then uses a search algorithm to find the parameter space location which minimizes the distance of the corresponding spectral space location from the image pixel reflectance, measured by a distance metric such as Euclidean distance. Two distinct approaches to the search algorithm can be identified: 1) pre-calculation of reflectance look-up tables (LUTs) by repeated runs of the forward model with systematically differing parameter values (Mobley et al., 2005), optionally with interpolation across tabulated points (Gastellu-Etchegorry et al., 2003; Liang et al., 2006) and 2) spectral matching by successive approximation using optimization techniques such as the Downhill Simplex or Levenberg–Marquardt (L–M) algorithm (Brando et al., 2009; Goodman & Ustin, 2007; Klonowski et al., 2007; Lee et al., 1999; Wolfe, 1978), which is feasible only if the forward model is sufficiently fast to permit many runs
⁎ Corresponding author. E-mail address: [email protected] (J. Hedley). 0034-4257/$ – see front matter © 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.rse.2009.07.008
for each image pixel. With either approach, for a complex radiative transfer model, both the forward modeling and the inversion may take substantial computational effort. Routine or large-scale operational image analysis by physics-based methods therefore demands the development of efficient approaches for both modeling and inversion. Often the forward model is simplified or approximated to facilitate inversion. In contrast, we develop a generic inversion technique which allows the use of any forward model and hence model accuracy need not be compromised. A key problem in the efficient representation of the parameter to reflectance mapping function is that it is not generally clear a priori how to subdivide the parameter space in order to evenly populate spectral space from the forward model. For example, modeled shallow water reflectance typically has an inverse exponential relationship with depth (Mobley, 1994), so a regular subdivision of a depth parameter would over-sample the deep water spectral space region with similar reflectance spectra produced from unnecessary forward model runs. At the same time the shallow water spectral space region would be relatively under-sampled with a higher discretization error in the tabulated reflectances (Fig. 1b). The interaction between multiple parameters (e.g. depth vs. water clarity) makes the general problem of efficient and accurate LUT construction very hard to tackle by analytical means. To address this problem we present an automated construction procedure for adaptive look-up trees (ALUTs) that evenly distributes and minimizes the discretization error in spectral space and facilitates fast inversion (Fig. 1). The construction technique can be applied to any forward radiative transfer model that is parameterized by a series of real and integer values. The following sections describe 1) the requirements of a forward model to which the ALUT construction
2528
J. Hedley et al. / Remote Sensing of Environment 113 (2009) 2527–2532
a set of set of m real-valued parameters, α1, α2…αm, and any number of discrete integer parameters, i1, i2…, etc. As a first step, any discrete valued parameters are combined into a single integer parameter by i = i1 + i2N1 + i3N1N2 + … so that i = 1, 2…N, where N1, N2… are the maximum value each integer parameter can take and N is the product of the maximum values. This gives a set N of mapping functions from the real-valued parameter space to remote sensing reflectance in spectral space, fi : Rm YRn hα 1 ; α 2 N α m iiRrs
Fig. 1. Adaptive LUT construction for a forward model of one continuous parameter, α, mapping to two spectral bands, hb1 ; b2 i. (a) The mapping function from parameter space to spectral space which the LUT seeks to accurately represent. (b) Discretization by regular steps in parameter space over-samples some regions of spectral space and under-samples others. (c) and (d) iterative adaptive construction fills in the gaps to approximate an even distribution of points in spectral space and efficiently represent the function.
ð1Þ
where fi is the forward model for particular discrete value of i and each f i consists of an assumed continuous mapping from mdimensional parameter space to n-dimensional spectral space, with Rrs being the spectral remote sensing reflectance in n-vector form. The separation of discrete and continuous parameters is necessary because the adaptive subdivision can only be applied to the continuous parameters. In essence, a series of separate ALUTs are constructed, one for each value of i, this accommodates forward models that contain discrete parameters, for example where a model component is based on a choice of one from a set reflectance spectra. The only additional requirement is that a priori limits for the continuous parameters are defined, so that α MIN V α i V α MAX , this is the i i same as is required for any LUT methodology and is also often used in successive approximation techniques to constrain the mapping function domain. 2.2. Adaptive look-up-table construction
method can be applied 2) the adaptive construction algorithm, and 3) the methodology for building and searching the tree representation of the ALUT. The subsequent sections present an illustrative example of ALUT construction using a shallow water radiative transfer model (Lee et al., 1999) and its application to hyperspectral airborne imagery of a coral reef. 2. Methods 2.1. Mapping from parameter space to spectral space The adaptive LUT construction can be applied to any forward model for n-band remote sensing reflectance modeled dependent on
ALUT construction begins with an initial regular subdivision of parameter space (as shown for 1-dimensional parameter space in Fig. 1b), then an iterative algorithm repeatedly assesses which regions in parameter space correspond to currently under-sampled regions of spectral space, the most effective parameter is subdivided in that region, the new points are added to the ALUT and the iteration continues (Fig. 1c,d). The key to implementing the subdivision algorithm for mdimensional parameter space is that the ALUT structure is represented as a hierarchical voxelation of the parameter space, where voxels are m-dimensional parameter bounding boxes, the 2m vertices of which are current points in the LUT (Fig. 2). To assess the local point distribution in spectral space, change in Rrs across the voxel is assessed by taking the mean distance in spectral space when changing
Fig. 2. Hierarchical voxelation of parameter space for a forward model of three continuous parameters hα 1 ; α 2 ; α 3 i. Black dots show the parameter space location of points in the ALUT where Rrs has been calculated, open circles show where subdivision is occurring and points need to be calculated for the next iteration.
J. Hedley et al. / Remote Sensing of Environment 113 (2009) 2527–2532
a parameter from its lower bound on one side of a voxel to its upper bound on the other side. This occurs along 2m − 1 edges of the voxel, one edge for each high-low combination of the parameters other than the one being assessed. Using Euclidean distance as the distance metric, the mean change in Rrs for parameter i is
DðiÞ =
1 m−1
2
m−1 2X
jR rs ði; j; 0Þ − Rrs ði; j; 1Þj
ð2Þ
j=1
with Rrs (i, j, 0|1) being the remote sensing reflectance corresponding to the vertices at the low and high parameter value end of the jth voxel edge which describes a change in parameter i. The mean change in Rrs for each parameter, D(i), is calculated for each parameter in each voxel when the voxel is created. The values are stored in a list which is continually updated and sorted during the construction process. Voxels at the top of the list represent the parameter space region which currently has the greatest discretization error in spectral space. The LUT is grown by taking the voxel at the top of the list and bisecting it into two new voxels at the midpoint of the parameter value that represents the greatest mean spectral change across the voxel (Fig. 2). When a voxel is subdivided 2m − 1 new vertices are required, their parameter locations are the midpoints of the subdivided edges, and 2m − 1 forward model calculations are required to generate the corresponding Rrs spectra. The mean spectral space distances for each parameter, D(i), are calculated for the new voxels (Eq. (2)) and added to the sorted list, and the construction procedure continues. To initialize the ALUT construction in m-dimensional space a regular grid of voxels is created with corners at the extremes of the parameter limits, typically this grid can just be a single voxel (Fig. 2a). When the forward model contains a discrete parameter, i.e. N N 1, the ALUT construction algorithm is repeated independently N times, with the forward mapping function being fi, i = 1,2,…N. Essentially N separate ALUTs are constructed, these are merged into a single treestructured LUT as described in the next section. 2.3. Search optimization by binary space partitioning trees
2529
when all tabulated points below the current node have been checked or ruled out. During the search, the Euclidean distance to the current best found Rrs vector is maintained. At each node the algorithm will descend down a branch only if the image reflectance to match is on the same side of the node's dividing plane as the branch, or is closer to the dividing plane than the current best matched Euclidean distance. Only then can points on the far side of the dividing plane possibly be a better match than the current best match. Note that the choice of np, maximum sample size for the PCA, may have a small effect on inversion time but has no effect on accuracy. The suggested value of 30,000 is probably suitable for all applications. 3. Example application The following sections illustrate an example ALUT inversion of a 5-parameter shallow water radiative transfer model applied to a hyperspectral airborne image of a tropical coral reef (Fig. 3). Speed and accuracy of ALUT bathymetric retrievals is compared to two implementations of the Levenberg–Marquardt successive approximation algorithm, with tide-corrected sonar depth transects used for ground-truth. 3.1. Imagery and forward model The adaptive LUT construction and inversion algorithm was applied to a 11278 × 5236 pixel 17-band CASI image (430–710 nm, FWHM 10–20 nm) acquired in 2002 of Heron Reef, Australia (23.44° S, 151.91° E, Fig. 3) with the ALUT constructed to represent a variant of Lee's semianalytic forward model for shallow water remote sensing reflectance (Lee et al., 1998, 1999) where subsurface remote sensing reflectance is estimated from subsurface solar zenith angle (θw =32°) subsurface view angle (θ=0°), and depth in meters, H, by ( "
# )! C 1 Du ðλÞ 1 − exp − + κ ðλÞH cosθw cosθ ( " # ) 1 1 DBu ðλÞ + + ρðλÞexp − κ ðλÞH π cosθw cosθ
dp ðλÞ rrs ðλÞ≈rrs
ð3Þ
To provide a spectral matching search algorithm that is more efficient than an exhaustive search the storage of all the reflectance vectors is post-processed into the structure of a binary space partitioning (BSP) tree. This is achieved by the following algorithm: 1. A list is initialized containing all Rrs vectors in the ALUT linked to their corresponding generating parameter values. 2. A principal components analysis (Manly, 1994) is performed on a random selection of np of the Rrs vectors in the list (example: np = 30,000), or the entire list if the current number of vectors in the list is less than np. 3. A dividing plane in spectral space is established, such that the normal to the plane is the first principal component axis and the plane passes through the mean of the selected Rrs vectors. 4. The list of points is subdivided into two lists dependent on which the side of the dividing plane the Rrs points lie on. 5. The procedure is applied recursively to the two new lists, starting at Step 2, unless the number of points in the remaining lists is below a predetermined limit (example: 20 points) or the tree depth reaches maximum limit (example: 200 subdivisions). In the resulting tree every node has an associated dividing plane and subdivides into two branches, the actual tabulated reflectances are stored in groups of 20 or less in the terminal nodes (leaves). Searching the BSP tree to find the closest match to an image reflectance is efficient since the global optimum can be found without checking every tabulated Rrs vector. The search traverses the tree recursively downwards to the leaves, stepping back up one level
Fig. 3. CASI Image of Heron Reef and location of sonar transects. The full image (top panel) is 11278 × 5236 pixels of 1 m × 1 m in 17 bands 430–719 nm with variable band widths from 10 to 20 nm. The images here are monochrome versions of an RGB image based on bands centered at 665, 564 and 480 nm. The central dark region in (1) is the island itself.
2530
J. Hedley et al. / Remote Sensing of Environment 113 (2009) 2527–2532
where the subsurface remote sensing reflectance for optically deep C B water, rdp rs (λ), and Du(λ) and Du(λ), can be estimated from dp
rrs ðλÞ ≈ ½0:084 + 0:170uðλÞuðλÞ
ð4Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C Du ðλÞ ≈ 1:03 1 + 2:4uðλÞ
ð5Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi B Du ðλÞ ≈ 1:04 1 + 5:4uðλÞ
ð6Þ
tance using the relation R rs = 0:54R 0þ = kQ where Q = 4 (Mobley, 1994) and k is a spectrally independent factor to convert above-sur face to below-surface diffuse reflectance, k ≈ R 0þ = Rð0 − Þ ≈ 0:7, estimated from multiple runs of PlanarRad, an open-source planeparallel radiative transfer model for directional radiance similar to the commercial software Hydrolight (Hedley, 2008; Mobley, 1994). 3.2. ALUT and Levenberg–Marquardt inversion The structure of the forward model to be encoded by the ALUT in the form of Eq. (1) is therefore
where uðλÞ = bb ðλÞ = ½aðλÞ + bb ðλÞ;
κ ðλÞ = aðλÞ + bb ðλÞ
ð7Þ
where a(λ) and bb(λ) are the total absorption and backscattering coefficients and estimated as the sum of the effects of water constituents by, aðλÞ = aw ðλÞ + ½a0 ðλÞ + a1 ðλÞ lnP P + G exp ½−0:015ðλ − 440Þ ð8Þ Y
bb ðλÞ = bbw ðλÞ + X ð400= λÞ
ð9Þ
The pure water absorption and backscatter aw(λ) and bbw(λ) are considered known (Buiteveld et al., 1994; Pope & Fry, 1997), P, G, X and H are free parameters to be solved by inversion, these describe absorption due to phytoplankton, CDOM absorption, particulate backscatter and depth respectively. The particulate backscatter slope is fixed, Y=1, and the phytoplankton spectral absorption parameters, a0(λ) and a1(λ), are tabulated in (Lee et al., 1998). For the substrate cover reflectance in Eq. (3), ρ(λ), a linear mixture model is used(Wettle & Brando, 2006; Brando et al., 2009) where ρðλÞ = Ei1 ðλÞm + Ei2 ðλÞð1 − mÞ is a linear mix, in proportion m, of any two benthic reflectances, Ei1 (λ), Ei2 (λ), taken from a set of thirteen typical reef benthic reflectances including several of sand, live and dead coral, algae and seagrass, so i1 = 2…13 and i2 = 1… i1 − 1 giving 78 pair combinations, describable by a single integer parameter, i = 1…78. The forward model outputs (Eq. (3)) were translated to vectors of above-surface remote sensing reflectance, Rrs using Rrs ð jÞ = 0:5rrs λj = 1 − 1:5rrs λj where λj is the CASI band j wavelength center (Lee et al., 1999). The CASI image had been previously corrected to units of above-surface diffuse reflectance, R(0+) Joyce (2004) and so was converted to above-surface remote sensing reflec-
fi : R5 YR17 hP; G; X; H; miiRrs for i = 1 N 78
ð10Þ
The lower and upper limits of the parameters used to establish the single initial parameter voxel (Fig. 2a) were: 0 ≤ P ≤ 1.5, 0 ≤ G ≤ 0.08, 0 ≤ X ≤ 0.007, 0 ≤ H ≤ 20, 0 ≤ m ≤ 1. Hence the maximum retrievable depth is 20 m, the ranges of G and X were chosen for relatively clear waters while the upper bound of P, the component of a(440) due to phytoplankton, was set substantially higher than would normally be expected in these waters to illustrate features of the ALUT construction (Fig. 4). Two ALUTs of 5 × 106 spectral space points were generated from the above parameterization, this size of ALUT being the maximum the computer's memory could accommodate. The first, henceforth called “Automatic ALUT”, followed the previously described procedure exactly with automatic subdivision of parameters based on the distribution of points in spectral space. However, the strong effect of water absorption on reflectance above 700 nm results in an extreme treatment of depth parameter subdivision: shallow depths are generated with steps as fine as 5 mm, while other parts of the ALUT contain no depths between 10 and 20 m (Fig. 4a). This is an accurate picture of the effect of depth on the spectral shape and magnitude of reflectance, but raises the question as to whether accuracy could be improved by a more practical discretization of the depth parameter. To answer this, a second ALUT was constructed where depth subdivision was forced to a minimum of five subdivisions and constrained to eight, so throughout the “Detailed Depth ALUT” depth is at worst represented to an accuracy of 0.63 m all the way to 20 m, but no finer than 0.08 m in shallowest regions (Fig. 4b).
Fig. 4. Example slices through (a) the automatically generated ALUT and (b) ALUT generated with constraints on the depth parameter subdivision. The diagrams show two axes in parameter space, depth and the phytoplankton parameter, P. The values for other parameters are fixed at G = 0, X = 0, m = 0 and the bottom reflectance is pure sand. Boxes represent cross-sections of voxels (Fig. 2) and line intersections are the point locations in the ALUT where Rrs has been calculated.
J. Hedley et al. / Remote Sensing of Environment 113 (2009) 2527–2532
2531
Fig. 5. Depth estimates versus sonar derived bathymetry, n = 5057. Dashed lines are the regression slope plus and minus one standard error, Δy is the mean of the absolute value of the errors between estimated and sonar bathymetry.
To compare the ALUT inversion approach to alternate methods the model was also inverted by two publicly available C code implementations of the L–M algorithm, the Netlib Cephes package1 and Levmar2.
4. Results and discussion The Levmar L–M optimization and both ALUT inversions resulted in almost identical accuracies for bathymetric retrieval assessed over 5057 sonar data points (Fig. 5). In all three cases r2 = 0.91, the slight difference between the L–M and ALUT inversions with respect to regression slopes (1.078 vs. ~0.93) and in the mean of absolute depth error (0.90 vs. ~0.77) can be partly accounted for by the lack of an upper depth parameter constraint on the L–M inversion. Constraining ALUT inversion to depths less than 20 m is essentially the incorporation of a priori information, since it was known all sonar depths were less than 20 m. However, incorporating the same upper bounds on the L–M inversion actually produced worse results, due to algorithm convergence issues in the vicinity of the constraints. Since the ALUT BSP tree search algorithm returns the global optimum it is free from convergence issues. While the L–M successive approximation was able to produce slightly closer spectral matches to the image reflectances than the ALUT (Fig. 6) this did not translate into higher bathymetric accuracy. Both methods were able to obtain spectral matches close to or within image noise levels, estimated from the covariance matrix of pixels over deep homogeneous water (Fig. 6). Hence the ability to produce a spectral match was not limiting for either approach. The Automatic ALUT inversion displayed clear discretization of depth retrievals (Fig. 5b) due to the depth parameter subdivision structure of the ALUT (Fig. 4a). However while the Detailed Depth ALUT produced a cosmetically superior scatterplot (Fig. 5c), accuracy assessed by the three metrics of r2, regression slope and mean absolute error was virtually identical to the Automatic ALUT (Fig. 5b). The distribution of spectral match errors was also virtually identical (Fig. 6). This indicates that the automatically generated structure of the ALUT does give an efficient representation of the models capability to retrieve environmental parameters. The very similar spread of the L–M and ALUT bathymetric retrievals points to a deficiency in the underlying forward radiative transfer model or image noise as the ultimate accuracy limiting factor. 1 2
URL: http://www.netlib.org/cephes. URL: http://www.ics.forth.gr/˜lourakis/levmar.
Fig. 6. Spectral matching errors over pixels containing sonar data points, expressed as distribution of Euclidean distance between image and matched spectral reflectances, n = 4146. Environmental and image noise is estimated as the noise-equivalent deltareflectance, NEΔRrs, by the covariance of spectral reflectance over a homogeneous deep water region, and then by calculating the radius of the bounding sphere of 100 simulated reflectances with noise component derived from the covariance matrix. Match errors below this limit are thus within the error that image noise alone may produce.
The computation time for ALUT pixel inversion compares favorably with the two L–M implementations, and the initial ALUT construction time is negligible compared to the full image analysis (Table 1). Due to the BSP tree search algorithm, the ALUT inversion time is dependent on how closely the forward model represents the spectra to be inverted. Hence the ALUT inversion is an order of magnitude faster when inverting synthetic spectra generated from the forward model with random parameter values (Table 1).
Table 1 Per-pixel inversion time from the two L–M implementations and the ALUT algorithm. Algorithm
Levmar L–M Cephes L–M ALUT
Table construction
Per-pixel inversion time
Whole image
Image
Synthetic
Analysis
(min)
(ms)
(ms)
(h)
n/a n/a 6
580 25 1.1
– 28 0.14
– – ~ 24
Inversion times for the image transect pixels and synthetic reflectances generated randomly from the forward model are shown. Hyphens indicate analyses that were not done. The whole image consists of approx. 60 × 106 pixels.
2532
J. Hedley et al. / Remote Sensing of Environment 113 (2009) 2527–2532
Increased inversion efficiency with an improved forward model is a desirable property and implies the ALUT algorithm could attain even better whole-image performance than the results presented here. Faster successive approximation algorithms may exist, while the comparison with two freely available implementations is not definitive the ALUT approach clearly offers competitive inversion times. Beyond efficient inversion, the primary advantage of the ALUT methodology is ease of use. The ALUT construction algorithm has only one user parameter, the number of spectral space points to generate. In fact no decision on this needs to be made, since the hierarchical structure of the ALUT means an existing ALUT can always have more points added to it later on and the final structure will be no different to having built an ALUT of that size in the first place. This is invaluable when using a computationally expensive forward model such as Hydrolight, since a small basic ALUT can be built initially for testing image analysis and then refined later without making any previously modeled spectra redundant. Building a conventional look-up table by regular parameter step sizes may result in an unpredictable total construction time, and early termination will result in a partially constructed LUT of no use for inversion. In contrast, ALUTs are “fully constructed” at all times during construction, so can be built to specific processing time limit, and early termination yields a fully usable ALUT for image inversion. The primary weakness of the ALUT algorithm is that in a worst case scenario the number of points required to describe the shape of the mapping function (Eq. (1)) to a specific level of discretization error in spectral space may rise exponentially as the dimensionality increases. This means the ALUT approach may be inefficient for forward models with a large number of free parameters to be inverted. However, note that it is the dimensionality and space-filling properties in spectral space that are relevant, not the number of parameters per se. A 10-parameter forward model may not necessarily fill out 10-dimensions in spectral space very effectively, due to correlation between spectral bands for example. This potential limitation will therefore be application-specific.
5. Conclusion A generic method that can invert any radiative transfer model that calculates remote sensing reflectance from a series of integer and realvalued parameters has been presented. The method is based on the hierarchical construction of a look-up-table which is post-processed into a binary space partitioning tree. The method has the following key advantages 1) is simple to use with no configuration parameters requiring expert knowledge 2) minimizes the number of forward runs required offering high efficiency for computationally expensive models 3) offers fast inversions, and 4) permits efficient timemanagement by allowing preview analysis before the full table is constructed. These features make the approach highly suitable for operational implementation of image analysis methods. Straightforward and robust inversion allows work to focus on the more important objective of improving the forward radiative transfer model employed.
Acknowledgments This work was funded by the UK's Natural Environment Research Council under grant NE/E015654/1, and also benefited from discussions held during meetings of the Office of Naval Research and Australian Research Council funded Working Group for Optically Shallow Water Remote Sensing. Funding for the Heron Reef CASI dataset and field validation work was provided by Australian Research Council grants to S. Phinn. Image corrections and pre-processing were conducted by Karen Joyce. References Brando, V. E., Anstee, J. M., Wettle, M., Dekker, A. G., Phinn, S. R., & Roelfsema, C. (2009). A physics based retrieval and quality assessment of bathymetry from suboptimal hyperspectral data. Remote Sensing of Environment, 113, 755−770. Buiteveld, H., Hakvoort, J., & Donze, M. (1994). The optical properties of pure water. In J. S. Jaffe (Ed.), Ocean Optics XIIProc. SPIE, Vol. 2258. (pp. 174−183). Darvishzadeh, R., Skidmore, A., Schlerf, M., & Atzberger, C. (2008). Inversion of a radiative transfer model for estimating vegetation lai and chlorophyll in a heterogeneous grassland. Remote Sensing of Environment, 112(5), 2592−2604. Gastellu-Etchegorry, J. P., Gascon, F., & Estève, P. (2003). An interpolation procedure for generalizing a look-up table inversion method. Remote Sensing of Environment, 87 (1), 55−71. Goodman, J., & Ustin, S. L. (2007). Classification of benthic composition in a coral reef environment using spectral unmixing. Journal of Applied Remote Sensing, 1(1), 011501. Hedley, J. (2008). A three-dimensional radiative transfer model for shallow water environments. Optics Express, 16(26), 21887−21902. Joyce, K., 2004. A method for mapping live coral cover using remote sensing. Ph.D. thesis, The University of Queensland. Klonowski, W. M., Fearns, P. R., & Lynch, M. J. (2007). Retrieving key benthic cover types and bathymetry from hyperspectral imagery. Journal of Applied Remote Sensing, 1(1), 011505. Lee, Z., Carder, K., Mobley, C., Steward, R., & Patch, J. (1998). Hyperspectral remote sensing for shallow waters. I. A semianalytical model. Applied Optics, 37(27), 6329−6338. Lee, Z., Carder, K., Mobley, C., Steward, R., & Patch, J. (1999). Hyperspectral remote sensing for shallow waters: 2. Deriving bottom depths and water properties by optimization. Applied Optics, 38(18), 3831−3843. Liang, S. (2007). Recent developments in estimating land surface biogeophysical variables from optical remote sensing. Progress in Physical Geography, 31(5), 501−516. Liang, S., Zhong, B., & Fang, H. (2006). Improved estimation of aerosol optical depth from MODIS imagery over land surfaces. Remote Sensing of Environment, 104(4), 416−425. Manly, B. F. J. (1994). Multivariate statistical methods, 2nd Edition. Chapman and Hall. Mobley, C. D. (1994). Light and Water. San Diego: Academic Press. Mobley, C. D., Sundman, L. K., Davis, C. O., Bowles, J. H., Downes, T. V., Leathers, R. A., Montes, M. J., Bissett, W. P., Kohler, D. D. R., Reid, R. P., Louchard, E. M., & Gleason, A. (2005). Interpretation of hyperspectral remote-sensing imagery by spectrum matching and look-up tables. Applied Optics, 44(17), 3576−3592. Pope, R. M., & Fry, E. S. (1997). Absorption spectrum (380–700 nm) of pure water. II. Integrating cavity measurements. Applied Optics, 36(33), 8710−8723. Wettle, M., & Brando, V. E. (2006). Sambuca: Semi-analytical model for bathymetry, unmixing and concentration assessment. Tech. Rep. CSIRO Land and Water Science Report 22/06 Canberra, Australia: CSIRO. Wolfe, M. A. (1978). Numerical methods for unconstrained optimization. New York: Van Nostrand Reinhold Company. Zhang, Y., Chen, J. M., Miller, J. R., & Noland, T. L. (2008). Leaf chlorophyll content retrieval from airborne hyperspectral remote sensing imagery. Remote Sensing of Environment, 112(7), 3234−3247.
Contents lists available at ScienceDirect
Remote Sensing of Environment j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / r s e
Efficient radiative transfer model inversion for remote sensing applications John Hedley a,⁎, Chris Roelfsema b, Stuart R. Phinn b a b
School of Biosciences, University of Exeter, EX4 4PS, UK Center for Remote Sensing and Spatial Information Science, School of Geography, Planning and Environmental Management, University of Queensland, Brisbane, Australia
a r t i c l e
i n f o
Article history: Received 13 May 2009 Received in revised form 16 July 2009 Accepted 18 July 2009 Keywords: Radiative transfer Inversion Look-up table (LUT) Successive approximation
a b s t r a c t A simple method for efficient inversion of arbitrary radiative transfer models for image analysis is presented. The method operates by representing the shape of the function that maps model parameters to spectral reflectance by an adaptive look-up tree (ALUT) that evenly distributes the discretization error of tabulated reflectances in spectral space. A post-processing step organizes the data into a binary space partitioning tree that facilitates an efficient inversion search algorithm. In an example shallow water remote sensing application, the method performs faster than an implementation of previously published methodology and has the same accuracy in bathymetric retrievals. The method has no user configuration parameters requiring expert knowledge and minimizes the number of forward model runs required, making it highly suitable for routine operational implementation of image analysis methods. For the research community, straightforward and robust inversion allows research to focus on improving the radiative transfer models themselves without the added complication of devising an inversion strategy. © 2009 Elsevier Inc. All rights reserved.
1. Introduction Inversion of physics-based radiative transfer models is an area of rapid development in remote sensing of both aquatic and terrestrial environments (Liang, 2007). In these approaches a parameterized forward model for reflectance such as Hydrolight (for aquatic applications) or PROSPECT, PROSAIL or DART (for terrestrial canopies) takes a series of parameters describing the optical properties of participating media or canopy structure and the model defines a mapping from parameter space to spectral space (Fig. 1a) (Darvishzadeh et al., 2008; Gastellu-Etchegorry et al., 2003; Mobley et al., 2005; Zhang et al., 2008). Image analysis then uses a search algorithm to find the parameter space location which minimizes the distance of the corresponding spectral space location from the image pixel reflectance, measured by a distance metric such as Euclidean distance. Two distinct approaches to the search algorithm can be identified: 1) pre-calculation of reflectance look-up tables (LUTs) by repeated runs of the forward model with systematically differing parameter values (Mobley et al., 2005), optionally with interpolation across tabulated points (Gastellu-Etchegorry et al., 2003; Liang et al., 2006) and 2) spectral matching by successive approximation using optimization techniques such as the Downhill Simplex or Levenberg–Marquardt (L–M) algorithm (Brando et al., 2009; Goodman & Ustin, 2007; Klonowski et al., 2007; Lee et al., 1999; Wolfe, 1978), which is feasible only if the forward model is sufficiently fast to permit many runs
⁎ Corresponding author. E-mail address: [email protected] (J. Hedley). 0034-4257/$ – see front matter © 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.rse.2009.07.008
for each image pixel. With either approach, for a complex radiative transfer model, both the forward modeling and the inversion may take substantial computational effort. Routine or large-scale operational image analysis by physics-based methods therefore demands the development of efficient approaches for both modeling and inversion. Often the forward model is simplified or approximated to facilitate inversion. In contrast, we develop a generic inversion technique which allows the use of any forward model and hence model accuracy need not be compromised. A key problem in the efficient representation of the parameter to reflectance mapping function is that it is not generally clear a priori how to subdivide the parameter space in order to evenly populate spectral space from the forward model. For example, modeled shallow water reflectance typically has an inverse exponential relationship with depth (Mobley, 1994), so a regular subdivision of a depth parameter would over-sample the deep water spectral space region with similar reflectance spectra produced from unnecessary forward model runs. At the same time the shallow water spectral space region would be relatively under-sampled with a higher discretization error in the tabulated reflectances (Fig. 1b). The interaction between multiple parameters (e.g. depth vs. water clarity) makes the general problem of efficient and accurate LUT construction very hard to tackle by analytical means. To address this problem we present an automated construction procedure for adaptive look-up trees (ALUTs) that evenly distributes and minimizes the discretization error in spectral space and facilitates fast inversion (Fig. 1). The construction technique can be applied to any forward radiative transfer model that is parameterized by a series of real and integer values. The following sections describe 1) the requirements of a forward model to which the ALUT construction
2528
J. Hedley et al. / Remote Sensing of Environment 113 (2009) 2527–2532
a set of set of m real-valued parameters, α1, α2…αm, and any number of discrete integer parameters, i1, i2…, etc. As a first step, any discrete valued parameters are combined into a single integer parameter by i = i1 + i2N1 + i3N1N2 + … so that i = 1, 2…N, where N1, N2… are the maximum value each integer parameter can take and N is the product of the maximum values. This gives a set N of mapping functions from the real-valued parameter space to remote sensing reflectance in spectral space, fi : Rm YRn hα 1 ; α 2 N α m iiRrs
Fig. 1. Adaptive LUT construction for a forward model of one continuous parameter, α, mapping to two spectral bands, hb1 ; b2 i. (a) The mapping function from parameter space to spectral space which the LUT seeks to accurately represent. (b) Discretization by regular steps in parameter space over-samples some regions of spectral space and under-samples others. (c) and (d) iterative adaptive construction fills in the gaps to approximate an even distribution of points in spectral space and efficiently represent the function.
ð1Þ
where fi is the forward model for particular discrete value of i and each f i consists of an assumed continuous mapping from mdimensional parameter space to n-dimensional spectral space, with Rrs being the spectral remote sensing reflectance in n-vector form. The separation of discrete and continuous parameters is necessary because the adaptive subdivision can only be applied to the continuous parameters. In essence, a series of separate ALUTs are constructed, one for each value of i, this accommodates forward models that contain discrete parameters, for example where a model component is based on a choice of one from a set reflectance spectra. The only additional requirement is that a priori limits for the continuous parameters are defined, so that α MIN V α i V α MAX , this is the i i same as is required for any LUT methodology and is also often used in successive approximation techniques to constrain the mapping function domain. 2.2. Adaptive look-up-table construction
method can be applied 2) the adaptive construction algorithm, and 3) the methodology for building and searching the tree representation of the ALUT. The subsequent sections present an illustrative example of ALUT construction using a shallow water radiative transfer model (Lee et al., 1999) and its application to hyperspectral airborne imagery of a coral reef. 2. Methods 2.1. Mapping from parameter space to spectral space The adaptive LUT construction can be applied to any forward model for n-band remote sensing reflectance modeled dependent on
ALUT construction begins with an initial regular subdivision of parameter space (as shown for 1-dimensional parameter space in Fig. 1b), then an iterative algorithm repeatedly assesses which regions in parameter space correspond to currently under-sampled regions of spectral space, the most effective parameter is subdivided in that region, the new points are added to the ALUT and the iteration continues (Fig. 1c,d). The key to implementing the subdivision algorithm for mdimensional parameter space is that the ALUT structure is represented as a hierarchical voxelation of the parameter space, where voxels are m-dimensional parameter bounding boxes, the 2m vertices of which are current points in the LUT (Fig. 2). To assess the local point distribution in spectral space, change in Rrs across the voxel is assessed by taking the mean distance in spectral space when changing
Fig. 2. Hierarchical voxelation of parameter space for a forward model of three continuous parameters hα 1 ; α 2 ; α 3 i. Black dots show the parameter space location of points in the ALUT where Rrs has been calculated, open circles show where subdivision is occurring and points need to be calculated for the next iteration.
J. Hedley et al. / Remote Sensing of Environment 113 (2009) 2527–2532
a parameter from its lower bound on one side of a voxel to its upper bound on the other side. This occurs along 2m − 1 edges of the voxel, one edge for each high-low combination of the parameters other than the one being assessed. Using Euclidean distance as the distance metric, the mean change in Rrs for parameter i is
DðiÞ =
1 m−1
2
m−1 2X
jR rs ði; j; 0Þ − Rrs ði; j; 1Þj
ð2Þ
j=1
with Rrs (i, j, 0|1) being the remote sensing reflectance corresponding to the vertices at the low and high parameter value end of the jth voxel edge which describes a change in parameter i. The mean change in Rrs for each parameter, D(i), is calculated for each parameter in each voxel when the voxel is created. The values are stored in a list which is continually updated and sorted during the construction process. Voxels at the top of the list represent the parameter space region which currently has the greatest discretization error in spectral space. The LUT is grown by taking the voxel at the top of the list and bisecting it into two new voxels at the midpoint of the parameter value that represents the greatest mean spectral change across the voxel (Fig. 2). When a voxel is subdivided 2m − 1 new vertices are required, their parameter locations are the midpoints of the subdivided edges, and 2m − 1 forward model calculations are required to generate the corresponding Rrs spectra. The mean spectral space distances for each parameter, D(i), are calculated for the new voxels (Eq. (2)) and added to the sorted list, and the construction procedure continues. To initialize the ALUT construction in m-dimensional space a regular grid of voxels is created with corners at the extremes of the parameter limits, typically this grid can just be a single voxel (Fig. 2a). When the forward model contains a discrete parameter, i.e. N N 1, the ALUT construction algorithm is repeated independently N times, with the forward mapping function being fi, i = 1,2,…N. Essentially N separate ALUTs are constructed, these are merged into a single treestructured LUT as described in the next section. 2.3. Search optimization by binary space partitioning trees
2529
when all tabulated points below the current node have been checked or ruled out. During the search, the Euclidean distance to the current best found Rrs vector is maintained. At each node the algorithm will descend down a branch only if the image reflectance to match is on the same side of the node's dividing plane as the branch, or is closer to the dividing plane than the current best matched Euclidean distance. Only then can points on the far side of the dividing plane possibly be a better match than the current best match. Note that the choice of np, maximum sample size for the PCA, may have a small effect on inversion time but has no effect on accuracy. The suggested value of 30,000 is probably suitable for all applications. 3. Example application The following sections illustrate an example ALUT inversion of a 5-parameter shallow water radiative transfer model applied to a hyperspectral airborne image of a tropical coral reef (Fig. 3). Speed and accuracy of ALUT bathymetric retrievals is compared to two implementations of the Levenberg–Marquardt successive approximation algorithm, with tide-corrected sonar depth transects used for ground-truth. 3.1. Imagery and forward model The adaptive LUT construction and inversion algorithm was applied to a 11278 × 5236 pixel 17-band CASI image (430–710 nm, FWHM 10–20 nm) acquired in 2002 of Heron Reef, Australia (23.44° S, 151.91° E, Fig. 3) with the ALUT constructed to represent a variant of Lee's semianalytic forward model for shallow water remote sensing reflectance (Lee et al., 1998, 1999) where subsurface remote sensing reflectance is estimated from subsurface solar zenith angle (θw =32°) subsurface view angle (θ=0°), and depth in meters, H, by ( "
# )! C 1 Du ðλÞ 1 − exp − + κ ðλÞH cosθw cosθ ( " # ) 1 1 DBu ðλÞ + + ρðλÞexp − κ ðλÞH π cosθw cosθ
dp ðλÞ rrs ðλÞ≈rrs
ð3Þ
To provide a spectral matching search algorithm that is more efficient than an exhaustive search the storage of all the reflectance vectors is post-processed into the structure of a binary space partitioning (BSP) tree. This is achieved by the following algorithm: 1. A list is initialized containing all Rrs vectors in the ALUT linked to their corresponding generating parameter values. 2. A principal components analysis (Manly, 1994) is performed on a random selection of np of the Rrs vectors in the list (example: np = 30,000), or the entire list if the current number of vectors in the list is less than np. 3. A dividing plane in spectral space is established, such that the normal to the plane is the first principal component axis and the plane passes through the mean of the selected Rrs vectors. 4. The list of points is subdivided into two lists dependent on which the side of the dividing plane the Rrs points lie on. 5. The procedure is applied recursively to the two new lists, starting at Step 2, unless the number of points in the remaining lists is below a predetermined limit (example: 20 points) or the tree depth reaches maximum limit (example: 200 subdivisions). In the resulting tree every node has an associated dividing plane and subdivides into two branches, the actual tabulated reflectances are stored in groups of 20 or less in the terminal nodes (leaves). Searching the BSP tree to find the closest match to an image reflectance is efficient since the global optimum can be found without checking every tabulated Rrs vector. The search traverses the tree recursively downwards to the leaves, stepping back up one level
Fig. 3. CASI Image of Heron Reef and location of sonar transects. The full image (top panel) is 11278 × 5236 pixels of 1 m × 1 m in 17 bands 430–719 nm with variable band widths from 10 to 20 nm. The images here are monochrome versions of an RGB image based on bands centered at 665, 564 and 480 nm. The central dark region in (1) is the island itself.
2530
J. Hedley et al. / Remote Sensing of Environment 113 (2009) 2527–2532
where the subsurface remote sensing reflectance for optically deep C B water, rdp rs (λ), and Du(λ) and Du(λ), can be estimated from dp
rrs ðλÞ ≈ ½0:084 + 0:170uðλÞuðλÞ
ð4Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C Du ðλÞ ≈ 1:03 1 + 2:4uðλÞ
ð5Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi B Du ðλÞ ≈ 1:04 1 + 5:4uðλÞ
ð6Þ
tance using the relation R rs = 0:54R 0þ = kQ where Q = 4 (Mobley, 1994) and k is a spectrally independent factor to convert above-sur face to below-surface diffuse reflectance, k ≈ R 0þ = Rð0 − Þ ≈ 0:7, estimated from multiple runs of PlanarRad, an open-source planeparallel radiative transfer model for directional radiance similar to the commercial software Hydrolight (Hedley, 2008; Mobley, 1994). 3.2. ALUT and Levenberg–Marquardt inversion The structure of the forward model to be encoded by the ALUT in the form of Eq. (1) is therefore
where uðλÞ = bb ðλÞ = ½aðλÞ + bb ðλÞ;
κ ðλÞ = aðλÞ + bb ðλÞ
ð7Þ
where a(λ) and bb(λ) are the total absorption and backscattering coefficients and estimated as the sum of the effects of water constituents by, aðλÞ = aw ðλÞ + ½a0 ðλÞ + a1 ðλÞ lnP P + G exp ½−0:015ðλ − 440Þ ð8Þ Y
bb ðλÞ = bbw ðλÞ + X ð400= λÞ
ð9Þ
The pure water absorption and backscatter aw(λ) and bbw(λ) are considered known (Buiteveld et al., 1994; Pope & Fry, 1997), P, G, X and H are free parameters to be solved by inversion, these describe absorption due to phytoplankton, CDOM absorption, particulate backscatter and depth respectively. The particulate backscatter slope is fixed, Y=1, and the phytoplankton spectral absorption parameters, a0(λ) and a1(λ), are tabulated in (Lee et al., 1998). For the substrate cover reflectance in Eq. (3), ρ(λ), a linear mixture model is used(Wettle & Brando, 2006; Brando et al., 2009) where ρðλÞ = Ei1 ðλÞm + Ei2 ðλÞð1 − mÞ is a linear mix, in proportion m, of any two benthic reflectances, Ei1 (λ), Ei2 (λ), taken from a set of thirteen typical reef benthic reflectances including several of sand, live and dead coral, algae and seagrass, so i1 = 2…13 and i2 = 1… i1 − 1 giving 78 pair combinations, describable by a single integer parameter, i = 1…78. The forward model outputs (Eq. (3)) were translated to vectors of above-surface remote sensing reflectance, Rrs using Rrs ð jÞ = 0:5rrs λj = 1 − 1:5rrs λj where λj is the CASI band j wavelength center (Lee et al., 1999). The CASI image had been previously corrected to units of above-surface diffuse reflectance, R(0+) Joyce (2004) and so was converted to above-surface remote sensing reflec-
fi : R5 YR17 hP; G; X; H; miiRrs for i = 1 N 78
ð10Þ
The lower and upper limits of the parameters used to establish the single initial parameter voxel (Fig. 2a) were: 0 ≤ P ≤ 1.5, 0 ≤ G ≤ 0.08, 0 ≤ X ≤ 0.007, 0 ≤ H ≤ 20, 0 ≤ m ≤ 1. Hence the maximum retrievable depth is 20 m, the ranges of G and X were chosen for relatively clear waters while the upper bound of P, the component of a(440) due to phytoplankton, was set substantially higher than would normally be expected in these waters to illustrate features of the ALUT construction (Fig. 4). Two ALUTs of 5 × 106 spectral space points were generated from the above parameterization, this size of ALUT being the maximum the computer's memory could accommodate. The first, henceforth called “Automatic ALUT”, followed the previously described procedure exactly with automatic subdivision of parameters based on the distribution of points in spectral space. However, the strong effect of water absorption on reflectance above 700 nm results in an extreme treatment of depth parameter subdivision: shallow depths are generated with steps as fine as 5 mm, while other parts of the ALUT contain no depths between 10 and 20 m (Fig. 4a). This is an accurate picture of the effect of depth on the spectral shape and magnitude of reflectance, but raises the question as to whether accuracy could be improved by a more practical discretization of the depth parameter. To answer this, a second ALUT was constructed where depth subdivision was forced to a minimum of five subdivisions and constrained to eight, so throughout the “Detailed Depth ALUT” depth is at worst represented to an accuracy of 0.63 m all the way to 20 m, but no finer than 0.08 m in shallowest regions (Fig. 4b).
Fig. 4. Example slices through (a) the automatically generated ALUT and (b) ALUT generated with constraints on the depth parameter subdivision. The diagrams show two axes in parameter space, depth and the phytoplankton parameter, P. The values for other parameters are fixed at G = 0, X = 0, m = 0 and the bottom reflectance is pure sand. Boxes represent cross-sections of voxels (Fig. 2) and line intersections are the point locations in the ALUT where Rrs has been calculated.
J. Hedley et al. / Remote Sensing of Environment 113 (2009) 2527–2532
2531
Fig. 5. Depth estimates versus sonar derived bathymetry, n = 5057. Dashed lines are the regression slope plus and minus one standard error, Δy is the mean of the absolute value of the errors between estimated and sonar bathymetry.
To compare the ALUT inversion approach to alternate methods the model was also inverted by two publicly available C code implementations of the L–M algorithm, the Netlib Cephes package1 and Levmar2.
4. Results and discussion The Levmar L–M optimization and both ALUT inversions resulted in almost identical accuracies for bathymetric retrieval assessed over 5057 sonar data points (Fig. 5). In all three cases r2 = 0.91, the slight difference between the L–M and ALUT inversions with respect to regression slopes (1.078 vs. ~0.93) and in the mean of absolute depth error (0.90 vs. ~0.77) can be partly accounted for by the lack of an upper depth parameter constraint on the L–M inversion. Constraining ALUT inversion to depths less than 20 m is essentially the incorporation of a priori information, since it was known all sonar depths were less than 20 m. However, incorporating the same upper bounds on the L–M inversion actually produced worse results, due to algorithm convergence issues in the vicinity of the constraints. Since the ALUT BSP tree search algorithm returns the global optimum it is free from convergence issues. While the L–M successive approximation was able to produce slightly closer spectral matches to the image reflectances than the ALUT (Fig. 6) this did not translate into higher bathymetric accuracy. Both methods were able to obtain spectral matches close to or within image noise levels, estimated from the covariance matrix of pixels over deep homogeneous water (Fig. 6). Hence the ability to produce a spectral match was not limiting for either approach. The Automatic ALUT inversion displayed clear discretization of depth retrievals (Fig. 5b) due to the depth parameter subdivision structure of the ALUT (Fig. 4a). However while the Detailed Depth ALUT produced a cosmetically superior scatterplot (Fig. 5c), accuracy assessed by the three metrics of r2, regression slope and mean absolute error was virtually identical to the Automatic ALUT (Fig. 5b). The distribution of spectral match errors was also virtually identical (Fig. 6). This indicates that the automatically generated structure of the ALUT does give an efficient representation of the models capability to retrieve environmental parameters. The very similar spread of the L–M and ALUT bathymetric retrievals points to a deficiency in the underlying forward radiative transfer model or image noise as the ultimate accuracy limiting factor. 1 2
URL: http://www.netlib.org/cephes. URL: http://www.ics.forth.gr/˜lourakis/levmar.
Fig. 6. Spectral matching errors over pixels containing sonar data points, expressed as distribution of Euclidean distance between image and matched spectral reflectances, n = 4146. Environmental and image noise is estimated as the noise-equivalent deltareflectance, NEΔRrs, by the covariance of spectral reflectance over a homogeneous deep water region, and then by calculating the radius of the bounding sphere of 100 simulated reflectances with noise component derived from the covariance matrix. Match errors below this limit are thus within the error that image noise alone may produce.
The computation time for ALUT pixel inversion compares favorably with the two L–M implementations, and the initial ALUT construction time is negligible compared to the full image analysis (Table 1). Due to the BSP tree search algorithm, the ALUT inversion time is dependent on how closely the forward model represents the spectra to be inverted. Hence the ALUT inversion is an order of magnitude faster when inverting synthetic spectra generated from the forward model with random parameter values (Table 1).
Table 1 Per-pixel inversion time from the two L–M implementations and the ALUT algorithm. Algorithm
Levmar L–M Cephes L–M ALUT
Table construction
Per-pixel inversion time
Whole image
Image
Synthetic
Analysis
(min)
(ms)
(ms)
(h)
n/a n/a 6
580 25 1.1
– 28 0.14
– – ~ 24
Inversion times for the image transect pixels and synthetic reflectances generated randomly from the forward model are shown. Hyphens indicate analyses that were not done. The whole image consists of approx. 60 × 106 pixels.
2532
J. Hedley et al. / Remote Sensing of Environment 113 (2009) 2527–2532
Increased inversion efficiency with an improved forward model is a desirable property and implies the ALUT algorithm could attain even better whole-image performance than the results presented here. Faster successive approximation algorithms may exist, while the comparison with two freely available implementations is not definitive the ALUT approach clearly offers competitive inversion times. Beyond efficient inversion, the primary advantage of the ALUT methodology is ease of use. The ALUT construction algorithm has only one user parameter, the number of spectral space points to generate. In fact no decision on this needs to be made, since the hierarchical structure of the ALUT means an existing ALUT can always have more points added to it later on and the final structure will be no different to having built an ALUT of that size in the first place. This is invaluable when using a computationally expensive forward model such as Hydrolight, since a small basic ALUT can be built initially for testing image analysis and then refined later without making any previously modeled spectra redundant. Building a conventional look-up table by regular parameter step sizes may result in an unpredictable total construction time, and early termination will result in a partially constructed LUT of no use for inversion. In contrast, ALUTs are “fully constructed” at all times during construction, so can be built to specific processing time limit, and early termination yields a fully usable ALUT for image inversion. The primary weakness of the ALUT algorithm is that in a worst case scenario the number of points required to describe the shape of the mapping function (Eq. (1)) to a specific level of discretization error in spectral space may rise exponentially as the dimensionality increases. This means the ALUT approach may be inefficient for forward models with a large number of free parameters to be inverted. However, note that it is the dimensionality and space-filling properties in spectral space that are relevant, not the number of parameters per se. A 10-parameter forward model may not necessarily fill out 10-dimensions in spectral space very effectively, due to correlation between spectral bands for example. This potential limitation will therefore be application-specific.
5. Conclusion A generic method that can invert any radiative transfer model that calculates remote sensing reflectance from a series of integer and realvalued parameters has been presented. The method is based on the hierarchical construction of a look-up-table which is post-processed into a binary space partitioning tree. The method has the following key advantages 1) is simple to use with no configuration parameters requiring expert knowledge 2) minimizes the number of forward runs required offering high efficiency for computationally expensive models 3) offers fast inversions, and 4) permits efficient timemanagement by allowing preview analysis before the full table is constructed. These features make the approach highly suitable for operational implementation of image analysis methods. Straightforward and robust inversion allows work to focus on the more important objective of improving the forward radiative transfer model employed.
Acknowledgments This work was funded by the UK's Natural Environment Research Council under grant NE/E015654/1, and also benefited from discussions held during meetings of the Office of Naval Research and Australian Research Council funded Working Group for Optically Shallow Water Remote Sensing. Funding for the Heron Reef CASI dataset and field validation work was provided by Australian Research Council grants to S. Phinn. Image corrections and pre-processing were conducted by Karen Joyce. References Brando, V. E., Anstee, J. M., Wettle, M., Dekker, A. G., Phinn, S. R., & Roelfsema, C. (2009). A physics based retrieval and quality assessment of bathymetry from suboptimal hyperspectral data. Remote Sensing of Environment, 113, 755−770. Buiteveld, H., Hakvoort, J., & Donze, M. (1994). The optical properties of pure water. In J. S. Jaffe (Ed.), Ocean Optics XIIProc. SPIE, Vol. 2258. (pp. 174−183). Darvishzadeh, R., Skidmore, A., Schlerf, M., & Atzberger, C. (2008). Inversion of a radiative transfer model for estimating vegetation lai and chlorophyll in a heterogeneous grassland. Remote Sensing of Environment, 112(5), 2592−2604. Gastellu-Etchegorry, J. P., Gascon, F., & Estève, P. (2003). An interpolation procedure for generalizing a look-up table inversion method. Remote Sensing of Environment, 87 (1), 55−71. Goodman, J., & Ustin, S. L. (2007). Classification of benthic composition in a coral reef environment using spectral unmixing. Journal of Applied Remote Sensing, 1(1), 011501. Hedley, J. (2008). A three-dimensional radiative transfer model for shallow water environments. Optics Express, 16(26), 21887−21902. Joyce, K., 2004. A method for mapping live coral cover using remote sensing. Ph.D. thesis, The University of Queensland. Klonowski, W. M., Fearns, P. R., & Lynch, M. J. (2007). Retrieving key benthic cover types and bathymetry from hyperspectral imagery. Journal of Applied Remote Sensing, 1(1), 011505. Lee, Z., Carder, K., Mobley, C., Steward, R., & Patch, J. (1998). Hyperspectral remote sensing for shallow waters. I. A semianalytical model. Applied Optics, 37(27), 6329−6338. Lee, Z., Carder, K., Mobley, C., Steward, R., & Patch, J. (1999). Hyperspectral remote sensing for shallow waters: 2. Deriving bottom depths and water properties by optimization. Applied Optics, 38(18), 3831−3843. Liang, S. (2007). Recent developments in estimating land surface biogeophysical variables from optical remote sensing. Progress in Physical Geography, 31(5), 501−516. Liang, S., Zhong, B., & Fang, H. (2006). Improved estimation of aerosol optical depth from MODIS imagery over land surfaces. Remote Sensing of Environment, 104(4), 416−425. Manly, B. F. J. (1994). Multivariate statistical methods, 2nd Edition. Chapman and Hall. Mobley, C. D. (1994). Light and Water. San Diego: Academic Press. Mobley, C. D., Sundman, L. K., Davis, C. O., Bowles, J. H., Downes, T. V., Leathers, R. A., Montes, M. J., Bissett, W. P., Kohler, D. D. R., Reid, R. P., Louchard, E. M., & Gleason, A. (2005). Interpretation of hyperspectral remote-sensing imagery by spectrum matching and look-up tables. Applied Optics, 44(17), 3576−3592. Pope, R. M., & Fry, E. S. (1997). Absorption spectrum (380–700 nm) of pure water. II. Integrating cavity measurements. Applied Optics, 36(33), 8710−8723. Wettle, M., & Brando, V. E. (2006). Sambuca: Semi-analytical model for bathymetry, unmixing and concentration assessment. Tech. Rep. CSIRO Land and Water Science Report 22/06 Canberra, Australia: CSIRO. Wolfe, M. A. (1978). Numerical methods for unconstrained optimization. New York: Van Nostrand Reinhold Company. Zhang, Y., Chen, J. M., Miller, J. R., & Noland, T. L. (2008). Leaf chlorophyll content retrieval from airborne hyperspectral remote sensing imagery. Remote Sensing of Environment, 112(7), 3234−3247.