1. Introduction

Two main directions determine science and technology applications of optical remote sensing data [1]: pattern recognition of the land surface objects on the images under processing and state parameter estimates for the selected object classes [2]. Photo pictures enable to interpret textures of the relevant objects. Multispectral and hyperspectral remote sensing systems are evolved towards using enhanced number of their spectral channels that serves to analyze a thin structure of the classes. Advances in optics and photonics facilitated remote sensing systems called hyperspectral imaging or imaging spectroscopy [3–5 ] as one of the perspective directions of remote sensing monitoring of the objects represented on an image as patterns to be recognized. This application area is considered as an extension of that given by the common-used digital cameras and has become a new tool of the land surface parameters inference from remote sensing images [6].

The characteristic feature of hyperspectral systems consists in the ability to increase the information content of registered data as compared with the multispectral systems. On the other hand, produced for the enhanced spectral characterization, the hyperspectral systems may have redundant channels for particular scenes due to the correlation between data in the neighboring channels. Thus a necessity emerged to reduce this redundancy by appropriate computational techniques [7]. The standard approach to the classification problem can be improved by preliminary feature optimization procedure [6]. Besides that, the optimization based on textural analysis implies finding the likelihood between any registered set of data and the theoretical distributions as well as to regularize the solution by employing the derivative functions characterizing the neighborhood of the pixels for the related classes [8].

Pattern recognition as a scientific discipline passed the way from the first techniques of optical imagery processing [9] to the updated techniques of computer vision to understand the nature of the phenomena generating the image in a particular subject area [10]. Any scene under processing is understood in the contextual constraints of the objects on it, which are recognized by object features at the lower level of abstraction [11]. The consideration of the image distortion problem due to the atmosphere as scattering and absorbing media was also among the computational approaches of imagery processing [12]. The fusion of these two disciplines has resulted in practical applications of remote sensing hyperspectral imagery processing for monitoring forest species [13], peat and forest fires identification [14], and for many other environmental problems.

The initial statements of the optical vision deal with a corrupted image recovery [15], texture analysis [16], perceptual grouping [17], object matching and recognition [18], pattern mining [19], retrieving the land surface parameters [20] using the atmospheric correction procedures [21]. The listed disciplines have created a basis for the pattern recognition of objects on remotely sensed images.

Machine-learning algorithms of pattern recognition have become indispensible means of remote sensing imagery processing. There are at least three stages of the available applications of optical data processing: 1) “a visual decoding” of separate airborne or space-borne pictures by the eye as an optical instrument of an experienced operator [22]; 2) to represent digital images in the form of such combinations of the spectral channels like Normalized Difference Vegetation Index (NDVI) and similar other intuitive products [23]; 3) to gain an automation of imagery processing including optimized sets of hyperspectral channels based on the machine-learning algorithms with cross-validation and independent estimates of the accuracy [24].

In this paper, we present some results of airborne hyperspectral imagery processing with the aim to retrieve detailed forest structure attributes and productivity parameters. Instead of the limited standard approach based on the vegetation indices concept [25, 26 ], we employed inverse modeling technique, modified Bayesian classification taking into account analysis of texture and spectral patterns and the optimization of the feature space robust to small changes in training data.

2. Measurements and instrumentation

The campaign of airborne and ground-based measurements was carry out in the summer of 2011 for the test areas in the Savvatyevskoe forestry (Tver, Russia). The purpose of ground-based measurements was the correction and widening of the available standard forest inventory data and collecting referencing data for the studied objects: areas with different types of terrain, water objects, vegetation associations, and anthropogenic zones. The airborne measurements covered the area of the size about 4x10 km by 13 overlapped hyperspectral images. The most part of the test area is covered by the forest stands of different ages mainly of the pine and birch species. The detailed description of data of this campaign can be found in [6, 10 ].

The airborne measuring system illustrated below consists of the hyperspectral camera (HSC), the power module, the notebook with the corresponding software, the aircraft gyro-stabilized platform and the professional photo camera for obtaining simultaneously high resolution images of the test area in the nadir direction. The technical equipment passed an obligatory preflight check and routine maintenance between the flights. The used HSC operating in visible and near-infrared (VNIR) spectral ranges is produced by Science and Production Enterprise (SPE) “Lepton” (Zelenograd, Moscow, Russia). The external appearance of the hyperspectrometer is given by Figs. 1(a) and 1(b) . Basic technical parameters can be summarized as follows:

 

Fig. 1 Airborne hyperspectral imaging spectrometer (SPE “Lepton”, Russia). a) – external appearance, b) – the spectrometer installed on the gyro-stabilized platform during the airborne survey, c) –schematic representation.

Download Full Size | PDF

The schematic composition together with the parameters of basic components of HSC can be seen in the Fig. 1(c). An image formed by the input objective is truncated by the spectral slit and come to the dispersive device. The dispersing element is the complex direct-vision prism consisting of 7 optically glued simple prisms produced by SPE Lepton. The prism supplies the linearity of the dispersion curve in VNIR spectral range and the co-linearity of the direction of the mean wavelength light beam propagation to the instrument axis. This makes images to be rectangular for all spectral channels. The output objective projects the dispersed image to the photo-receiver. As a result, each row of the Charge Coupled Device (CCD) matrix provides spectral radiances corresponding to any pixel of hyperspectral image. Available dark channels allow us to investigate the characteristics of the intrinsic noise of the CCD matrix.

The width of spectral channels depends on the wavelength. Thus, for the equalization of sensitivity of the instrument, some of channels should be joined into pseudo channels. Details of this procedure can be found in [6]. The nominal spatial resolution is near to 1.1 m at the typical flight altitudes 2 km. The pixel size along the trajectory of the flight depends on the speed of the airplane and changes within 0.66-0.91 m. Having sufficiently small size and weight, the HCS can be potentially used for sensing from unmanned aerial vehicles.

3. Recognition of forest species and productivity parameters retrieval in optical imagery processing

We are following here the main ideas of the realization of optical imagery processing given in [2, 6, 10 ]. The process consists in two basic stages: recognition of the forest stand attributes and retrieval of productivity parameters.

For the first stage we proposed an improvement of the standard Bayesian classification principle which allows us to join the texture and spectral analysis for the solution of the recognition problem. The texture analysis procedures serve to identify the likelihood between any registered set of data and the theoretical distributions as well as to regularize the solution by employing the derivative functions characterizing the neighborhood of the pixels for the related classes. As the result we obtain estimates of prior probabilities for basic classes in the considered scene. The spectral analysis uses the Gaussian mixture model to approximate the real probability distribution of spectral features of the classified objects. The detailed description of the algorithm and corresponded training data for the considered test area is given by [10].

The high dimensionality of the spectral feature space is typical for hyperspectral data processing. Thus, at the fixed number of samples of the training set, the known curse of dimensionality problem dealing with instability of the training process usually appears. Instability means that even small changes in the training set lead to significant response of estimates of classifier parameters. The proposed classifier also subject to this problem since the Gaussian mixture model used contains the inverse covariance matrices for each of the classes. The solution of this problem demands the effective reduction of the feature space. In the paper [6] we proposed the new stable algorithm of the selection of the most informative channels which allows reducing redundancy of the spectral channels in the particular subject area. Here we employed the sequences of spectral channels indicated in [6].

Figure 2 gives some results of the recognition of the forest species composition based on processing the hyperspectral image for the selected area with prevailing pine and birch species. The map of the ground-based forest inventory is represented by Fig. 2(a) in the form of colors: orange – for prevailing pine species, blue – for prevailing birch species (in both cases the more the age of tree’s stands is, the darker color is for these plots), the plots with horizontal lines represent the logging places. The quarters on Fig. 2(a) are outlined by solid black lines and are denoted by bold numbers as it is typical for the ground-based forest inventory discipline. The plots within the quarters are reproduced by black dashed lines with normal numbers. Typical for each plot are the following four numbers: a conditional number of the plot, the average age of the forest stands, square meters of the plot, its wood quality (the more the number of ground-based forest inventory data, the less quality of the forest wood).

 

Fig. 2 Recognition of the forest tree species for the selected region with 7 plots: a) – the forest inventory map: orange – pine dominance, blue – birch dominance, the yellow frame indicates the location of the test area; b) – the photo image of the test area; c) – results of the recognition of the forest species composition taking into account the levels of the Sun illumination of the forest canopy (1 – shadowed, 2 – semi-sunlit, 3 – sunlit), magenta – unrecognized objects.

Download Full Size | PDF

Figure 2(b) reveals the picture of the area from the air-photo instrument installed on the same gyro-stabilized platform of the airborne carrier as the imaging spectrometer. Figure 2(c) shows the results of the forest pattern recognition using the imaging spectrometer data processing for the area represented by the yellow frame of Figs. 2(a) and 2(b). At least, 7 plots belonging to the prevailing pine species (plots 1, 3-5, 7) and birch species (plots 2, 6) are seen on Fig. 2(c). We can see the pure pine plots denoted as 10P (the species content (P – pine, B – birch) with the resolution of 10 percent) and the mixed forest plots (90% of the birch stands and 10% of the pine stands, both plots denoted as 9B1P) given by the pixel-by-pixel recognition procedure. We can also see the spatial distribution of pixels characterized by the tree’s categories of the Sun illumination conditions in accordance with three colors of the related pixels (1 – shadowed, 2 – semi-sunlit, 3 – sunlit).

The total weighted error amounts to 16.6%, for shadowed pixels it is 18.5%, for partly illuminated – 18.0% and for completely illuminated – 16.2%. We should have in mind that the accuracy of ground-based data for the forest species composition amounts to 5-20% taking into account a natural level of the errors due to the boundary pixels. We can thus conclude that we are able to make the forest species automation recognition using the airborne imaging spectrometer with the accuracy comparable to that given by the laborious inventory procedures.

Once the forest classes are recognized by remote sensing imaging spectrometer data processing, the next stage is to retrieve parameters characterizing the forest canopies of different species and ages. The green phytomass amount of leaves/needles and the wood biomass amount are among these parameters.

Projective characteristics of the forest canopies for different species and ages near the nadir angles of viewing are of major interest in our studies, once the main task is to process remote sensing images of high spectral and spatial resolution. Direct problem of atmospheric optics is solved to represent the radiances registered by the imaging spectrometer in terms of the incoming spectral flux (irradiance) for the canopy model from the direct solar radiation and the diffuse scattered flux from the sky under the related Sun zenith angle, the initial flux at the top of the atmosphere, the atmospheric transmittance function for direct solar radiation as well as the modeling inter-crown, intra-crown and multiple scattering inside the canopy of the radiation, which fall down within the view of a measuring instrument, and have the corresponding spectral reflectivities and the shading parts between the crowns and on the crowns, respectively [24]. In this multi-parametric approach, we represent the spectral radiances normalized on their integral values depending on the canopy density $D_{canopy}$ and tree’s crown density $D_{crown}$ for the related species $L_{norm}(D_{canopy}{D}_{crown})$. The product $D_{canopy}{D}_{crown}$ characterizes the projective cover of the forest canopy, and we assume once the 3D structure of the canopy is unacceptable using data only from imaging spectrometer, that the phytomass volume B of leaves/needles can be retrieved by the projective cover characteristics $B=B(D_{canopy},D_{crown})$. The retrieval procedures represent solutions of the inverse problem of atmospheric optics.

Any solution of the direct problem of atmospheric optics implies calculations of the modeling approach in multi-dimensional space of the indicated sets of parameters. The inverse problem of retrieving parameters of the forest canopy may not have a single solution and needs regularization of the ill-posed problem once the calculation database does not describe all variations of the spectra registered. A certain type of spatial resolution reduction to a particular scale is needed for the airborne hyperspectral data processing since under the resolution near to 1 m each tree stand can be represented by several pixels. Therefore, the parameter of the canopy density cannot be defined for a single pixel, and we have to use data of the direct model starting from the characteristic scale near to 10 m. Indeed, this scale must be larger (10-20 m), having in mind the swath width of the airborne imaging spectrometer that is near to 500 m.

In Fig. 3 , we show the main stages of data processing. Pattern recognition using spectral radiances for pixels of high spatial resolution is within the processing procedure. The reduction of the spatial resolution to the scales 10-20 m is conducted both for initial hyperspectral data (the binary classification – forests and remaining objects) and for the recognition results due to the necessity to determine a pure species contribution to the mixed forests. The reduction procedure consists in two stages: the possible mixed forest calculation and selection of the spectral channels used. The first procedure serves to characterize the mixed forest in the lower resolution representation that together with the direct problem of pure species enables to construct a transitional direct model for the mixed forest. The second procedure leads to finding the radiances of the lower resolution. Thus, the inverse problem is solved and the output information products are obtained from these two operational fluxes.

 

Fig. 3 Block diagram of the inverse problem module.

Download Full Size | PDF

The inverse problem of the forest stand parameters retrieval is solved for pixels referred on the hyperspectral images to the forest classes. The number of options in the direct problem is given by the following parameters: Sun elevations for different seasons, phyto-elements shading, end-member reflectivities of the sunlit crowns, intra- and inter-crown shades, particular forest species, number of the spectral channels used, etc. The densities of crowns and canopy of any plot are the main parameters to be retrieved. The steepest descent method known in mathematics is used to replace the pure running over techniques by this approach to enhance the efficiency of the computational treatments.

Partial solutions of the inverse problem are found while searching any correspondence between the measured spectra and modeling results of the direct problem. Figure 4 reproduces these techniques. Modeling radiances in each spectral channel represented by the relevant surface are subjected to the cross-sections by a plane on the level of the radiances measured for the pair of spectral channels (blue and red color at Fig. 4). The partial solution is given by a point or an area (green color) where the particular cross-section takes place between the curves in the coordinates given by the canopy density and tree’s crown density, each such point corresponds to the phytomass amount of forest vegetation.

 

Fig. 4 Scheme of searching a partial solution of the inverse problem.

Download Full Size | PDF

The phytomass implies for grasslands and agricultural crops that the leaves and stems can be mowed for a particular area (for example, of the 2x2 m size) and represented in terms of raw materials and dry weight in the laboratory. Forest researchers have their own estimates of the phytomass amounts (units: tonnes per hectar) using also the laboratory dry weight parameters of leaves/needles and crown spreading for a particular plot. Special formulas are used to re-calculate these phytomass estimates to the wood content (the total biomass amount) including tree’s trunks and the age of the forest canopy on a plot. The situation is as if all available leaves/needles are gathered for laboratory testing once we dwell on this model of the outgoing radiation formation and the forest parameters retrieval.

Additionally to the phytomass amount, various parameterizations are used to retrieve other parameters such as the absorbed photosynthetically active radiation (APAR), pure carbon production. APAR values are given by the wavelengths 400-700 nm taken by vegetation patterns from incoming solar fluxes during the vegetation growth, and only a small portion of these values goes on the photosynthesis and carbon sequestering.

The listed retrieval procedures of remote sensing hyperspectral data processing can be used for enhanced parameterization of forested environments in climate models. This is another case of the applications concerning terrestrial ecosystems function and bio-geographical consideration of forest problems [27]. The common-used applications in this particular subject area imply finding estimates from remotely sensed data relating to biological productivity of different types of forests. Net Primary Productivity (NPP) characterizing the income of carbon for a particular season or larger intervals of time as a part of the APAR fraction (fAPAR) extraction is the major factor of the forest growth [28]. The NPP values are measured in units of grams of carbon for an area per a year. There are models of the NPP formation originating from energy and water exchange between a canopy and the atmosphere including calculations of the evapotranspiration rates (a joint effect of the water vapor evaporation and transpiration of living plants).

Figure 5 depicts a test area with prevailing deciduous (birch, different tones of blue color) and coniferous (pine, different tones of orange color) species. We can see at Fig. 5(a) the forested area, roads, a settlement on the left side, etc. We can also see at Fig. 5(b) the quarters of forest inventory and some plots within them, logging places (plots with horizontal lines), etc. As usual, the darker the color is, the more is the age of the related stands within the selected plots.

 

Fig. 5 A test area of thematic interpretation of airborne hyperspectral data processing: RGB-synthesized image – a), the forest inventory map of the test area – b), the results of the recognition (dark-green –birch species, cyan – pine species, light-green – grass, black – roads, other objects) – c) and the forest species composition divided on several classes under coarse resolution (from pure conifer species – red color to pure deciduous species – blue color with the intermediate colors relating to the mixed forests) – d).

Download Full Size | PDF

Figure 5(c) shows the recognition results as well as the spatial distribution of the forest species (from coniferous – red color to deciduous – blue color), specially reduced to coarse resolution of 10 m for visual demonstration and interpretation of the direct modeling results Fig. 5(d). The high level of correspondence between these ground-based and remotely sensed illustrations is feasible.

Next stage is in searching the partial solutions of the inverse problem in accordance with Fig. 4 scheme. Data from all pure species are transferred to the transitional model and optimal parameters of the crowns density and the canopy density are obtained. The projective characteristics retrieval on this test site is given by Fig. 6 . The quantitative characteristics of the projective cover area are from the values near to 0.007 (red color) to values near to 0.97 (green color) at Fig. 6(a). Low projective cover values are prevailing in the lower part of this picture where in accordance with Fig. 5(b) data both types of coniferous and deciduous species are seen. The density of the canopy is represented by the values from 0.016 (red color) to 0.995 (green color). The intermediate and low values of these characteristics are prevailing at the lower part of this picture (Fig. 6(b)) while high values are typical at its upper part. The density of crowns is within (0.07-0.98), respectively, but again low values of this quantity (red color) are typical in the lower part of this picture and higher values at its upper part (Fig. 6(c)). Intuitively, we can feel from Fig. 6 images that the forests in the lower parts of these images are different (red colors are prevailing) as compared with their upper parts (green colors are prevailing). Not so much depending on the species composition (see Fig. 5), we can ascertain that the current state of forests in the lower parts of Fig. 6 images is worse than in their upper parts. This fact may be due to the soil moisture effects (the moisture deficit or contrary its surplus), enhanced air pollution as a result of proximity of the lower terrain to the road, etc.

 

Fig. 6 The results of the inverse problem solution represented as the projected cover area – a), cover density – b), crowns density – c).

Download Full Size | PDF

The estimates of the biological productivity of the forest stands are conducted using the related parameterizations depending on the species composition and projective characteristics. This discipline has had its relevant pre-requisites, but dry phytomass values of leaves/needles obtained by the laboratory treatments of the related samples have been maintained to be the most sensitive to the parameters of the crowns density and canopy density. The results of calculation of the leaves/needles phytomass for the test area are given by Fig. 7(a) . The highest phytomass corresponds to the coniferous species. These values can be seen to be low at the lower part of this picture (red color) while intermediate values are apparent at its upper part (green color). Similar characterization for NPP values (Fig. 7(b)) shows their low values for the lower part of the test area and their higher values for its upper part. These values are important for construction of carbon cycle models. In fact, we can see that the NDVI values are not sensible to the carbon cycle changes (Fig. 7(c)).

 

Fig. 7 Retrieval results of biological productivity of the forest stands: the volume of dry leaves phytomass – a), NPP values – b), NDVI values – c).

Download Full Size | PDF

4. Conclusion

Optical remote sensing data processing problem is highlighted from the point of view of machine-learning algorithms and optimization procedures in computer vision. Main attention is paid to the hyperspectral imagery applications to recognize the forest attributes as well as to retrieve the productivity parameters for the related classes. The proposed technique allows the simultaneous use of texture and spectral features, the stability of training is provided by an improved feature selection algorithm.

The retrieval results of the forest productivity parameters are based on solving the direct problem of atmospheric optics and the inverse problem of these parameters retrieval using hyperspectral imagery processing and a modeling scheme of searching partial solutions of the ill-posed problem. The inverse modeling employs the links between projective characteristics and such parameters as the biomass of the stand fractions and the net primary production. The technique is alternative to the standard approaches based on vegetation indices.

The future improvements can be obtained from the fusion of the imaging spectrometer installation on the same gyro-stabilized airborne platform with the lidar, which provides obtaining 3D structure of the forest canopy using both these passive and active instruments [29]. Such technique was already applied in practice for estimation of the above ground biomass in tropical forest [30] using the restricted vegetation indices concept. Dealing with wider computational procedures of the recognition and parameters retrieval while hyperspectral imagery processing, our approach serves to extend the relevant application areas.