Abstract
The recent development of Earth observation satellites with multiangular capabilities and enhanced spectral resolution has led to preliminary attempts at determining the height of atmospheric scatterers, in particular, of top-cloud heights and smoke plumes originating from forest fires. Inspired by these previous studies, the present work presents an original methodology for the determination of the three-dimensional distribution of high-contrast atmospheric aerosols using multiangular images. The method starts with the approximately known geometry of image acquisition and a set of tie points and uses a linearized and regularized functional model to obtain the position of atmospheric scatterers identified by means of a semiassisted procedure on two or more images. A subsequent application to a CHRIS/PROBA-1 scene of Mount Etna following its eruption on June 14, 2014, allows determining the volcanic plume three-dimensional structure with a precision in the 100–200 m level.
1. Introduction
In the last two decades, Earth observation programmes have experienced a substantial development with the usage of sensors including not only increasing geometric resolution and number of radiometric bands but also multiangular capabilities. Among the missions and sensors providing multiangular images, one may find the Advanced Along-Track Scanning Radiometer (AATSR) onboard ENVISAT (active until 2012) and its predecessor ATSR-2, which provided images with 0° and +55° along-track angles in seven radiometric bands (from 550 nm to 1200 nm) to a spatial resolution of 1 km. The Multiangle Imaging Spectroradiometer (MISR) onboard TERRA/SAR is capable of registering 0°, 26.1°, 45.6°, 60.0°, and 70.5° images in four spectral bands of the visible and near infrared with a spatial resolution of 275 m at nadir. Further, onboard the currently most recent mission, Sentinel-3 (launched in February 2016), the sensor Sea and Land Surface Temperature Radiometer (SLSTR) is designed to collect 0° and −55° images with a spatial resolution of 300 m in the highest resolution mode. All of these sensors, however, are overcome in spatial resolution and number of radiometric bands by the Compact High Resolution Imaging Spectrometer (CHRIS) onboard ESA’s PROBA-1 mission. CHRIS acquires images to a spatial resolution of 17 m (in acquisition mode 5) or 34 m (modes 1 to 4) in a large number of bands (62 in mode 1, 18 in modes 2, 3, and 4, and 37 in mode 5) for the spectral range of 400–1050 nm within nominal along-track angles of 0°, ±36°, and ±55° (see [1]).
The most significant benefit obtained from the use of multiangular images versus the traditional use of only nadiral images derives from a better characterization of the spectral response of the surface appearing in the scene without the use of additional prior hypothesis [2–4]. Some few works, however, use this multiangular capability to retrieve information on the vertical distribution of scatterers in the atmosphere, albeit with little detailed results (i.e., poor vertical resolution) and just for significantly dense scatterers. Hence, Horváth and Davies [5] and Moroney et al. [6] determined top-cloud heights, and Kahn et al. [7] obtained heights of smoke plumes due to forest fires. Theoretical radiative transfer and diffusion models have also been the topic of recent research (e.g., [8, 9]). Since knowledge of the distribution of aerosols in the atmosphere has been pointed out by the Intergovernmental Panel on Climate Change (IPCC) as one of the key aspects for a better understanding of the climate evolution mechanisms [10], it seems that research on the application of remote sensing imagery to this question must be fostered. There is no doubt that satellite imagery has been increasingly used in recent years for aerosol retrieval; however, its main application has remained in the context of determination of aerosol general parameters (e.g., Aerosol Optical Depths) having left the determination of vertical distribution, with very little exceptions, to other techniques like LIDAR from the CALIPSO (Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation) programme [11, 12].
In this paper, we present a methodology for the determination of the three-dimensional structure of high-contrast aerosols in the atmosphere using multiangular images of the CHRIS/PROBA-1 sensor along with its application to the determination of the vertical distribution of Mount Etna’s volcanic plume. In the following section, we explain the data preparation and the functional and stochastic models that are input to the adjustment by the least-squares method. Then we will present and discuss the results and draw the corresponding conclusions, including the limitations of the proposed methodology and possible lines for future research.
2. Materials and Methods
2.1. Data Preparation
We use the VISAT-BEAM open-source toolbox [13] that ESA released for public analysis, processing, and exploitation of their products to work with the CHRIS images of Mount Etna dated June 28, 2014—subsequent to its eruption on June 14, 2014—which are available at the ESA download site (ESA Earth Observation Users’ Single Sign On, https://eo-sso-idp.eo.esa.int/idp/umsso20/admin), where the volcanic plume is still clearly visible. The image files show the scenes corresponding to along-track nominal angles of 0°, +36°, +55°, and −55° (the image file corresponding to −36° is missing). Further, it has to be noted that the −55° scene contains little information useful for the determination of the three-dimensional structure of the plume: a small portion of the plume with too little features whose homologous points could be identified in the other scenes; therefore, we decided not to include it in the analysis. A note of caution is worth here regarding the angular value: +55°, +36°, and 0° are only the nominal along-track angle values of CHRIS/PROBA, while the real ones may be actually different due to the across-track pointing of the platform. While we will be continuously referring to +55°, +36°, and 0° scenes, these values cannot be entered in the subsequent computations as true values of the image geometry, which indeed involves along-track as well as across-track angular values. At any rate, knowledge of the actual values for these angles will not be needed for our calculation.
For every image, its metadata include an estimation of the image center time (here 5:06:42, 5:05: 55, and 5:05:07 for 0°, +36°, and +55° scenes, resp.). We use these times to interpolate in the corresponding telemetry file that can be retrieved from ESA’s REDU Center (http://194.78.233.110/products/data/CHRIS_Additional_data/) and compute the WGS84 Cartesian coordinates for every image center.
CHRIS images usually suffer from a misalignment problem of unknown origin, which makes the image centers differ considerably from their targeted positions [4, 14]. Besides the inclusion into VISAT-BEAM of the corresponding telemetry file, some ground control points (GCPs) also have to be supplied in order to correct image misalignments. Apart from these, some additional ground points with known coordinates were used to evaluate the quality of the resulting transformation. We used Google Earth to retrieve the terrain coordinates of all ground points. We found the resulting typical errors to be around 100 m for the 0° and 36° images, whereas 200 to 300 m typical errors were found for the 55° image. This is likely to be caused by greater point identification difficulties in a very inclined image. As an example, the distribution of the selected GCPs over the 55° scene, and their corresponding coordinates, is given in Figure 1 and Table 1.
Next we identified homologous points in the different images of the plume using one radiometric band, where the contrast is high (band 26, 682.2 nm). It is worth noting that other studies of volcanic plumes have used spectral bands of very similar wavelengths (e.g., [15]). Zooming in and out over the images by means of the VISAT-BEAM, we can be quite accurate in identifying some homologous points in the different images and obtaining their coordinates, especially for the plume outline. For other parts in the plume, however, the identification of homologous points turned out to be almost impossible. Figures 2–4 show general layouts of the set of homologous points (accurate identification was done in zoomed-in views).
The larger limitation in identifying homologous points over the 55° image is self-evident: some of the homologous points in the two other images are out of scene now; other points were included albeit with a limited degree of certainty about their correctness, while some others were more reliable.
It has to be noted that these coordinates correspond to the plume projection onto the terrain for every image. As appearing in Figure 5 in the next subsection, these coordinates will be denoted with lowercase letters , , and .
2.2. Functional and Stochastic Models
For the sake of simplicity, we derive now the functional model based on two images. Its generalization to three or more images, to be used later, is straightforward.
After interpolation in the telemetry file, the satellite positions for the two scene shots are assumed to be known (at least approximately) in a global geodetic system like the WGS84. A right-hand set of Cartesian coordinates is defined in this system with -axis in direction of the (conventional) Earth rotation axis and -axis and -axis in the equatorial plane (with -axis pointing along the direction of 0° longitude). Transformation from geodetic coordinates latitude φ, longitude λ, and height over the reference ellipsoid to Cartesian coordinates x, y, and z is obtained applying the well-known formulae [16–18]where ν is the curvature radius of the ellipsoid in the prime vertical:
and is the semimajor axis of the ellipsoid, whose eccentricity is denoted by . The relationship between ellipsoidal height and orthometric height (i.e., above the geoid or mean sea level) is also straightforward,
through knowledge of the geoid undulation N, here assumed to be constant for the entire study area m as taken from the EGM2008 geoid [19].
For a particular scatterer in the atmosphere of unknown coordinates (), we have the projected point onto the terrain of coordinates () for image 1 and () for image 2. These coordinates are known after direct measurement of image coordinates and terrain georeferencing under VISAT-BEAM and later transformation from geodetic coordinates to Cartesian coordinates with the formulas above. We also have the corresponding satellite coordinates () and () obtained after interpolation using the telemetry file, and the coordinates of the scene target (). The corresponding geometry is shown in Figure 5.
Projective rays fulfill the straight-line parametric equations,
for every image with the respective parameter . The coordinates of the scatterer () are common to both sets of equations and therefore can be determined with two or more images.
For simplicity, we can work with the first of these equations only and denote it after a slight rearrangement as
with being the functional model of the problem at hand. This is an equation that is exactly fulfilled with the exact values of satellite coordinate in scene , , straight-line parameter , terrain coordinate , and scatterer coordinate . We could be naïve enough to think that we truly know three of these parameters (the satellite position is known by means of the telemetry, the image coordinate is measured, and parameter of the straight line is defined through the coordinates of these two points belonging to the line) so that there is only one unknown, coordinate . However, these three parameters are not known with total accuracy so that, if simply fixed, the inaccuracies in their values would be transferred to inaccuracies in the final coordinate with no possible way to separate the errors corresponding to each of the error sources in addition to causing an incontrollable bias in the solution. We must acknowledge, however, that they have been obtained only approximately: the satellite coordinate was computed after interpolation for the image central time; terrain coordinate is the result of a measurement and it is affected by the inevitable random error that is present in every measurement process (including also the previous measurements needed for the image georeferencing); and, finally, parameter can be roughly determined taking into account the approximate scene inclination (0°, 36°, or 55°) or more accurately from the use of these (approximate) satellite and terrain coordinates. To remedy this, we will use system regularization that accounts for the estimated error size in the (inaccurately) known parameters and has the effect of fixing the known parameters albeit with a loose margin, that is, not as exact values but as statistical values each with a given most probable value and a given standard deviation so that the inaccuracies in the approximate parameters are not transferred to the final solution.
Consider, first, the fact that the functional model (see (5)) is not linear in all variables; second, consider the fact that approximate values for all these variables may be obtained (they will be denoted as , , , and ), and, third, consider that the disparate absolute degree of accuracy of these approximate values will require the system to be regularized; linearization of the model is therefore necessary. After first-degree Taylor’s expansion, the functional model readswhere partial derivatives are to be evaluated for the approximate values, , , and denote corrections to the corresponding approximate values , , and , and denotes the residual to be determined after the adjustment for the observation . After linearization, we have transformed the problem into the determination of the corrections δ, , and to approximate values , , and as well as the residual of the observed value . We will then obtain the final adjusted values as
Computing the derivatives and rearranging the equation to leave the unknowns in the left-hand side, except for the residual term, we can write
Analogous expressions can be obtained for the and coordinates. For two images (denoted with subscripts 1 and 2), considering the coordinates of several scatterers (denoted with superscripts , , etc.), we can finally write the corresponding system of equations in matrix form as
withwhere the columns refer, respectively, to the unknowns δ, δ, δ, δ, δ, δ, , , , , , , and so forth. That is, columns 12 to 16 repeat the structure of columns 7 to 11 now for scatterer number 2 and so on. Similarly, rows are generated in blocks of 6 for every scatterer.
Vectors , , and have the following structure:
Rectangular matrix and column vector are made of numeric values that can be computed from the approximate or measured values, whereas column vectors (observation residuals) and (correction to approximate values) are the unknowns to be determined in the adjustment. Before regularization, to be explained later on, we have, for scatterers, equations with unknowns for the case of a model based on 2 images, equations with unknowns for a model based on 3 images, and so forth.
The resulting Cartesian coordinates can be transformed back to geodetic coordinates by making use of the well-known formulas [17],
in terms of the following auxiliary variables and the curvature radius already defined in (2): where is the semiminor axis of the ellipsoid. We can therefore obtain the ellipsoidal height or, by means of (3), the height above mean sea level .
Now the system of equations in (9) has to be adjusted by least-squares method, since we assume the measurands to be statistical variables following normal distributions. A reasonable a priori estimate for the precision of each terrain coordinate can be 100 m for the 0° and 36° images and 200 m for the 55° image, in accordance with the georeferencing analysis previously mentioned. This also goes along the lines of Alonso and Moreno [14] who found georeferencing errors to be around 100 m for flat areas and up to 300 m for abrupt reliefs. Once defined in this fashion, the standard error σ of each measured terrain coordinate, a diagonal weight matrix W (with dimension equal to the number of rows in A) is created with weights in its diagonal, with each of them defined as
Further, the system of equations in (9) has unknowns of very different nature for which we expect very different values: corrections to the approximate straight-line parameters, , and so forth; satellite coordinates, δ, and so forth; and position of the scatterers, , and so forth. The application of parameter regularization to the system of equations is interesting then, as some other authors [3] have already pointed out for the case of multiangular images.
We pay attention now to the computation of the different approximate values and their corresponding degrees of accuracy, which will define the figures to use in the regularization. First, we have assumed that all pixels of each image are simultaneously acquired; obviously, this is only a working assumption, especially since CHRIS/PROBA uses a nonstandard push-broom acquisition technique with a nonuniform rotation speed of the platform in order to maximize image quality. However, we can explain the validity of this simplified hypothesis as a working assumption for the final 3D reconstruction of the plume as follows. Taking into account the fact that the time difference with respect to the central time in each scene can amount to 10 s [14] and the satellite is orbiting the Earth 600 km above its surface [1] with an orbital period of 100 min [1], the differences in the satellite positions during the image acquisition can amount to 70 km. Previous studies indicate that there are already differences of the order of several km between satellite coordinates obtained by different techniques [14]. All in all, as they say, these discrepancies usually have a negligible effect on the final solution. As an example, we can reckon that an error of 70 km in the satellite position may represent in the worst case only some 115 m of error in the position of a scatterer located 1 km above the Earth’s surface. We therefore estimate as 50 km the typical error of the initial approximate coordinates of the satellite (assumed to be static during acquisition) and understand that the impact on the final results (typically of the order of tens of meters) is acceptable for the definition of the plume 3D structure.
Another working assumption is the static character assumed for the plume during the successive image acquisitions. Obviously, the effect of the wind modifies the shape of the plume, although it can be argued that the relative displacement between points belonging to the plume will be more than one degree of magnitude less than the total wind speed. For instance, a wind speed of 10 km/h (2.8 m/s) represents some 130 m of wind displacement between two consecutive images acquired with some 47 s time separation, surely leading to a relative displacement of points in the plume below the level of some 10–20 m. There are also other intrinsic movements in the plume, due to thermal turbulent convection, which would be difficult to accurately model; however, we can also assume them to be below the level of a few tens of meters for our time intervals. A simulation using a mesoscale numerical model [20] agrees with the order of magnitude of this velocity.
Regarding the coordinates of a particular scatterer in the atmosphere, we can use as a first approximation their planimetric coordinates in the nominal along-track 0° image along with a reasonable average height (e.g., the target height). We estimate the typical corrections to this value as some 500 m.
Finally, straight-line parameters can be computed from the scene shot geometry or, more precisely, in terms of the satellite and scatterer approximate coordinates. Their values result from the order of unity and the corresponding corrections to be obtained may be of the order of 0.05.
We will therefore perform the regularization of the system of equations by adding one pseudoobservation equation in matrix A for every unknown. These regularization equations will be of the type (in this example for coordinate )
and will be accompanied by the corresponding weight in the diagonal matrix W:where = 50000 m is the estimated error for satellite coordinates, as mentioned previously.
The resulting system of equations is solved by least-squares method as
The a posteriori unit weight standard deviation is obtained [21] aswhere and are, respectively, the number of equations and the number of unknowns. The proximity of the value to unity is an indicator of the quality of the adjustment (correctness of functional and stochastic models). We can also apply the data snooping test [21] to determine the homologous points that were incorrectly identified.
3. Results and Discussion
After adjustment of the model based on two images (0° and 36°), we obtain residuals typically below 100 m with a few exceptions for points 140, 142, 144, 145, and 146 reaching 150–200 m, which can be regarded as acceptable results (they are consistent with the prior precision). The a posteriori unit weight standard deviation is 0.914, and its closeness to unity confirms the validity of the adjustment. Iteration of the computations produces significant changes in neither the residuals nor the coordinates of the scatterers.
Regarding the adjustment of the model based on three images (0°, 36°, and 55°), we start by noticing an undesired result: the a posteriori unit weight standard deviation is now 2.998, which indicates that the existing errors are 3 times, on average, the expected value. By inspecting the corresponding residuals, we find out that some of them have been affected by gross errors. They correspond to some of the doubtful observations we made over the third image (55°). Starting from the largest residual (corresponding to point 142 on image 3), we iteratively eliminate the point with the worst residual in the successive adjustments as it is customary in the data snooping procedure until all the remaining residuals are acceptable. The a posteriori unit weight standard deviation is now 1.870 and the adjustment can be accepted. We observe, however, that the final residuals are mostly below 100 m for the 0° and 36° images, whereas they have a typical size of 400–500 m for the 55° images. This leads us to conclude that the error estimation for 0° and 36° images (σ = 100 m) was correct, whereas the error estimation for the 55° image (σ = 200 m) was too optimistic. At any rate, the resulting coordinates of the scatterers do not differ much from those obtained with the two-image model so that both adjustments can be regarded as valid. These differences are mostly below 60 m, being smaller in the vertical components than in the horizontal ones, amounting to some 200 m for the worst case. It all can be regarded as a sensible result if we take into account the different uncertainties involved in the model.
The three-dimensional distribution of the particular scatterers defined in the plume contour is shown over a Google Earth view in Figure 6. Although these scatterers constitute only a small sample of the volcanic plume, they offer the possibility to characterize its three-dimensional distribution including height, shape, and direction (towards southeast).
As a final comment, we can recall that the best observational geometry for the determination of coordinates by intersection is the one that is made at a right angle in the unknown object [22]. This means that, to optimally define the vertical structure of the plume by using only two images, we would need images at 45° and −45°. If more than two images are to be acquired, the corresponding angles should be more or less evenly separated from this value on both sides. The design of the CHRIS/PROBA-1 nominal angles (with 45° almost halfway from 36° and 55°) closely conforms to this idea and must be regarded as nearly optimal, at least from the purely geometrical point of view, obviating the fact that identification of homologous points is more difficult the larger the view angle difference is.
4. Conclusions
An original methodology for the determination of the three-dimensional distribution of high-contrast atmospheric scatterers using multiangular images was derived and subsequently applied to the case of Mount Etna’s volcanic plume from high-resolution CHRIS/PROBA-1 images. The final results proved to be reliable within the 100–200 m level.
The methodology presents, however, some limitations in its different steps. Regarding data acquisition, a high degree of contrast of aerosols in the atmosphere is needed to perform the necessary identification of homologous points in the different inclination scenes. The correct identification of these homologous points is a critical factor limiting the precision of the final results. Gross error testing and identification are therefore unavoidable. Another limiting factor comes from the relatively small sample of identifiable features appearing on all the different views of the plume. These limited numbers of homologous points may suffice, however, to roughly define the three-dimensional distribution of the plume. With regard to the different simplifications in the computations, we note the assumption that the plume does not move or change its shape during the total acquisition time, as well as the simplification that for every image all pixels have been acquired in a common time. Both may lead to errors assumed to be within the indicated total error budget (100–200 m). Finally, from the algebraic point of view, we may note that the problem is ill-conditioned and requires the use of regularization. It is expected that application to other cases where more images are available (the five of them or at least images on both sides of the nadir view) produces better results than those obtained here only with 0°, +36°, and +55° images and would allow drawing conclusions on the measuring precision achievable using different sets of angles.
Further research may include application to other areas and types of dense aerosols (e.g., wildfire plumes), automated measuring of homologous points over the different images, inclusion of a more detailed characterization of the platform movement during the acquisition, and the use of robust estimation to detect incorrect measurements.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.