Retrieval of forest stand attributes using optical airborne remote sensing data Vladimir V. Kozoderov,1,* Timofei V. Kondranin,2 Egor V. Dmitriev,2,3 and Anton A. Sokolov4 2
1 M.V. Lomonosov Moscow State University, Leninskiye Gory, GSP-1, Moscow, 119991, Russia Moscow Institute of Physics and Technology (State University), Institutskiy per., 9, Dolgoprudny, Moscow Region, 141700, Russia 3 Institute of Numerical Mathematics of Russian Academy of Science, ul. Gubkina, 8, Moscow, Moscow, 119333, Russia 4 Laboratoire de PhysicoChimie de l'Atmosphère, Université du Littoral Côte d’Opale & Université Lille Nord de France, Avenue Maurice Schumann, 189A, Dunkerque, 59140, France * [email protected]
Abstract: Optical remote sensing data processing is proposed for the airborne images of high spectral and spatial resolution. Optimization techniques are undertaken to gain information about spatial distribution of the pixels on the hyperspectral images and the texture of the forest stands of different species and ages together with reducing redundancy of the spectral channels used. The category of neighborhood of pixels for particular forest classes and the step up method of selecting optimal spectral channels are employed in the relevant processing procedures. We present examples of pattern recognition for the forests as a result of separating pixels, which characterize the sunlit tops, shaded space and intermediate cases of the Sun illumination conditions on the hyperspectral images. ©2014 Optical Society of America OCIS codes: (280.1415) Biological sensing and sensors; (280.4788) Optical sensing and sensors; (280.4991) Passive remote sensing.
References and links 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16.
S. Z. Li, Markov Random Field Modeling in Computer Vision (Springer-Verlag, 1995). V. V. Kozoderov and E. V. Dmitriev, “Remote sensing of soils and vegetation: pattern recognition and forest stand structure assessment,” Int. J. Remote Sens. 32(3), 5699–5717 (2011). A. K. Jain, “Advances in mathematical models in image processing,” Proc. IEEE 69(5), 502–528 (1981). N. S. Friedland and A. Rosenfeld, “Compact object recognition using energy-based optimization,” IEEE Trans. Pattern Anal. Mach. Intell. 14(7), 770–777 (1992). B. Tso and R. C. Olsen, “A contextual classification scheme based on MRF model with improved parameter estimation and multiscale fuzzy line process,” Remote Sens. Environ. 97(1), 127–136 (2005). J. Bolton and P. Gader, “Random set framework for context-based classification with hyperspectral imagery,” IEEE Trans. Geosci. Rem. Sens. 47(11), 3810–3821 (2009). V. V. Kozoderov, “Application of optical remote sensing data to study natural and climatic processes,” Climate and Nature 2(3), 3–16 (2012). F. G. Hall, D. E. Knapp, and K. F. Huemmrich, “Physically based classification and satellite mapping of biophysical characteristic in the southern boreal forest,” J. Geophys. Res. 102(D24), 29567–29580 (1997). T. Hakala, J. Suomalainen, S. Kaasalainen, and Y. Chen, “Full waveform hyperspectral LiDAR for terrestrial laser scanning,” Opt. Express 20(7), 7119–7127 (2012). E. V. Dmitriev and V. V. Kozoderov, “Optimization of spectral bands for hyperspectral remote sensing of forest vegetation,” Proc. SPIE 8887, 888705 (2013). J. Suomalainen, T. Hakala, H. Kaartinen, E. Räikkönen, and S. Kaasalainen, “Demonstration of a virtual active hyperspectral LiDAR in automated point cloud classification,” ISPRS J. Photogramm. Remote Sens. 66(5), 637– 641 (2011). M. Dalponte, L. Bruzzone, and D. Gianelle, “Fusion of hyperspectral and LIDAR remote sensing data for classification of complex forest areas,” IEEE Trans. Geosci. Rem. Sens. 113, 1416–1427 (2009). V. V. Kozoderov and E. V. Dmitriev, “Remote sensing of soils and vegetation: regional aspects,” Int. J. Remote Sens. 29(9), 2733–2748 (2008). V. V. Kozoderov, T. V. Kondranin, E. V. Dmitriev, and V. P. Kamentsev, “Mapping forest and peat fires using hyperspectral airborne remote-sensing data,” Izv., Atmos. Ocean. Phys. 48(9), 941–948 (2012). J. Besag and P. J. Green, “Spatial statistics and Bayesian computation,” J. R. Stat. Soc., B 55(1), 25–37 (1993). M. Bertero, T. A. Poggio, and V. Torre, “Ill-posed problems in early vision,” Proc. IEEE 76(8), 869–889 (1988).
#210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15410
17. V. V. Kozoderov, E. V. Dmitriev, and V. P. Kamentsev, “A system of airborne remote sensing data processing of high spectral and spatial resolution,” Issledovanie Zemli Iz Kosmosa. 6, 57–64 (2013). 18. V. V. Kozoderov, T. V. Kondranin, and E. V. Dmitriev, Thematic processing of multispectral and hyperspectral airspace images (MIPT Publication, 2013). 19. S. Z. Li, “Object recognition from range data prior to segmentation,” Image Vis. Comput. 10(8), 566–576 (1992). 20. R. Kohavi, “A study of cross-validation and bootstrap for accuracy estimation and model selection,” Int. Joint Conf. on Artificial Intell. (IJCAI), 528–535 (1995). 21. T. Hastie, R. Tibshirani, and J. Friedman, The Elements of Statistical Learning (Springer, 2001). 22. M. Dalponte, L. Bruzzone, L. Vescovo, and D. Gianelle, “The role of spatial resolution and classifier complexity in the analysis of hyperspectral images of forest areas,” Remote Sens. Environ. 113(11), 2345–2355 (2009). 23. V. V. Kozoderov, T. V. Kondranin, and E. V. Dmitriev, “ An apparatus and programmatic system of hyperspectral airspace imagery processing,” Proc. Int. Symp. Atm. Radiation and Dynamics (ISARD), 41 (2013).
1. Introduction The vision theory of the optical airborne data processing implies an optimal computational procedure realization to retrieve parameters from registered radiances for the land surface objects called “patterns” [1]. Pattern recognition techniques are a valuable approach to the processing of remote sensing data for images of the land surface – atmosphere system [2]. The related cognitive technologies originate from the mathematical models of artificial intelligence and data mining [3]. This discipline has passed through the computational procedures of the objects classification on the processed images [4] to the updated techniques of the contextual recognition of the objects as fuzzy sets [5] and the random set framework [6]. The contextual recognition for the forested environments assumes that the neighboring pixels are taken into account for the particular classes of forests of different species and ages [7] since these classes are formed by the “end-members” (sunlit forest canopy, sunlit background and shaded background for a particular view and solar illumination angle) [8]. Each pixel views collections of the individual forest trees and the pixel-level reflectance can thus be computed as a linear mixture of sunlit tree tops, sunlit background (or understory) and shadows. Applications of optical remote sensing include the use of airborne photo-pictures, hyperspectral images and the laser scanning (lidar, light detection and ranging) data. The hyperspectral imagery processing enables to gain information about the distribution of spectral and spatial characteristics of the forest stands of different species composition and age classes. Data from lidars and imaging spectrometers with Charge Coupled Devices (CCD) serve to combine the priorities of the active and passive remote sensing instruments to carry out both these types of measurements from the same airborne platform [9]. The idea can be realized by using high-productive computers and optimization techniques of a neighborhood measure extraction for selected forest classes and of accounting for the optimal set of the spectral channels, thus reducing the possible redundancy of the entire hyperspectral data sets [10]. The combined activity contributes to obtaining the 3D products for the forest canopy [11, 12]. The related mathematical treatments of using different types of classifiers combined with the parameter retrieval procedures for coniferous and deciduous forest classes are maintained to have much wider applications as compared with the image features and object shapes extraction, which relates to photometry and geometry in pixel-level reflectance representation of the forested land cover [13]. The pixel fraction and reflectance of the end-members are only a part in the listed techniques. Instead of the photometry and geometry constraints, the improved models are developed of the functional description of outgoing spectral radiation, in which such parameters of the forest canopy like the vegetation biomass density for particular forest species and ages are embedded [14]. This permits us to calculate the relationships between the registered radiances and the forest biomass densities (the direct problem of atmospheric optics). The next stage is to find solutions of this problem as cross-sections of the related curves in the multi-dimensional space given by the parameters of these models (the inverse problem).
#210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15411
The typical solutions of the inverse problem may not be mathematically unique and the computational procedure is undertaken to their regularization by finding minima of the functional called “the energy for the particular class of forests”. In general, the energy category is considered as a linking chain between the image pixels as random fields and our understanding what class of objects each pixel belongs to, in the known Bayesian paradigm of statistical decision making [15]. The optimization 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 [16]. This article presents some results of the forest species and age recognition for a test area using the listed mathematical framework for ensembles of airborne data from an imaging spectrometer produced in Russia operating in near to 200 spectral channels from the wavelengths 400 nm to 1000 nm. Our main purpose was to arrange a mutual airborne and ground-based field campaigns on a test area with geo-referencing pixels of the hyperspectral images to the standard description of the relevant “pure stands” with the ultimate goal to answer the question, whether we can have a practical tool of the pattern recognition of the forest species and age classes or not. Our prior assumption was as follows. Should we recognize the texture of the young, mature and old forest that is perceived as an alternation of the pixels characterizing the sunlit tops and the shaded background by the computer vision techniques, and we can find such a tool within perceptual grouping of the pixels using their spectral features. As usual, a confusion matrix will enable to compare the predicted ages and the ages given by a particular map of the forest inventory for the test area. 2. Instruments and processing methods The measuring system used for obtaining simultaneous photo and hyperspectral high resolution images of the test area in the nadir direction consists of the hyperspectral camera (HSC), the power module, the notebook with the available software, the gyro-stabilized platform and the professional photo camera. The light-weight airborne prismatic HSC is produced and calibrated by Science and Production Enterprise “Lepton” (Zelenograd, Moscow, Russia). The complex direct-vision prism consisting of 7 optically glued simple prisms supplies the linearity of the dispersion curve in the visible and near infrared (VNIR) region 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 basic parameters of HSC are as follows (Fig. 1). The camera has 287 active and 3 noise channels in the spectral range 401-1017 nm. The spectral resolution of the channels significantly depends on the central wavelengths and changes from 0.36 nm to 14 nm. The radiometric resolution of the digitized signals is 212 − 31 = 4065 levels. The central wave lengths are presented at Fig. 1(a). The main channels belong to the visible region. The typical feature of such prism instruments is in their different spectral resolution for the short-wave and long-wave regions. In particular, the width of the channels increases more than on order in the near infrared region (Fig. 1(b)). The measuring system passed an obligatory preflight check and routine maintenance between the flights. The instrument is calibrated by the Enterprise “Lepton” every year. Before the flight, this imaging spectrometer takes special tests to study the linearity of the present calibration and the relationship between signal and noise. An optimal storage time, an effective number of gray-scale levels and spectral channels are maintained during each flight.
#210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15412
Fig. 1. The central wavelengths (a) and the spectral resolution (b) of the imaging spectrometer.
The algorithms used for the retrieval of forest stand attributes are based on the Bayesian classifier, which is optimal in the sense of minimization of the misclassification probability. According to the Bayesian theory of decision making, when both the prior distribution and the likelihood function of a pattern are known, the best that can be estimated from these sources of knowledge is given by this classifier. The learning process of the Bayesian classification algorithm consists in the estimation of prior probabilities of the classes and corresponding likelihood functions. The prior probabilities are obtained from results of texture analysis in terms of Markov random fields (MRF). We follow here the account of this problem given in [1] concerning the labeling problem in computer vision with the relevant contextual constraints and the neighborhood system for pixels of a particular class of objects. MRF theory provides a consistent way of modeling context depending entities such as image pixels. This is achieved through characterizing mutual influences among such entities using MRF probabilities. The local properties of MRFs lead to parallel computational algorithms. The conditional probability density functions of features are modeled in terms of Gaussian mixtures. The known techniques are considered to be hard to implement for hyperspectral data processing because the number of learning samples necessary for providing the stable probability density function estimate increases too fast with the number of features. In practice, even if we have a set of learning samples large for all considered classes, the learning process would be very expensive in computational terms. The maximum a posteriori (MAP) solution is typical in this framework. The MAP-MRF labeling formalism assumes that a search of minima of a specified function called the energy of the image scene or of the particular class of objects is undertaken [17]. The Bayesian classification principle is applicable for different cases of optical image processing [18]. The classification schemes are based on the assumption that similar ground objects have similar spectral and texture characteristics. The shape of registered spectra depends on many factors and may have significant variations even for objects which seem to be homogeneous. The reflected radiation spectra can significantly change according to the volume distribution of leaves and branches, the state of leaves and underlying surface even in the case when such basic attributes of forest stands as species composition, ages and heights are the same. Modeling of spectral reflectivity of vegetation and other ground objects demands the availability of many parameters which cannot be easily obtained in practice. Applicability of the Bayesian classifiers for processing of airborne images of high spatial resolution (1-2 m) is restricted by the fact that the imaging spectrometer producers can ensure very high signal-to-noise ratio for the pixels relating to the completely illuminated by the Sun tree’s crowns. Meanwhile, the ratio for the background pixels is not so high. The MAP-MRF labeling concerns a supposition of the pixels distribution function. If we have an image or its separate segments which can be described by a constant radiance function for all pixels, then we deal with the known in mathematics string problem (the membrane in the two-dimensional #210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15413
case) [19]. The energy for this class of objects with very low variety of radiances (or the biological diversity of vegetation) has its zero minimum and is positive otherwise. If the prior supposition is in constant gradient of the registered radiances for the land surface, then we deal the rod problem, which energy is represented by the integral from the second derivative function to the second power. It is possible to make a prior supposition about the constant curvature for the class of objects under consideration and to obtain the formula for the energy in the form of the integral from the third order derivative function to the second power, etc. The posteriori energy minimization for any class of objects results in smoothing down the differences between the radiances for the listed three categories of the related pixels (sunlit tops, backgrounds and the intermediate cases of partly illuminated and partly shaded endmembers). Additionally to these smoothing constraints, we have elaborated a rigorous approach to optimize spectral channels of HSC based on searching their most informative sets. The high dimensionality of the feature space at the fixed number of samples of the learning set is typical for hyperspectral data processing. It makes the learning process unstable because of the known problem of the “curse of dimensionality”. The parametric Bayesian classifiers also suffer from this problem to a lesser extent than the nonparametric ones. One of the frequently used techniques for the feature space optimization is the stepwise forward selection method. For the classification algorithm we can apply this method as follows. At each step the features are subdivided into two groups: features already included in the classifier and other features which could be potentially included. After that, features from the second group are added one-by-one to the classifier and corresponding classification errors are estimated. The feature corresponding to the maximum decrease of error is selected. If the error decrease is statistically significant, the feature is included in the classifier and the process continues. Otherwise, the process is stopped. The standard form of the stepwise forward selection has the following two known disadvantages. At first, the classification errors are estimated in the dependent manner because the test set of samples is used for learning the classifier. Secondly, the optimal sequence of the features obtained by this method is highly sensitive to small variations of the learning set. Especially, it concerns the last members of the sequence. We propose to use a modified method allowing us to avoid these problems. The scheme of the method is demonstrated by Fig. 2(a). The stages presented at Fig. 2(a) serve to find the most probable sequence in the forward selection process by the series of launches to estimate the classification error and finding the most frequent features on the different levels of their extraction from the ensembles of the selected features. This calculation process is necessary to select the most informative features, which lead to the minimal classification error. The number of the learning samples for each of the classes is often large for the estimation of the classification error using the holdout cross-validation method [20, 21]. In this method, the initial set of samples is randomly, with the given probability divided on two subsets – training and testing. We used the probability 0.5 and thus the subsets obtained are almost equal. The classifier is always learned using the training subset and the classification errors are estimated in independent manner only for the testing subset. Thus, employing this method for the stepwise forward selection, we avoid the overtraining problem that is more important than the decrease of accuracy in the related statistical estimates due to the reduction of the number of training samples. If the number of the learning samples is not sufficiently large, it is possible to use the similar algorithm based on the leave-one-out cross-validation [22], however this method is much more expensive in computational terms. Random partitioning of the learning data set using the holdout cross-validation is repeated many times. Then the stepwise forward selection is launched for each pair of the training and testing sets obtained. As a result, we have a locally optimal sequence of features. The sequences obtained as a rule are different in composition and length. The unique regularized solution of the feature selection problem can be obtained as the most probable sequence.
#210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15414
Fig. 2. Optimization of the feature space: a) – scheme of the selection of the subset of the most informative features; b) – an example of the selection of the most probable sequence of features. Data corresponding to different levels are highlighted by the grey scale.
In the standard technique, the leading feature corresponds to the minimum total probability of misclassification. However, it must be taken into account that the dissimilarity of the classified objects can be very different. In this case, the choice of the leading feature is mostly defined by the presence of similar objects corresponding to high misclassification errors. On the other hand, some well-distinguished objects, for instance, water, soil and vegetation areas can be classified exactly using just one spectral channel. The channels allowing the exact classification can be not unique and change for different groups. Thus, in our study we define the leading feature as the spectral channel allowing the exact classification of the maximum number of different objects available in the learning set. The algorithm of finding the most probable sequence of features is demonstrated by Fig. 2(b). At the beginning (level 0), we form the set of the members of the sequences coming
#210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15415
after the leading feature and calculate the mode. Then, we form a subset of the sequences corresponding to the mode (level 1) and repeat the process with the next members. The process continues (levels 2, 3, 4) up to selecting a unique sequence. In the case of a few equivalent sequences, the random choice is used. The regularized solution is found to be much less sensitive to variations of the learning samples. The robustness of this method increases with the number of the launches. As a result of the listed optimization constraints, we can select the information layers for a particular class of forests, which pixels belong to the sunlit tops, shaded space and intermediate cases of the Sun illumination conditions on the hyperspectral images. Depending on the species and age of the forests, these pixels for the same class have different spectral distributions and the information layers serve to improve the recognition accuracy. 3. Results and discussion The described techniques of hyperspectral imagery processing were tested on data of the airborne campaign conducted in Tver region of Russia. The test region of the size about 4x10 km was encompassed by 13 overlapped images obtained from the airplane equipped on the same gyro-stabilized platform by the imaging spectrometer and photo-camera. The location of the region and corresponding tracks of the flights are represented in Fig. 3. Each hyperspectral image contains 500 pixels across the track and approximately from 10000 to 14000 pixels along the track. The spatial resolution across the track is stable and amounts to 1.1 m at the flight altitude about 2 km above the ground level. The pixel size along the track depends on the flight speed and changes within 0.66-0.91 m. The landscape of the test region contains different types of natural and artificial objects: the most easily distinguishable objects in accordance with their spectral and texture features are water bodies, open soils, asphalt roads, buildings, forests and grasslands. As we can see from Fig. 3, about half of the area is covered by the forest stands for which the ground-based forest inventory maps are available. Besides that, an extensive field campaign was conducted on the test area to arrange its geo-botanical, forest typology and other descriptions for georeferencing the selected “pure stands” to make the validation and ground truth estimates. Analyzing the forest inventory information, we selected regions with pure species and age composition. This allowed us to form the data set for training the classifier. Table 1 gives information about the ensembles obtained, which contain the measured spectral radiances for the indicated above homogeneous forest stands. The total number of spectra and the time of measurement are presented. We can see the best representation of the pine trees within these ensembles (14 gradations from the young forest to the age near to 140 years). Two examples of the same age (76 years) for the lines 8 and 9 of Table 1 should be perceived as the only plot of the pine species that was encompassed by two different tracks with the total numbers of the pixels equal to 1557 and 4019, respectively. The procedure of binning channels usually serves to enhance the signal-to-noise ratio of the imaging spectrometer. We constructed the confidence intervals for the intrinsic noise of CCD and the signal on the basis of the measurements obtained. The ability of the classification on the basis of spectral features is doubtful in the case if the variability of the noise is significant relative to the signal variability. The relative contribution of the noise for different spectral channels is represented by Fig. 4. As we can see from Fig. 1 data, the characteristic feature of the HSC instrument is that it has the spectral channels with the higher resolution in the short-wave region and with the lower resolution in the long-wave region. The noise component of the instrument is much higher for the short-wave lengths than for the long-wave lengths. We have revealed that some of the short-wave channels are pure atmospheric and carry no information about the underlying forests in our attempts to test the classifier, but the problem emerges to account for joining up the most informative channels to enhance their signal-to-noise ratio. Let us presume the HSC measurements are suitable for solving the considered problem if the noise contribution is about 5-10%. In Fig. 4(a) we can see that the first 184 channels up to the yellow region do not mount the 5% threshold. The noise is maximal for the channels of
#210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15416
the violet region where the noise contribution is about 70%. Channels of the near infrared region have the noise contribution about 1-2%.
Fig. 3. The test area investigated during the airborne campaign of hyperspectral measurements in Tver region (Russia). Locations of the tracks obtained are indicated by the red lines. Table 1. Parameters of the Spectral Radiance Ensembles for Homogeneous Forest Areas Species
Pine
Birch Aspen Elm
Age
Number of spectra
13 16 26 36 47 56 66 76 76 86 96 106 116 126 136 16 51 71 11 –
1046 2428 2061 447 5025 1807 7551 1557 4019 2055 3156 2191 644 1932 695 1729 1634 5656 2545 534
Local time of the airplane tracks 11-10-57 11-31-32 12-05-23 11-10-57 11-17-58 11-46-25 11-52-48 11-52-48 11-59-25 11-46-25 11-38-02 11-59-25 11-24-25 11-04-17 11-46-25 11-31-32 11-38-02 11-46-25 11-38-02 10-47-53
The procedure of binning improves the situation. The binning is performed by the method of running window which has some fixed spectral interval. The reduction of the spectral
#210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15417
resolution to 2 nm (Fig. 4(b)) leads to mounting of the 5% threshold by the channels of the green region. The significant part of the blue region channels can be used at the binning up to 5 nm (Fig. 4(c)). All joined channels mount 10% threshold in the case of binning up to 30 nm. However, it was found a useful information applicable for a fine classification (for instance, to retrieve the tree age) that can be lost under such spectral resolution. For that reason we have stopped on the binning up to 5 nm. After the binning procedure, we obtained available 87 aggregated spectral channels in the visible and near infrared ranges. Before solving the recognition problem, we found and eliminated the bands with significant radiometric distortion which are shaped by the stripes along the trajectory of the flight appearing due to the systematic calibration error. The most part of the stripes are not visible by eyes on the channel images. The distorted channels can be detected from the results of the binary classification if the relevant scene contains a large size object comparable with the width of the airplane track. Using this method, we have found 11 corrupted channels to be eliminated. The distortion in these channels can't be completely corrected by standard statistical methods. More exact calibration is needed for the correction.
Fig. 4. Contribution of the noise variability to the signal variability for the hyperspectral measurements used. a) – no binning; b) – binning to 2 nm; c) – binning to 5 nm; d) – binning to 30 nm. The horizontal dashed line signifies 5% level of the noise variability.
The more complicated than the classification of the objects according to their spectral features is the recognition of the forest classes of different species and age on the test area. Some results were shown in [23] that the accuracy of the recognition is enhanced if textures of the forest canopies are taken into account. The texture is formed by the end-members of the sunlit tops, shaded space and intermediate parts of tree’s crowns. If to separate the pixels characterizing the relevant forest classes, then the textures can be modeled and the accuracy improved. The necessity to normalize the registered spectral radiances was mentioned to make the recognition algorithms more universal in extrapolating the created spectral database on other processed scenes. The normalized on their integral values radiances for these three information layers were maintained to be represented by the 3D graphs within the central wavelength and age axes [17]. This form of normalization leads to a decrease of variations in the spectral curves for the
#210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15418
forest stands of all ages and to an inversion of the spectra in the short-wave channels without additional changes in the long-wave channels. The inversion is due to the normalization procedure where the amplitudes of the normalized radiances for the shaded background are increased in the short-wave channels as compared with tree’s sunlit tops. As a result, the short-wave layers of these graphs for different ages were shown to be sensitive to the shaded parts of the pine and birch canopies and the long-wave layers to their illuminated parts. Figure 5 demonstrates a distribution of the normalized radiances for the pixels corresponding to the sunlit tops versus the age of the forest trees for the samples given by Table 1 of the pine (a) and birch (b) species. All spectra are presented at the left parts of these graphs and only for the channels of the visible region at the right parts. These are crosssections of the mentioned 3D graphs for the normalized radiances.
Fig. 5. The along-channel distribution of the normalized spectral radiances for the sunlit endmembers relating to the pine (a) and birch (b) forest classes of different ages. The colors presented correspond to central wavelengths of HSC spectral channels. All spectra are on the left side and spectra in the visible region only are on the right side.
We can see that the radiances for the completely illuminated by the Sun pine tree forests (Fig. 5(a)) behave differently in the visible and near infrared regions. Oscillations are apparent in the visible region with a period near to 40 years and they are even more pronounced in the near infrared region. Not understandable is the alternation of the colors in the visible region: the amplitudes of the spectra are increased for the red, orange, yellow, dark-green and lightgreen curves, correspondingly, being the smallest for the curves represented in the blue color. The more monotonic is the behavior of the curves on Fig. 5(b) for the birch tree forests where no such oscillations can be remarked, but the color distribution is similar to that on Fig. 5(a). These characteristic features of the coniferous and deciduous species depending on their ages need to be thoroughly investigated with the forest researchers.
#210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15419
To make the pattern recognition algorithms more universal it is possible to use as characteristic features not only the normalized radiances, but also the relative levels of the registered radiances. The information layers were formed according to the ensembles of the spectra. The ensembles obtained model the projective characteristics and the density of the forest canopies. Both normalized radiances and their relative levels were stored by computer means as particular cells for the selected species for further recognition of different types of the forest classes. We have obtained that the best classification accuracy corresponds to meadow vegetation. In our data set we have spectra for 3 types of the distinguished grasses. The total probability of misclassification is about 1%. Water objects include central and coastal parts of Volga river, Orsha river, the sand pit and the pond on the selected test area. The misclassification errors are less than 3%. The test area contains 4 different types of open sandy soils distinguished by different admixtures (peat, gravel, clay, etc). For this subclass the error is about 3-4%. Different types of the asphalt and concrete roads are discriminated with errors 520%. The classification accuracy of the forest vegetation depends on the problem. In general, we can conclude that the main errors of the forest species recognition are related to the shaded background since the imaging spectrometer does not always ensure the high values of the signal-to-noise ratio for these scenes as was told above. We have proven for Table 1 data that the mature forest stand composed of pines and birches can be recognized with the errors less than 2% using this optimization of the spectral channels. For the young forest the error becomes greater (about 5%). We also used data for the aspen species which appears to distinguish from the pine better than from the birch of the same age.
Fig. 6. Confusion matrix (probabilities of true classification) of the broad age - class recognition for the pine and birch species with the fully illuminated parts of their crowns.
The situation is worse if to consider the forest stands of different ages. The confusion matrix was constructed for the available and predicted ages of the pine species. The related probabilities of misclassification for some ages are less or equal to 50%. Thus, the age of the forest stands can hardly be recognized by the selected stable procedures with low number of the spectral channels. Such number is insufficient, but we present Fig. 6 as an example of the broad age-class recognition: young, ripening and mature forest stands for all data of the birch species (Table 1) and for sampling data of the pine species for the ages 16, 76, 136. We can see that the most part of misclassification is observed for the birch stands (the probability error reaches 20%). The pine stands of the considered age-classes are discriminated with the much higher accuracy. Our opinion is that these differences in the accuracy of the stands #210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15420
recognition are due to the vertical structure of the birch and pine trees resulting in the effect of more sunlit tops appeared in the ensembles for the pine forest as compared with the birch forest. Some example of the forest species recognition is shown by Fig. 7. The test area corresponds to the airplane track 11-59-25 encompassing samples for the pine forest stands of ages 76 and 106 years as represented in Table 1. The size is 654 pixels along the track and 270 pixels across the track. The location on the forest inventory map of the Savvatyevskoe forestry is represented in Fig. 7(a) (marked by the yellow frame).
Fig. 7. Results of the forest tree species recognition for the selected 20 test regions. Location of the test area on the forest inventory map of the Savvatyevskoe forestry (a): orange – pine dominance, blue – birch dominance, magenta – spruce dominance; the RGB-synthesized image (b) of the area; the forest species composition of the area (c): 1 – for the pixels corresponding to the completely shaded background, 2 – for the pixels with half of the forested environment illuminated by the Sun and half of the shaded background, 3 – for the pixels referenced as the sunlit tops.
The typical forest inventory in Russia usually deals with the ground-based map schemes for a particular area as separated quarters and plots within them. The schemes are represented #210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15421
by specific colors (the plots for prevailing birch species by blue: light blue for young forest, dark blue for aged forest; the same for prevailing pine species by orange; the former logging places are denoted by the horizontal lines, etc.). Each quarter and plot has its own number on these maps. 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). Not all plots are presented by this full set of numbers. The part of the forest inventory map corresponding to the test area is zoomed in Fig. 7(a). The quarters are denoted by bold numbers. The RGB-synthesized image is presented on Fig. 7(b). Figure 7(c) depicts the results of the species recognition for the illustrated area in accordance with the three information layers of the forest characterization: the sunlit tops, the shaded background, the intermediate pixels of partly illuminated and partly shaded conditions. Such separation of the pixels may be seen to detail the texture of the image. Numbers on Fig. 7(c) for the relevant plots denote the percentage cover of the prevailing species for the ground-based mapping scheme: P – pine, S – spruce, B – birch. A rather good agreement can be marked between the recognition results and the forest inventory data. Though the signal-tonoise ratio is different for the pixels representing the illuminated and background parts of the related end-members, the ability of the classifier to describe the texture of the selected classes contributes to the total accuracy enhancement of the pattern recognition procedures. The classification error of the species composition for the selected 20 regions of Fig. 7(c) was calculated in accordance with the outlined 3 gradations of pixels by the formula
ε (i ) =
( p (i ) pine − pˆ (i ) pine ) 2 + ( p(i )birch − pˆ (i )birch ) 2 2
,
where i is the number of the corresponding region, p pine and pbirch are the percentage areas of the pine and birch, respectively, according to the forest inventory, pˆ pine and pˆ birch are the similar values given by the recognition. The errors may be large for some regions, but were found to be minimal for the extended regions. Figure 8 gives the analysis of the classification errors of the species composition for the regions. There is a natural level of the errors due to the boundary pixels (the curve of black color on Fig. 8). The errors might be significant for small plots due to their possible georeferencing problem. In general, we can see the contribution of the illuminated, shaded and partly illuminated and shaded pixels for the selected regions. The errors are higher for the smaller regions, but the weighted errors seem to be acceptable for the total test area and different gradations of the sun illuminated pixels. We have calculated their contribution by the formula 20
ε total = ε (i )W (i ), i =1
Area of region are the weighting functions of the plots. The total Total area of regions weighted error amounts to 12.2%, for shadowed pixels it is 13.4%, for partly illuminated – 11.1% and for completely illuminated – 12.3%. We should have in mind that the accuracy of ground-based data for the forest species composition is estimated as near to 10%.
where W (i ) =
#210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15422
Fig. 8. The color curves represent the regional classification errors of the species composition (measured in percents) for the pixels corresponding to the shaded, semi-sunlit and sunlit areas and for all pixels. The black curve corresponds to the natural uncertainty because of the pixels on the boundaries between different regions.
4. Conclusion We have addressed the problem of the forest species and age recognition using airborne hyperspectral imagery processing. The mathematical formalism of data processing includes the optimization procedures of the Bayesian classifier constraints for the pattern recognition by the spectral and texture features as well as accounting for the spectral channel sets to reduce their possible redundancy in the particular subject area. We have tested only the improved Bayesian classifier understanding that other machine-learning algorithms might also exist for the proposed purposes. This testing is due to our intention to realize the processing procedures using the parallel algorithms. Other known classifiers may not satisfy the demands of parallelization. We have obtained the higher accuracy of the pattern recognition by selecting pixels characterizing sunlit tree tops, sunlit background and shadows for particular classes of coniferous and deciduous forests. We consider that in the particular case study the vertical structure of the pine and birch trees results in appearance of larger statistical ensembles for the sunlit tops of the coniferous species of the same age as compared with the similar ensembles for the deciduous species. This leads to the higher accuracy in the recognition of the pine species of different age, though this problem needs additional consideration in the context of sampling materials obtained from the airborne imaging spectrometer and ground-based campaigns. The separation of the illuminated, shaded, and partly shaded tree’s crown regions serves to understand the accuracy of the pattern recognition, once the imaging spectrometer producers ensure high values of signal-to-noise ratio mainly for the sunlit tops, but not for the remaining background pixels. The relevant automation procedures of data processing are designed to improve the efficiency of the proposed techniques by high-productive computer systems. These techniques can be used by researchers and/or managers in the forestry discipline. Acknowledgments This work is supported by the Russian Fund for Basic Research (# 13-01-00185, 14-0500598, 14-07-00141).
#210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15423
1 M.V. Lomonosov Moscow State University, Leninskiye Gory, GSP-1, Moscow, 119991, Russia Moscow Institute of Physics and Technology (State University), Institutskiy per., 9, Dolgoprudny, Moscow Region, 141700, Russia 3 Institute of Numerical Mathematics of Russian Academy of Science, ul. Gubkina, 8, Moscow, Moscow, 119333, Russia 4 Laboratoire de PhysicoChimie de l'Atmosphère, Université du Littoral Côte d’Opale & Université Lille Nord de France, Avenue Maurice Schumann, 189A, Dunkerque, 59140, France * [email protected]
Abstract: Optical remote sensing data processing is proposed for the airborne images of high spectral and spatial resolution. Optimization techniques are undertaken to gain information about spatial distribution of the pixels on the hyperspectral images and the texture of the forest stands of different species and ages together with reducing redundancy of the spectral channels used. The category of neighborhood of pixels for particular forest classes and the step up method of selecting optimal spectral channels are employed in the relevant processing procedures. We present examples of pattern recognition for the forests as a result of separating pixels, which characterize the sunlit tops, shaded space and intermediate cases of the Sun illumination conditions on the hyperspectral images. ©2014 Optical Society of America OCIS codes: (280.1415) Biological sensing and sensors; (280.4788) Optical sensing and sensors; (280.4991) Passive remote sensing.
References and links 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16.
S. Z. Li, Markov Random Field Modeling in Computer Vision (Springer-Verlag, 1995). V. V. Kozoderov and E. V. Dmitriev, “Remote sensing of soils and vegetation: pattern recognition and forest stand structure assessment,” Int. J. Remote Sens. 32(3), 5699–5717 (2011). A. K. Jain, “Advances in mathematical models in image processing,” Proc. IEEE 69(5), 502–528 (1981). N. S. Friedland and A. Rosenfeld, “Compact object recognition using energy-based optimization,” IEEE Trans. Pattern Anal. Mach. Intell. 14(7), 770–777 (1992). B. Tso and R. C. Olsen, “A contextual classification scheme based on MRF model with improved parameter estimation and multiscale fuzzy line process,” Remote Sens. Environ. 97(1), 127–136 (2005). J. Bolton and P. Gader, “Random set framework for context-based classification with hyperspectral imagery,” IEEE Trans. Geosci. Rem. Sens. 47(11), 3810–3821 (2009). V. V. Kozoderov, “Application of optical remote sensing data to study natural and climatic processes,” Climate and Nature 2(3), 3–16 (2012). F. G. Hall, D. E. Knapp, and K. F. Huemmrich, “Physically based classification and satellite mapping of biophysical characteristic in the southern boreal forest,” J. Geophys. Res. 102(D24), 29567–29580 (1997). T. Hakala, J. Suomalainen, S. Kaasalainen, and Y. Chen, “Full waveform hyperspectral LiDAR for terrestrial laser scanning,” Opt. Express 20(7), 7119–7127 (2012). E. V. Dmitriev and V. V. Kozoderov, “Optimization of spectral bands for hyperspectral remote sensing of forest vegetation,” Proc. SPIE 8887, 888705 (2013). J. Suomalainen, T. Hakala, H. Kaartinen, E. Räikkönen, and S. Kaasalainen, “Demonstration of a virtual active hyperspectral LiDAR in automated point cloud classification,” ISPRS J. Photogramm. Remote Sens. 66(5), 637– 641 (2011). M. Dalponte, L. Bruzzone, and D. Gianelle, “Fusion of hyperspectral and LIDAR remote sensing data for classification of complex forest areas,” IEEE Trans. Geosci. Rem. Sens. 113, 1416–1427 (2009). V. V. Kozoderov and E. V. Dmitriev, “Remote sensing of soils and vegetation: regional aspects,” Int. J. Remote Sens. 29(9), 2733–2748 (2008). V. V. Kozoderov, T. V. Kondranin, E. V. Dmitriev, and V. P. Kamentsev, “Mapping forest and peat fires using hyperspectral airborne remote-sensing data,” Izv., Atmos. Ocean. Phys. 48(9), 941–948 (2012). J. Besag and P. J. Green, “Spatial statistics and Bayesian computation,” J. R. Stat. Soc., B 55(1), 25–37 (1993). M. Bertero, T. A. Poggio, and V. Torre, “Ill-posed problems in early vision,” Proc. IEEE 76(8), 869–889 (1988).
#210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15410
17. V. V. Kozoderov, E. V. Dmitriev, and V. P. Kamentsev, “A system of airborne remote sensing data processing of high spectral and spatial resolution,” Issledovanie Zemli Iz Kosmosa. 6, 57–64 (2013). 18. V. V. Kozoderov, T. V. Kondranin, and E. V. Dmitriev, Thematic processing of multispectral and hyperspectral airspace images (MIPT Publication, 2013). 19. S. Z. Li, “Object recognition from range data prior to segmentation,” Image Vis. Comput. 10(8), 566–576 (1992). 20. R. Kohavi, “A study of cross-validation and bootstrap for accuracy estimation and model selection,” Int. Joint Conf. on Artificial Intell. (IJCAI), 528–535 (1995). 21. T. Hastie, R. Tibshirani, and J. Friedman, The Elements of Statistical Learning (Springer, 2001). 22. M. Dalponte, L. Bruzzone, L. Vescovo, and D. Gianelle, “The role of spatial resolution and classifier complexity in the analysis of hyperspectral images of forest areas,” Remote Sens. Environ. 113(11), 2345–2355 (2009). 23. V. V. Kozoderov, T. V. Kondranin, and E. V. Dmitriev, “ An apparatus and programmatic system of hyperspectral airspace imagery processing,” Proc. Int. Symp. Atm. Radiation and Dynamics (ISARD), 41 (2013).
1. Introduction The vision theory of the optical airborne data processing implies an optimal computational procedure realization to retrieve parameters from registered radiances for the land surface objects called “patterns” [1]. Pattern recognition techniques are a valuable approach to the processing of remote sensing data for images of the land surface – atmosphere system [2]. The related cognitive technologies originate from the mathematical models of artificial intelligence and data mining [3]. This discipline has passed through the computational procedures of the objects classification on the processed images [4] to the updated techniques of the contextual recognition of the objects as fuzzy sets [5] and the random set framework [6]. The contextual recognition for the forested environments assumes that the neighboring pixels are taken into account for the particular classes of forests of different species and ages [7] since these classes are formed by the “end-members” (sunlit forest canopy, sunlit background and shaded background for a particular view and solar illumination angle) [8]. Each pixel views collections of the individual forest trees and the pixel-level reflectance can thus be computed as a linear mixture of sunlit tree tops, sunlit background (or understory) and shadows. Applications of optical remote sensing include the use of airborne photo-pictures, hyperspectral images and the laser scanning (lidar, light detection and ranging) data. The hyperspectral imagery processing enables to gain information about the distribution of spectral and spatial characteristics of the forest stands of different species composition and age classes. Data from lidars and imaging spectrometers with Charge Coupled Devices (CCD) serve to combine the priorities of the active and passive remote sensing instruments to carry out both these types of measurements from the same airborne platform [9]. The idea can be realized by using high-productive computers and optimization techniques of a neighborhood measure extraction for selected forest classes and of accounting for the optimal set of the spectral channels, thus reducing the possible redundancy of the entire hyperspectral data sets [10]. The combined activity contributes to obtaining the 3D products for the forest canopy [11, 12]. The related mathematical treatments of using different types of classifiers combined with the parameter retrieval procedures for coniferous and deciduous forest classes are maintained to have much wider applications as compared with the image features and object shapes extraction, which relates to photometry and geometry in pixel-level reflectance representation of the forested land cover [13]. The pixel fraction and reflectance of the end-members are only a part in the listed techniques. Instead of the photometry and geometry constraints, the improved models are developed of the functional description of outgoing spectral radiation, in which such parameters of the forest canopy like the vegetation biomass density for particular forest species and ages are embedded [14]. This permits us to calculate the relationships between the registered radiances and the forest biomass densities (the direct problem of atmospheric optics). The next stage is to find solutions of this problem as cross-sections of the related curves in the multi-dimensional space given by the parameters of these models (the inverse problem).
#210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15411
The typical solutions of the inverse problem may not be mathematically unique and the computational procedure is undertaken to their regularization by finding minima of the functional called “the energy for the particular class of forests”. In general, the energy category is considered as a linking chain between the image pixels as random fields and our understanding what class of objects each pixel belongs to, in the known Bayesian paradigm of statistical decision making [15]. The optimization 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 [16]. This article presents some results of the forest species and age recognition for a test area using the listed mathematical framework for ensembles of airborne data from an imaging spectrometer produced in Russia operating in near to 200 spectral channels from the wavelengths 400 nm to 1000 nm. Our main purpose was to arrange a mutual airborne and ground-based field campaigns on a test area with geo-referencing pixels of the hyperspectral images to the standard description of the relevant “pure stands” with the ultimate goal to answer the question, whether we can have a practical tool of the pattern recognition of the forest species and age classes or not. Our prior assumption was as follows. Should we recognize the texture of the young, mature and old forest that is perceived as an alternation of the pixels characterizing the sunlit tops and the shaded background by the computer vision techniques, and we can find such a tool within perceptual grouping of the pixels using their spectral features. As usual, a confusion matrix will enable to compare the predicted ages and the ages given by a particular map of the forest inventory for the test area. 2. Instruments and processing methods The measuring system used for obtaining simultaneous photo and hyperspectral high resolution images of the test area in the nadir direction consists of the hyperspectral camera (HSC), the power module, the notebook with the available software, the gyro-stabilized platform and the professional photo camera. The light-weight airborne prismatic HSC is produced and calibrated by Science and Production Enterprise “Lepton” (Zelenograd, Moscow, Russia). The complex direct-vision prism consisting of 7 optically glued simple prisms supplies the linearity of the dispersion curve in the visible and near infrared (VNIR) region 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 basic parameters of HSC are as follows (Fig. 1). The camera has 287 active and 3 noise channels in the spectral range 401-1017 nm. The spectral resolution of the channels significantly depends on the central wavelengths and changes from 0.36 nm to 14 nm. The radiometric resolution of the digitized signals is 212 − 31 = 4065 levels. The central wave lengths are presented at Fig. 1(a). The main channels belong to the visible region. The typical feature of such prism instruments is in their different spectral resolution for the short-wave and long-wave regions. In particular, the width of the channels increases more than on order in the near infrared region (Fig. 1(b)). The measuring system passed an obligatory preflight check and routine maintenance between the flights. The instrument is calibrated by the Enterprise “Lepton” every year. Before the flight, this imaging spectrometer takes special tests to study the linearity of the present calibration and the relationship between signal and noise. An optimal storage time, an effective number of gray-scale levels and spectral channels are maintained during each flight.
#210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15412
Fig. 1. The central wavelengths (a) and the spectral resolution (b) of the imaging spectrometer.
The algorithms used for the retrieval of forest stand attributes are based on the Bayesian classifier, which is optimal in the sense of minimization of the misclassification probability. According to the Bayesian theory of decision making, when both the prior distribution and the likelihood function of a pattern are known, the best that can be estimated from these sources of knowledge is given by this classifier. The learning process of the Bayesian classification algorithm consists in the estimation of prior probabilities of the classes and corresponding likelihood functions. The prior probabilities are obtained from results of texture analysis in terms of Markov random fields (MRF). We follow here the account of this problem given in [1] concerning the labeling problem in computer vision with the relevant contextual constraints and the neighborhood system for pixels of a particular class of objects. MRF theory provides a consistent way of modeling context depending entities such as image pixels. This is achieved through characterizing mutual influences among such entities using MRF probabilities. The local properties of MRFs lead to parallel computational algorithms. The conditional probability density functions of features are modeled in terms of Gaussian mixtures. The known techniques are considered to be hard to implement for hyperspectral data processing because the number of learning samples necessary for providing the stable probability density function estimate increases too fast with the number of features. In practice, even if we have a set of learning samples large for all considered classes, the learning process would be very expensive in computational terms. The maximum a posteriori (MAP) solution is typical in this framework. The MAP-MRF labeling formalism assumes that a search of minima of a specified function called the energy of the image scene or of the particular class of objects is undertaken [17]. The Bayesian classification principle is applicable for different cases of optical image processing [18]. The classification schemes are based on the assumption that similar ground objects have similar spectral and texture characteristics. The shape of registered spectra depends on many factors and may have significant variations even for objects which seem to be homogeneous. The reflected radiation spectra can significantly change according to the volume distribution of leaves and branches, the state of leaves and underlying surface even in the case when such basic attributes of forest stands as species composition, ages and heights are the same. Modeling of spectral reflectivity of vegetation and other ground objects demands the availability of many parameters which cannot be easily obtained in practice. Applicability of the Bayesian classifiers for processing of airborne images of high spatial resolution (1-2 m) is restricted by the fact that the imaging spectrometer producers can ensure very high signal-to-noise ratio for the pixels relating to the completely illuminated by the Sun tree’s crowns. Meanwhile, the ratio for the background pixels is not so high. The MAP-MRF labeling concerns a supposition of the pixels distribution function. If we have an image or its separate segments which can be described by a constant radiance function for all pixels, then we deal with the known in mathematics string problem (the membrane in the two-dimensional #210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15413
case) [19]. The energy for this class of objects with very low variety of radiances (or the biological diversity of vegetation) has its zero minimum and is positive otherwise. If the prior supposition is in constant gradient of the registered radiances for the land surface, then we deal the rod problem, which energy is represented by the integral from the second derivative function to the second power. It is possible to make a prior supposition about the constant curvature for the class of objects under consideration and to obtain the formula for the energy in the form of the integral from the third order derivative function to the second power, etc. The posteriori energy minimization for any class of objects results in smoothing down the differences between the radiances for the listed three categories of the related pixels (sunlit tops, backgrounds and the intermediate cases of partly illuminated and partly shaded endmembers). Additionally to these smoothing constraints, we have elaborated a rigorous approach to optimize spectral channels of HSC based on searching their most informative sets. The high dimensionality of the feature space at the fixed number of samples of the learning set is typical for hyperspectral data processing. It makes the learning process unstable because of the known problem of the “curse of dimensionality”. The parametric Bayesian classifiers also suffer from this problem to a lesser extent than the nonparametric ones. One of the frequently used techniques for the feature space optimization is the stepwise forward selection method. For the classification algorithm we can apply this method as follows. At each step the features are subdivided into two groups: features already included in the classifier and other features which could be potentially included. After that, features from the second group are added one-by-one to the classifier and corresponding classification errors are estimated. The feature corresponding to the maximum decrease of error is selected. If the error decrease is statistically significant, the feature is included in the classifier and the process continues. Otherwise, the process is stopped. The standard form of the stepwise forward selection has the following two known disadvantages. At first, the classification errors are estimated in the dependent manner because the test set of samples is used for learning the classifier. Secondly, the optimal sequence of the features obtained by this method is highly sensitive to small variations of the learning set. Especially, it concerns the last members of the sequence. We propose to use a modified method allowing us to avoid these problems. The scheme of the method is demonstrated by Fig. 2(a). The stages presented at Fig. 2(a) serve to find the most probable sequence in the forward selection process by the series of launches to estimate the classification error and finding the most frequent features on the different levels of their extraction from the ensembles of the selected features. This calculation process is necessary to select the most informative features, which lead to the minimal classification error. The number of the learning samples for each of the classes is often large for the estimation of the classification error using the holdout cross-validation method [20, 21]. In this method, the initial set of samples is randomly, with the given probability divided on two subsets – training and testing. We used the probability 0.5 and thus the subsets obtained are almost equal. The classifier is always learned using the training subset and the classification errors are estimated in independent manner only for the testing subset. Thus, employing this method for the stepwise forward selection, we avoid the overtraining problem that is more important than the decrease of accuracy in the related statistical estimates due to the reduction of the number of training samples. If the number of the learning samples is not sufficiently large, it is possible to use the similar algorithm based on the leave-one-out cross-validation [22], however this method is much more expensive in computational terms. Random partitioning of the learning data set using the holdout cross-validation is repeated many times. Then the stepwise forward selection is launched for each pair of the training and testing sets obtained. As a result, we have a locally optimal sequence of features. The sequences obtained as a rule are different in composition and length. The unique regularized solution of the feature selection problem can be obtained as the most probable sequence.
#210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15414
Fig. 2. Optimization of the feature space: a) – scheme of the selection of the subset of the most informative features; b) – an example of the selection of the most probable sequence of features. Data corresponding to different levels are highlighted by the grey scale.
In the standard technique, the leading feature corresponds to the minimum total probability of misclassification. However, it must be taken into account that the dissimilarity of the classified objects can be very different. In this case, the choice of the leading feature is mostly defined by the presence of similar objects corresponding to high misclassification errors. On the other hand, some well-distinguished objects, for instance, water, soil and vegetation areas can be classified exactly using just one spectral channel. The channels allowing the exact classification can be not unique and change for different groups. Thus, in our study we define the leading feature as the spectral channel allowing the exact classification of the maximum number of different objects available in the learning set. The algorithm of finding the most probable sequence of features is demonstrated by Fig. 2(b). At the beginning (level 0), we form the set of the members of the sequences coming
#210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15415
after the leading feature and calculate the mode. Then, we form a subset of the sequences corresponding to the mode (level 1) and repeat the process with the next members. The process continues (levels 2, 3, 4) up to selecting a unique sequence. In the case of a few equivalent sequences, the random choice is used. The regularized solution is found to be much less sensitive to variations of the learning samples. The robustness of this method increases with the number of the launches. As a result of the listed optimization constraints, we can select the information layers for a particular class of forests, which pixels belong to the sunlit tops, shaded space and intermediate cases of the Sun illumination conditions on the hyperspectral images. Depending on the species and age of the forests, these pixels for the same class have different spectral distributions and the information layers serve to improve the recognition accuracy. 3. Results and discussion The described techniques of hyperspectral imagery processing were tested on data of the airborne campaign conducted in Tver region of Russia. The test region of the size about 4x10 km was encompassed by 13 overlapped images obtained from the airplane equipped on the same gyro-stabilized platform by the imaging spectrometer and photo-camera. The location of the region and corresponding tracks of the flights are represented in Fig. 3. Each hyperspectral image contains 500 pixels across the track and approximately from 10000 to 14000 pixels along the track. The spatial resolution across the track is stable and amounts to 1.1 m at the flight altitude about 2 km above the ground level. The pixel size along the track depends on the flight speed and changes within 0.66-0.91 m. The landscape of the test region contains different types of natural and artificial objects: the most easily distinguishable objects in accordance with their spectral and texture features are water bodies, open soils, asphalt roads, buildings, forests and grasslands. As we can see from Fig. 3, about half of the area is covered by the forest stands for which the ground-based forest inventory maps are available. Besides that, an extensive field campaign was conducted on the test area to arrange its geo-botanical, forest typology and other descriptions for georeferencing the selected “pure stands” to make the validation and ground truth estimates. Analyzing the forest inventory information, we selected regions with pure species and age composition. This allowed us to form the data set for training the classifier. Table 1 gives information about the ensembles obtained, which contain the measured spectral radiances for the indicated above homogeneous forest stands. The total number of spectra and the time of measurement are presented. We can see the best representation of the pine trees within these ensembles (14 gradations from the young forest to the age near to 140 years). Two examples of the same age (76 years) for the lines 8 and 9 of Table 1 should be perceived as the only plot of the pine species that was encompassed by two different tracks with the total numbers of the pixels equal to 1557 and 4019, respectively. The procedure of binning channels usually serves to enhance the signal-to-noise ratio of the imaging spectrometer. We constructed the confidence intervals for the intrinsic noise of CCD and the signal on the basis of the measurements obtained. The ability of the classification on the basis of spectral features is doubtful in the case if the variability of the noise is significant relative to the signal variability. The relative contribution of the noise for different spectral channels is represented by Fig. 4. As we can see from Fig. 1 data, the characteristic feature of the HSC instrument is that it has the spectral channels with the higher resolution in the short-wave region and with the lower resolution in the long-wave region. The noise component of the instrument is much higher for the short-wave lengths than for the long-wave lengths. We have revealed that some of the short-wave channels are pure atmospheric and carry no information about the underlying forests in our attempts to test the classifier, but the problem emerges to account for joining up the most informative channels to enhance their signal-to-noise ratio. Let us presume the HSC measurements are suitable for solving the considered problem if the noise contribution is about 5-10%. In Fig. 4(a) we can see that the first 184 channels up to the yellow region do not mount the 5% threshold. The noise is maximal for the channels of
#210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15416
the violet region where the noise contribution is about 70%. Channels of the near infrared region have the noise contribution about 1-2%.
Fig. 3. The test area investigated during the airborne campaign of hyperspectral measurements in Tver region (Russia). Locations of the tracks obtained are indicated by the red lines. Table 1. Parameters of the Spectral Radiance Ensembles for Homogeneous Forest Areas Species
Pine
Birch Aspen Elm
Age
Number of spectra
13 16 26 36 47 56 66 76 76 86 96 106 116 126 136 16 51 71 11 –
1046 2428 2061 447 5025 1807 7551 1557 4019 2055 3156 2191 644 1932 695 1729 1634 5656 2545 534
Local time of the airplane tracks 11-10-57 11-31-32 12-05-23 11-10-57 11-17-58 11-46-25 11-52-48 11-52-48 11-59-25 11-46-25 11-38-02 11-59-25 11-24-25 11-04-17 11-46-25 11-31-32 11-38-02 11-46-25 11-38-02 10-47-53
The procedure of binning improves the situation. The binning is performed by the method of running window which has some fixed spectral interval. The reduction of the spectral
#210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15417
resolution to 2 nm (Fig. 4(b)) leads to mounting of the 5% threshold by the channels of the green region. The significant part of the blue region channels can be used at the binning up to 5 nm (Fig. 4(c)). All joined channels mount 10% threshold in the case of binning up to 30 nm. However, it was found a useful information applicable for a fine classification (for instance, to retrieve the tree age) that can be lost under such spectral resolution. For that reason we have stopped on the binning up to 5 nm. After the binning procedure, we obtained available 87 aggregated spectral channels in the visible and near infrared ranges. Before solving the recognition problem, we found and eliminated the bands with significant radiometric distortion which are shaped by the stripes along the trajectory of the flight appearing due to the systematic calibration error. The most part of the stripes are not visible by eyes on the channel images. The distorted channels can be detected from the results of the binary classification if the relevant scene contains a large size object comparable with the width of the airplane track. Using this method, we have found 11 corrupted channels to be eliminated. The distortion in these channels can't be completely corrected by standard statistical methods. More exact calibration is needed for the correction.
Fig. 4. Contribution of the noise variability to the signal variability for the hyperspectral measurements used. a) – no binning; b) – binning to 2 nm; c) – binning to 5 nm; d) – binning to 30 nm. The horizontal dashed line signifies 5% level of the noise variability.
The more complicated than the classification of the objects according to their spectral features is the recognition of the forest classes of different species and age on the test area. Some results were shown in [23] that the accuracy of the recognition is enhanced if textures of the forest canopies are taken into account. The texture is formed by the end-members of the sunlit tops, shaded space and intermediate parts of tree’s crowns. If to separate the pixels characterizing the relevant forest classes, then the textures can be modeled and the accuracy improved. The necessity to normalize the registered spectral radiances was mentioned to make the recognition algorithms more universal in extrapolating the created spectral database on other processed scenes. The normalized on their integral values radiances for these three information layers were maintained to be represented by the 3D graphs within the central wavelength and age axes [17]. This form of normalization leads to a decrease of variations in the spectral curves for the
#210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15418
forest stands of all ages and to an inversion of the spectra in the short-wave channels without additional changes in the long-wave channels. The inversion is due to the normalization procedure where the amplitudes of the normalized radiances for the shaded background are increased in the short-wave channels as compared with tree’s sunlit tops. As a result, the short-wave layers of these graphs for different ages were shown to be sensitive to the shaded parts of the pine and birch canopies and the long-wave layers to their illuminated parts. Figure 5 demonstrates a distribution of the normalized radiances for the pixels corresponding to the sunlit tops versus the age of the forest trees for the samples given by Table 1 of the pine (a) and birch (b) species. All spectra are presented at the left parts of these graphs and only for the channels of the visible region at the right parts. These are crosssections of the mentioned 3D graphs for the normalized radiances.
Fig. 5. The along-channel distribution of the normalized spectral radiances for the sunlit endmembers relating to the pine (a) and birch (b) forest classes of different ages. The colors presented correspond to central wavelengths of HSC spectral channels. All spectra are on the left side and spectra in the visible region only are on the right side.
We can see that the radiances for the completely illuminated by the Sun pine tree forests (Fig. 5(a)) behave differently in the visible and near infrared regions. Oscillations are apparent in the visible region with a period near to 40 years and they are even more pronounced in the near infrared region. Not understandable is the alternation of the colors in the visible region: the amplitudes of the spectra are increased for the red, orange, yellow, dark-green and lightgreen curves, correspondingly, being the smallest for the curves represented in the blue color. The more monotonic is the behavior of the curves on Fig. 5(b) for the birch tree forests where no such oscillations can be remarked, but the color distribution is similar to that on Fig. 5(a). These characteristic features of the coniferous and deciduous species depending on their ages need to be thoroughly investigated with the forest researchers.
#210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15419
To make the pattern recognition algorithms more universal it is possible to use as characteristic features not only the normalized radiances, but also the relative levels of the registered radiances. The information layers were formed according to the ensembles of the spectra. The ensembles obtained model the projective characteristics and the density of the forest canopies. Both normalized radiances and their relative levels were stored by computer means as particular cells for the selected species for further recognition of different types of the forest classes. We have obtained that the best classification accuracy corresponds to meadow vegetation. In our data set we have spectra for 3 types of the distinguished grasses. The total probability of misclassification is about 1%. Water objects include central and coastal parts of Volga river, Orsha river, the sand pit and the pond on the selected test area. The misclassification errors are less than 3%. The test area contains 4 different types of open sandy soils distinguished by different admixtures (peat, gravel, clay, etc). For this subclass the error is about 3-4%. Different types of the asphalt and concrete roads are discriminated with errors 520%. The classification accuracy of the forest vegetation depends on the problem. In general, we can conclude that the main errors of the forest species recognition are related to the shaded background since the imaging spectrometer does not always ensure the high values of the signal-to-noise ratio for these scenes as was told above. We have proven for Table 1 data that the mature forest stand composed of pines and birches can be recognized with the errors less than 2% using this optimization of the spectral channels. For the young forest the error becomes greater (about 5%). We also used data for the aspen species which appears to distinguish from the pine better than from the birch of the same age.
Fig. 6. Confusion matrix (probabilities of true classification) of the broad age - class recognition for the pine and birch species with the fully illuminated parts of their crowns.
The situation is worse if to consider the forest stands of different ages. The confusion matrix was constructed for the available and predicted ages of the pine species. The related probabilities of misclassification for some ages are less or equal to 50%. Thus, the age of the forest stands can hardly be recognized by the selected stable procedures with low number of the spectral channels. Such number is insufficient, but we present Fig. 6 as an example of the broad age-class recognition: young, ripening and mature forest stands for all data of the birch species (Table 1) and for sampling data of the pine species for the ages 16, 76, 136. We can see that the most part of misclassification is observed for the birch stands (the probability error reaches 20%). The pine stands of the considered age-classes are discriminated with the much higher accuracy. Our opinion is that these differences in the accuracy of the stands #210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15420
recognition are due to the vertical structure of the birch and pine trees resulting in the effect of more sunlit tops appeared in the ensembles for the pine forest as compared with the birch forest. Some example of the forest species recognition is shown by Fig. 7. The test area corresponds to the airplane track 11-59-25 encompassing samples for the pine forest stands of ages 76 and 106 years as represented in Table 1. The size is 654 pixels along the track and 270 pixels across the track. The location on the forest inventory map of the Savvatyevskoe forestry is represented in Fig. 7(a) (marked by the yellow frame).
Fig. 7. Results of the forest tree species recognition for the selected 20 test regions. Location of the test area on the forest inventory map of the Savvatyevskoe forestry (a): orange – pine dominance, blue – birch dominance, magenta – spruce dominance; the RGB-synthesized image (b) of the area; the forest species composition of the area (c): 1 – for the pixels corresponding to the completely shaded background, 2 – for the pixels with half of the forested environment illuminated by the Sun and half of the shaded background, 3 – for the pixels referenced as the sunlit tops.
The typical forest inventory in Russia usually deals with the ground-based map schemes for a particular area as separated quarters and plots within them. The schemes are represented #210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15421
by specific colors (the plots for prevailing birch species by blue: light blue for young forest, dark blue for aged forest; the same for prevailing pine species by orange; the former logging places are denoted by the horizontal lines, etc.). Each quarter and plot has its own number on these maps. 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). Not all plots are presented by this full set of numbers. The part of the forest inventory map corresponding to the test area is zoomed in Fig. 7(a). The quarters are denoted by bold numbers. The RGB-synthesized image is presented on Fig. 7(b). Figure 7(c) depicts the results of the species recognition for the illustrated area in accordance with the three information layers of the forest characterization: the sunlit tops, the shaded background, the intermediate pixels of partly illuminated and partly shaded conditions. Such separation of the pixels may be seen to detail the texture of the image. Numbers on Fig. 7(c) for the relevant plots denote the percentage cover of the prevailing species for the ground-based mapping scheme: P – pine, S – spruce, B – birch. A rather good agreement can be marked between the recognition results and the forest inventory data. Though the signal-tonoise ratio is different for the pixels representing the illuminated and background parts of the related end-members, the ability of the classifier to describe the texture of the selected classes contributes to the total accuracy enhancement of the pattern recognition procedures. The classification error of the species composition for the selected 20 regions of Fig. 7(c) was calculated in accordance with the outlined 3 gradations of pixels by the formula
ε (i ) =
( p (i ) pine − pˆ (i ) pine ) 2 + ( p(i )birch − pˆ (i )birch ) 2 2
,
where i is the number of the corresponding region, p pine and pbirch are the percentage areas of the pine and birch, respectively, according to the forest inventory, pˆ pine and pˆ birch are the similar values given by the recognition. The errors may be large for some regions, but were found to be minimal for the extended regions. Figure 8 gives the analysis of the classification errors of the species composition for the regions. There is a natural level of the errors due to the boundary pixels (the curve of black color on Fig. 8). The errors might be significant for small plots due to their possible georeferencing problem. In general, we can see the contribution of the illuminated, shaded and partly illuminated and shaded pixels for the selected regions. The errors are higher for the smaller regions, but the weighted errors seem to be acceptable for the total test area and different gradations of the sun illuminated pixels. We have calculated their contribution by the formula 20
ε total = ε (i )W (i ), i =1
Area of region are the weighting functions of the plots. The total Total area of regions weighted error amounts to 12.2%, for shadowed pixels it is 13.4%, for partly illuminated – 11.1% and for completely illuminated – 12.3%. We should have in mind that the accuracy of ground-based data for the forest species composition is estimated as near to 10%.
where W (i ) =
#210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15422
Fig. 8. The color curves represent the regional classification errors of the species composition (measured in percents) for the pixels corresponding to the shaded, semi-sunlit and sunlit areas and for all pixels. The black curve corresponds to the natural uncertainty because of the pixels on the boundaries between different regions.
4. Conclusion We have addressed the problem of the forest species and age recognition using airborne hyperspectral imagery processing. The mathematical formalism of data processing includes the optimization procedures of the Bayesian classifier constraints for the pattern recognition by the spectral and texture features as well as accounting for the spectral channel sets to reduce their possible redundancy in the particular subject area. We have tested only the improved Bayesian classifier understanding that other machine-learning algorithms might also exist for the proposed purposes. This testing is due to our intention to realize the processing procedures using the parallel algorithms. Other known classifiers may not satisfy the demands of parallelization. We have obtained the higher accuracy of the pattern recognition by selecting pixels characterizing sunlit tree tops, sunlit background and shadows for particular classes of coniferous and deciduous forests. We consider that in the particular case study the vertical structure of the pine and birch trees results in appearance of larger statistical ensembles for the sunlit tops of the coniferous species of the same age as compared with the similar ensembles for the deciduous species. This leads to the higher accuracy in the recognition of the pine species of different age, though this problem needs additional consideration in the context of sampling materials obtained from the airborne imaging spectrometer and ground-based campaigns. The separation of the illuminated, shaded, and partly shaded tree’s crown regions serves to understand the accuracy of the pattern recognition, once the imaging spectrometer producers ensure high values of signal-to-noise ratio mainly for the sunlit tops, but not for the remaining background pixels. The relevant automation procedures of data processing are designed to improve the efficiency of the proposed techniques by high-productive computer systems. These techniques can be used by researchers and/or managers in the forestry discipline. Acknowledgments This work is supported by the Russian Fund for Basic Research (# 13-01-00185, 14-0500598, 14-07-00141).
#210926 - $15.00 USD (C) 2014 OSA
Received 25 Apr 2014; revised 9 Jun 2014; accepted 10 Jun 2014; published 17 Jun 2014 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015410 | OPTICS EXPRESS 15423
