Introduction

In recent years, machine learning methods have been widely used for various tasks of automation and increasing the information content of procedures for thematic processing and analysis of aerospace images in the visible and near infrared spectral ranges. Multispectral satellite images of low and medium spatial resolution are traditionally used for survey of the soil and vegetation cover and the construction of large-scale thematic maps [1]. With the increase in the spatial and spectral resolution of satellite equipment, a number of novel tasks associated with remote sensing of natural and anthropogenic objects have arisen.

High (1-4 m) and very high (< 1 m) spatial resolution of panchromatic satellite images, forms the basis of methods for solving new tasks of monitoring land, forest and water resources, searching for mineral deposits and assessing ecological situation, which are more complex from the point of view of increased consumer requirements. There is a need to develop special approaches for analyzing large amounts of information and obtaining remote estimates of the characteristics of the examined objects with a given accuracy. Improving the efficiency of thematic processing of aerospace images of high spatial and spectral resolution is in high demand in many applications in the fields of natural resource management, agriculture, forestry and environmental monitoring [2].

The current trend in the development of methods for thematic processing of high resolution images is the combined use of spectral and texture features. So for example, a method of spectral-texture processing of aerial hyperspectral images of a forest canopy was presented in [3]. Analysis of the results of test calculations for selected areas of the Savvatyevskoe forestry (Russia, Tver) showed that the proposed approach provides a significant increase in the accuracy of classification of the species composition and age groups, in comparison with the averaged spectral characteristics. It should be noted that for synthesized multispectral images, taking into account spectral features showed an increase in accuracy by more than 10%.

The presented results on improving the accuracy due to the use of texture features are also confirmed by the comparison with the previously obtained results of thematic processing of hyperspectral images of nearby territories presented in [4]. Both effective nonparametric methods of cluster analysis [5] and optimized ensemble machine learning algorithms [6] can be successfully used for spectral-texture classifications.

The work [7] shows new possibilities of using statistical texture analysis of satellite images of very high spatial resolution to retrieve the structural parameters of forest stands, characterizing the variety of sizes and density of crowns, as well as the relative position of individual trees. The presented technique is based on the parameterization of linear relationships between the Haralick texture features and the structural parameters of pine stands. The results obtained can be effectively used to provide more accurate estimates of the aboveground biomass of forest stand fractions.

The accuracy of texture analysis depends on the chosen feature extraction method [8]. In this paper, we discuss the possibilities of using various statistical methods for measuring textures based on the matrix representation.

Texture feature extraction and classification methods

The panchromatic satellite images presented in grayscale are considered as an object of texture analysis since they have the highest spatial resolution. The texture is formed by the spatial arrangement and the mutual combination of structural elements. Natural objects are characterized by a random arrangement of structural elements and significant variations in their parameters, for example, tone and size. Thus, the task of constructing parameters characterizing a particular texture is associated with obtaining statistical estimates.

Statistical methods of texture analysis are based on assessing the spatial distribution of local characteristics of structural elements for all possible locations in the image and extracting statistical parameters from the obtained distributions of local characteristics. Matrix methods assume that the desired distribution is discrete and has a finite number of elements. An image for which a texture extraction is made must contain a sufficiently large number of structural elements to obtain reliable estimate of probability mass function. Texture features obtained on the basis of matrix methods are subdivided into characteristics of the 1 st and 2 nd orders.

The construction of first order texture characteristics implies that the structural elements are individual pixels of the original image [9]. Gray-Level Matrix (GLM) is a vector of frequencies of gray level occurrence in the processed image 𝐼(𝑥, 𝑦) with the size 𝐿 𝑥 × 𝐿 𝑦 :

GLM(𝑘) = # {(𝑥, 𝑦)|𝐼(𝑥, 𝑦) = 𝑘, (𝑥, 𝑦) ∈ 𝐿 𝑥 × 𝐿 𝑦 }

where # means the number of elements in the set, 𝑘 = 1, . . . , 𝑁 𝑔𝑙 , and 𝑁 𝑔𝑙 is the number of gray scales.

For building the Gray Level Difference Matrix (GLDM), the original image 𝐼 is converted into a difference image 𝐷𝐼:

𝐷𝐼 Δ𝑥,Δ𝑦 = |𝐼(𝑥, 𝑦) − 𝐼(𝑥 + Δ𝑥, 𝑦 + Δ𝑦)|

where parameters Δ𝑥 and Δ𝑦 are displacements along the horizontal and vertical directions, respectively. GLDM is a vector of frequencies of occurrence of absolute values of differences in gray levels at the given displacement:

GLDM(𝑘) = # {(𝑥, 𝑦)|𝐷𝐼(𝑥, 𝑦) = 𝑘, (𝑥, 𝑦) ∈ (𝐿 𝑥 − Δ𝑥)×(𝐿 𝑦 − Δ𝑦)} , 𝑘 = 0, . . . , 𝑁 𝑔𝑙 −1.

An example of constructing GLM and GLDM matrices is shown in Figure 1. For calculating extract texture features, GLM and GLDM are converted into corresponding probability mass function estimates

𝐹 GLM (𝑘) = GLM(𝑘) 𝑁 𝑔𝑙 ∑︀ 𝑘=1 GLM(𝑘)

, 𝐹 GLDM (𝑘) = GLDM(𝑘)

𝑁 𝑔𝑙 −1 ∑︀ 𝑘=0 GLDM(𝑘)

.

The corresponding 1 st order texture characteristics are presented in Table 1.

The extraction of second order texture characteristics is primarily associated with the construction of two-dimensional distributions. Structural elements in this case consist of two pixels or two groups of pixels, for each of which a corresponding characteristic is determined. The best known is the method proposed in [10]. The method uses structural elements consisting of two

𝐼(𝑥, 𝑦)

Expectation

𝑁 𝑔𝑙 −1 ∑︀ 𝑘=0 𝑘𝐹 GLDM (𝑘) Mean-square 1 𝐿 𝑥 𝐿 𝑦 𝐿𝑥 ∑︀ 𝑥=1 𝐿𝑦 ∑︀ 𝑦=1 𝐼 2 (𝑥, 𝑦) Contrast 𝑁 𝑔𝑙 −1 ∑︀ 𝑘=0 𝑘 2 𝐹 GLDM (𝑘) Entropy − 𝑁 𝑔𝑙 ∑︀ 𝑘=1 𝐹 GLM (𝑘) log 𝐹 GLM (𝑘) Angular Second Moment 𝑁 𝑔𝑙 −1 ∑︀ 𝑘=0 𝐹 2 GLDM (𝑘) Energy 𝑁 𝑔𝑙 ∑︀ 𝑘=1 𝐹 2 GLM (𝑘) Entropy − 𝑁 𝑔𝑙 −1 ∑︀ 𝑘=0 𝐹 GLDM (𝑘) log 𝐹 GLDM (𝑘) Variance 𝑁 𝑔𝑙 ∑︀ 𝑘=1 (𝑘 − 𝜇) 2 𝐹 GLM (𝑘)

pixels at a certain specified distance (adjacency distance). One of these pixels is called reference.

For the reference pixel, the neighboring one is selected in a given direction of adjacency. To describe the spatial relationship between the reference and neighboring pixels, the frequencies of occurrence of the corresponding gray-scale pairs are calculated for all possible positions of the reference pixel in the original image. Based on these frequencies, we can form a matrix known as the Gray-Level Co-occurrence Matrix (GLCM) or Spatial Gray-Level Dependency Matrix (SGLDM). An example of constructing the GLCM is shown in Figure 2. The GLCM is a square matrix containing integer values. The size of GLCM is determined by the number of gray levels. So in the example presented in Figure 2, the original image has 8 gray levels and, respectively, the GLCM has a size of 8 × 8. The matrix has a symmetrical appearance if you do not take into account the order of the grayscale in the reference and neighboring pixels. The original image 𝐼(𝑥, 𝑦) is a function of two spatial coordinates, so for each pixel we can calculate the differential estimate of the gradient function. The magnitude of the gradient

𝑔(𝑥, 𝑦) = √︃ (︂ 𝜕𝐼 𝜕𝑥 )︂ 2 + (︂ 𝜕𝐼 𝜕𝑦 )︂ 2

characterizes the rate of changing tone of the image in reference pixels. To obtain a difference estimate, we use the Sobel operator:

𝑆(𝑥, 𝑦) = √︁ 𝑆 2 𝑥 (𝑥, 𝑦) + 𝑆 2 𝑦 (𝑥, 𝑦) ≃ 𝑔(𝑥, 𝑦)

where

𝑆 𝑥 (𝑥, 𝑦) = [𝐼(𝑥 + 1, 𝑦 − 1) + 2𝐼(𝑥 + 1, 𝑦) + 𝐼(𝑥 + 1, 𝑦 + 1)] − − [𝐼(𝑥 − 1, 𝑦 − 1) + 2𝐼(𝑥 − 1, 𝑦) + 𝐼(𝑥 − 1, 𝑦 + 1)] , 𝑆 𝑦 (𝑥, 𝑦) = [𝐼(𝑥 − 1, 𝑦 + 1) + 2𝐼(𝑥, 𝑦 + 1) + 𝐼(𝑥 + 1, 𝑦 + 1)] − − [𝐼(𝑥 − 1, 𝑦 − 1) + 2𝐼(𝑥, 𝑦 − 1) + 𝐼(𝑥 + 1, 𝑦 − 1)] .

Thus, by specifying the number of gradient gradations to be equal 𝑁 𝑔𝑙 , we can build an image of modules of brightness gradients

𝐺(𝑥, 𝑦) = int [︂ 𝑆(𝑥, 𝑦) − 𝑆 min 𝑆 max − 𝑆 min ]︂ • 𝑁 𝑔𝑙

in pixels of the original image. Gray Gradient Co-occurrence Matrix (GGCM) [11] is built in the same way as GLCM, only for the image 𝐺(𝑥, 𝑦).

The estimate of the probability mass function of the co-occurrence of a given number of gray levels can be obtained as the normalized GLCM 𝑝(𝑖, 𝑗) = GCLM(𝑖, 𝑗) The detailed description all of them is presented in [12]. It should also be noted that in work [13] it is stated that for most practical tasks, it is sufficient to use 5 of them, which are given in Table 2. The necessary marginal expectations and marginal STD can be calculated as:

𝑁 𝑔𝑙 ∑︀ 𝑖,𝑗=1𝜇 𝑖 = 𝑁 ∑︁ 𝑖=1 𝑁 ∑︁ 𝑗=1 𝑖 • 𝑝(𝑖, 𝑗), 𝜇 𝑗 = 𝑁 ∑︁ 𝑖=1 𝑁 ∑︁ 𝑗=1 𝑗 • 𝑝(𝑖, 𝑗), 𝜎 𝑖 = ⎯ ⎸ ⎸ ⎷ 𝑁 ∑︁ 𝑖=1 𝑁 ∑︁ 𝑗=1 (𝑖 − 𝜇 𝑖 ) • 𝑝(𝑖, 𝑗).

Texture segmentation of panchromatic satellite images is based on the moving window method. The moving window is a rectangular contour selecting the analyzed part of the image under processing. The size of the window is determined by the characteristic scale of recognized textures. If the window size is chosen too small, the result of the texture classification will represent the high frequency noise. On the other hand, too large size of the window leads to excessive smoothing of the contours of recognized objects. The center of the window runs through all the points of the panchromatic image. To reduce the amount of computation in practical tasks, it is sufficient to run only pixels whose coordinates correspond to the pixel centers of the joint multispectral image, when processing panchromatic and multispectral images together.

To carry out the supervised classification texture features we employed an ensemble algorithm known as Error Correcting Output Codes (ECOC). The algorithm is designed to formalize the responses of binary learners as the multiclass classifier based on some results of the theory of information and coding. Le the algorithm of binary classification is the Support Vector Machine (SVM) with the Gaussian kernel. The response of the SVM algorithm is the classification score which means the normalized distance from the classified sample to the discriminant surface in the area of the relevant class. The coding stage of the ECOC algorithm consists of calculating classification scores or each of SVM binary learners defined by one-versus-one coding design matrix and corresponded Hinge binary losses for each of the considered classes. The decoding stage consists of the selection of the class corresponding to the minimum average loss. For optimizing the feature space arises, the regularized forward selection method is used. The method has better stability in comparison with the standard greedy selection algorithm suffering from the high sensitivity of the selected optimal sequence of features to small changes in the training set. On the other hand, the regularized algorithm requires more computational costs.

The classification quality was assessed by confusion matrix (CM) and the related parameters, such as total error (TE), total omission error (TOE) and total commission errors (TCE). CM is the basic classification quality characteristic allowing a comprehensive visual analysis of different aspects of the classification method used. Rows of CM represent reference classes and columns represent predicted classes. TE is defined as the amount of incorrectly classified samples over the total number of samples. TOE is the mean omission error over all classes considered, where the omission error is the amount of false classified samples of the selected class over all samples of this class. TCE is the mean commission error over all possible responses of the classifier used, where the commission error is defined as the probability of false classification for each possible classification result.

Results and discussion

During the joint spectral-texture processing of satellite images, texture features are usually employed for solving the following two tasks: segmentation of the contours of natural and manmade objects (including the selection of building zones and forest areas), and classification of structural parameters of the forest canopy. Thus, for carrying out numerical experiments using the above methods, we selected two relevant test plots in WorldView-2 panchromatic images with the spatial resolution ∼ 0.5 m. The first test plot, which hereinafter referred to as Konstantinovsky, is located on the territory of the Savvatyevskoe forestry (Tver region) near the Domnikovo village. The plot contains several large zones corresponding to 5 types of objects of varying complexity having well differing texture: water surface (Konstantinovsky sand quarry), pine forest, building zone, field and peat swamp forest. The RGB image of the test plot, the corresponding expert map of objects and the results of texture analysis are shown in Figure 3.

The classifier used is sensitive to the difference in the number of training samples for the considered classes. Increasing the number of training samples for one of the classes leads to increasing the prior probability of its classification. Thus, in order to avoid this problem, we used a balanced training set containing 500 samples for each class. Remaining data (testing set) were used for independent validation. Also, since the accuracy of texture classification essentially depends on the size of the moving window, we carried out a series of calculations, which allowed us to determine the range of acceptable size values. Thus, we used the moving window with the size of 109 pixels in horizontal and vertical directions. The original image has been reduced to 64 gray levels. The described above feature optimization was used for the GLCM and GGCM methods to avoid the curse of dimensionality problem. For the classification of texture features obtained by the GLM and GLDM methods, we used a full set of features.

Estimates of total characteristics of classification quality obtained from training (resubstitution method) and testing sets (independent validation) are represented in the Table 3. We can see that the GLCM method provides the most accurate results. The total error is about 1%. At that the total omission and commission errors are very close. It should be noted that difference between dependent and independent estimates of the error is insignificant, which indicates a good generalization ability of the trained ECOC SVM classifier. The total errors of the GGCM and GLM methods are significantly higher, but remain at an acceptable level. It should be noted that a visual comparison of the classification results presented in Figure 3 shows that GGCM reproduces the expert map of objects much better and contains significantly less noise in comparison with GLM. Table 4 contains class-wise omission and commission classification errors for Konstantinovsky test plot. The best classification result corresponds to water surface. We can see that all the methods reveal a high accuracy. The building zone and field are classified with acceptable level of errors. The worst accuracies correspond to conifer and peat swamp forest stands, however it low enough for GLCM and GGCM methods. GLDM demonstrates worst results and it cannot be used texture segmentation of forest areas.

The second test plot, hereinafter referred to as GFP Dementyev, is located on the territory of the Bronnitskoe forestry near the Lubninka village. The plot is part of the territory of the geographical forest plantations (GFP) of the forester P.I. Dementyev. The RGB image of GFP Dementyev and the corresponding expert classification map are shown in Figure 4. The choice of this site is due to the large variety of plantations with different structures. By the variety of species, the stands of the Bronnitskoe forestry cover the main forest-forming species of Russia. From the 1950s to the present, various species and ecotypes of larch, which are grown here outside the natural habitat, have been tested in this area. The forest canopy of the test site contains 7 visually noticeable texture classes: 1 -mixed conifer stand (larch, pine and spruce) with a dense canopy and high values of density; 2 and 3 -agricultural areas of different structure; 4 -deciduous stand with a predominance of birch, dense canopy and the relative stocking of 0.9; 5 -mixed birch stand with the relative stocking of 0.8; 6 -mixed birch stand with a pronounced cluster structure of the canopy; 7 -cultivated plantations of larch with a regular structure. In this case GLCM and GGCM methods show similar results. Total classification error is about 4%. Accuracy of GLCM seems to be a little higher compare to GGCM, however we can see for this method that difference between resubstitution and independent estimates is also more significant than for GGCM. GLM and GLDM method demonstrate weak classification results. Analyzing class-wise errors presented in Table 5, we can see that the regular structure of larch stands corresponds to minimum errors, about 2% for GLCM and GGCM methods. This result is confirmed also by Figure 4. The highest errors correspond to the dense deciduous stand, however it is explained by small number of pixels corresponding to this object.