AUTOMATIC APPROACH TO VHR SATELLITE IMAGE CLASSIFICATION

In this paper, we present a proposition of a fully automatic classification of VHR satellite images. Unlike the most widespread approaches: supervised classification, which requires prior defining of class signatures, or unsupervised classification, which must be followed by an interpretation of its results, the proposed method requires no human intervention except for the setting of the initial parameters. The presented approach bases on both spectral and textural analysis of the image and consists of 3 steps. The first step, the analysis of spectral data, relies on NDVI values. Its purpose is to distinguish between basic classes, such as water, vegetation and non-vegetation, which all differ significantly spectrally, thus they can be easily extracted basing on spectral analysis. The second step relies on granulometric maps. These are the product of local granulometric analysis of an image and present information on the texture of each pixel neighbourhood, depending on the texture grain. The purpose of texture analysis is to distinguish between different classes, spectrally similar, but yet of different texture, e.g. bare soil from a built-up area, or low vegetation from a wooded area. Due to the use of granulometric analysis, based on mathematical morphology opening and closing, the results are resistant to the border effect (qualifying borders of objects in an image as spaces of high texture), which affect other methods of texture analysis like GLCM statistics or fractal analysis. Therefore, the effectiveness of the analysis is relatively high. Several indices based on values of different granulometric maps have been developed to simplify the extraction of classes of different texture. The third and final step of the process relies on a vegetation index, based on near infrared and blue bands. Its purpose is to correct partially misclassified pixels. All the indices used in the classification model developed relate to reflectance values, so the preliminary step of recalculation of pixel DNs to reflectance is required. Thanks to this, the proposed approach is in theory universal, and might be applied to different satellite system images of different acquisition dates. The test data consists of 3 Pleiades images captured on different dates. Research allowed to determine optimal indices values. Using the same parameters, we obtained a very good accuracy of extraction of 5 land cover/use classes: water, low vegetation, bare soil, wooded area and built-up area in all the test images (kappa from 87% to 96%). What constitutes important, even significant changes in parameter values, did not cause a significant declination of classification accuracy, which demonstrates how robust the proposed method is. * Corresponding author


INTRODUCTION
The most popular satellite image classification approaches are not fully automatic.Regardless of the specific method of classification: spectral, contextual or object-based, the intervention of a human operator is required at some stage of the process, either in the preliminary stagedefining class signatures in the supervised approach (and defining the segmentation process in object-based classification), or in the final stagein the interpretation of results in the unsupervised approach.This intervention is very often required due to the different conditions of the images which are classified (e.g.different terrain illumination conditions, or a different growing season, but also differences in scanner calibration).
In this paper, we present a proposition of a maximized automatic classification of VHR satellite images irrespective of time and of the conditions of image acquisition.Unlike the most widespread approaches, the proposed method requires no human intervention, except the setting of the initial parameters.It can be achieved by the use of indicesincluding some well known ones, such as the normalized difference vegetation index (Rouse, 1974), or others, developed specifically for this method: the index based on the result of granulometric analysis of the texture of an image (Kupidura, 2015) allowing to separate, distinctly, different land cover classes in regard to their spectral or textural characteristics.
Current and reliable information on land cover is still an important element in decision-making in terms of spatial planning.Each new industrial investment, e.g.transmission or transport lines covering a large area, necessitates the acquisition of the most updated data at the required thematic and geometric level of detail.Therefore, the motivation for this study was the need to obtain, in a short period of time, reliable information about the current land coverage for large-area development, with detail to allow for a preliminary estimate of the cost of the planned investment.Thus, the selection of classes was imposed, on the one hand by the diversity of the existing forms of coverage, and on the other hand, by the need for an overall presentation of information for a large area, acquired in a short period of time.Hence, the need for maximum automation, simplification and the development of a universal model for the separation of land cover classes: water, forest, low vegetation, bare soil and built-up area.

Fundamentals of the proposed model
The main assumption of the model is that a combination of both the spectral and textural characteristics of the classes mentioned may be sufficiently separative and allow to distinguish pixels representing these classes with satisfying efficiency.Another important assumption is to base the model on singular indices, rather than groups of parameters, to make the model as simple as possible and also to allow, if needed, to easily change the parameters of classification to obtain optimal results.
Basing on these initial terms, we propose the normalized difference vegetation index (NDVI), presented by Rouse et al. (1974) as the main spectral index.Other studies (Kupidura, 2002;Thanapura et al., 2006;Jia et al., 2014) confirm that NDVI values help to separate basic land cover types, such as water, vegetation and non-vegetation, and also, they are partially resistant to differences in illumination conditions.
As the source of textural information, we propose a granulometric analysis of an image.Its efficiency has been researched and proven previously by several researchers (Mering et al., 1996;Mering and Chopin, 2002;Mering et al., 2010;Mering, 2011, Kupidura, 2010, 2015;Kupidura et al., 2010, Szeszko, 2014).Granulometric analysis is based originally on a sequence of morphological openings (Haas et al., 1967) using structuring elements of sequentially increasing size.However, a closing operation might also be used in the same way; the research by Kupidura (2015) shows that both types of granulometry are equally useful, and often are complimentary, for an analysis of single band (monochromatic) satellite images.Granulometric analysis may be applied locally, producing series of granulometric maps (Dougherty et al., 1992) images containing information on the character of texture in a defined neighbourhood.One of the important advantages of granulometric maps (and granulometric analysis in general) is that it is not susceptible to the border effect (Ruiz et al., 2004;Kupidura, 2015), i.e. qualifying borders of objects in images as spaces of high texture, affecting other methods of texture analysis like GLCM statistics, discrete wavelet transformation or fractal analysis.This feature is one of the causes of the relatively high effectiveness of morphological granulometry as a source of textural information which might be used in image classification.Another important advantage of this method is its multiscalityeach granulometric map presents information on a different size of texture element.
The third index has been introduced as a result of the performed tests.It is based on the Simple Ratio (Rouse et al., 1974) formula, using near infrared and blue reflectance values (instead of a red one).According to Dupigny-Girou and Lewis (Dupigny-Girou et al., 1999), a combination of these spectral bands helps to distinguish some forms of land cover, such as urban areas, water bodies, and full-leaf versus leafless conditions, better than standard vegetation indices (based on near infrared and red bands), and its efficiency in this subject matter has been confirmed in our research.The results of research by Dupigny-Girou (Dupigny-Girou, 2007) indicate that the NIR/blue ratio discriminated among moisture-driven land cover / land use variations, wetland types and moisture stress extent.
To achieve a maximum universality of the model, all indices have been evaluated for satellite images subjected to radiometric correctionrecalculation of digital numbers to Top-of-atmosphere (TOA) reflectance values.

Test area, images and methodology of the model development
The presented model has been proposed for the VHR satellite images for 5 basic land cover classes: (1) "Water"inland water: rivers, lakes, water bodies, (2) "Forest"wooded areas, shrubs, orchards and plantations, (3) "Low vegetation"meadows and pastures, arable land with vegetation, (4) "Bare soil"arable land without vegetation, (5) "Built-up area"built-up areas and lines of communication.
The structure of the model being developed has been evaluated based on WorldView-2 and Pleiades satellite images (T1) and then tested on other Pleiades images (T2, T3 in table 1).The values of considered indices and also their formulas (likewise in the case of granulometric analysis) have been evaluated and tested, as well.The final form of the proposed model is the result of a somewhat iterative process considering the results of the test data.2.
The performed tests and analysis allowed to establish the final form of the classification model as well as to evaluate the optimal values of indices.

Classification model
Basing on the research performed on test images, we proposed the following model of classification, presented in figure 3. The model consists of 3 steps.The first step, an analysis of spectral data, relies on NDVI values.Its purpose is to distinguish basic classes, such as water, vegetation and non-vegetation, which all differ significantly spectrally, and can thus be easily extracted basing on spectral analysis.Previous analyses allowed to evaluate V1 and V2, values of NDVI: V1 distinguishing water and non-vegetation (i.e.bare soil and built-up area), and V2 distinguishing non-vegetation from vegetation (i.e.forest and low vegetation).This step allowed to extract 3 groups of classes (according to the model -"GROUPS").
The second step relies on granulometric maps.The goal of this step is to distinguish classes previously grouped based on NDVI values: low vegetation and forest, and bare soil or built-up area.
Due to the different texture form of high texture classes, i.e. forest and built-up area, significant textural information is provided by different granulometric maps for different classes.Thus, values G1 and G2 refer to two indices based on values of granulometric maps: from C1 to C7 (closing) and from O1 to O7 (opening).This step allowed to distinguish 5 preliminary classes within the groups (see figure 3).
The third and final step again relies on spectral values.The NIR/B index helped to distinguish certain parts of forest from built-up area.Because of low NDVI values (caused by low NIR values), the shaded parts of forest were assigned to the nonvegetation group in the first step of the process, and due to high texture, they were classified as built-up area.The NIR/B index and the B1 parameter allow to correct the assignment of these misclassified pixels.This step allowed to fix the final classes.As it is shown in the model scheme (figure 3), each class shall meet a set of conditions presented below.
1. 'Water"this is the only class in Group 1, and therefore all the pixels belonging to this group are assigned to Class 1.The condition for belonging to Class 1 can therefore be determined as follows: GROUPS = 1 2.
"Forest"it is a class of Group 2 (it also contains class 3 "low vegetation") which is characterized by a stronger, more heavy texture than Class 3 (see figure 2.2).In order to distinguish it from class 3, four granulometric maps are used after opening and closing with the lowest index C1, C2, C3, C4, O1, O2, O3, O4.The "Forest" Class is represented by pixels with high values on those granulometric maps, where the sum of all 8 maps exceeds the limit described as G1.The condition of pixels belonging to Class 2 can be determined as follows: GROUPS = 2 AND C1+C2+C3+C4+O1+O2+O3+O4 > G1 3. "Low vegetation"is a class in Group 2 characterized by a weaker texture and, therefore, lower than Class 2 values on the granulometric maps (compare figure 2.3).In order to distinguish it from Class 2, the same eight granulometric maps are used: C1, C2, C3, C4, O1, O2, O3, O4.The "low vegetation" class is represented by pixels with high granulometric values on these maps, where the sum of all 8 maps does not exceed the limit value, defined as G1.The condition of pixels belonging to Class 3 can be determined as follows: GROUPS = 2 AND C1+C2+C3+C4+O1+O2+O3+O4 <= G1 4.
"Bare soil"this is the class of Group 3 characterized by less heavy texture, and hence lower values on granulometric maps compared to the second class of this group -"Built-up area" (Class 5).Since the difference in the texture of these two similar spectral classes is most evident on a slightly smaller scale, granulometric maps with higher indexes must be used: C4, C5, C6, C7, O4, O5, O6, O7 .The "Bare soil" class is represented by pixels with low values on those granulometric maps, where the sum of all 8 maps is less than the limit defined as G2.The condition of pixels belonging to Class 4 can be determined as follows: GROUPS = 3 AND C4+C5+C6+C7+O4+O5+O6+O7 < G2 5. "Built-up area"this is another class of Group 3 which has a clear, strong texture, and consequently, high values on granulometric maps, due to the size of the texture grainthat is on those maps with higher indexes.According to the proposed solution, the "Built-up area" Class is represented by pixels, whose sum of the values in all 8 maps is not less than the limit defined as G2.The condition of pixels belonging to Class 5 can be determined as follows: GROUPS = 3 AND C4+C5+C6+C7+O4+O5+O6+O7 >= G2 As observed on the test images in the case of deciduous forest with a specific form i.e. not fully spring, non-consistent canopy or sparse forest, there may be tiny shaded areas.These surfaces belong, according to the specification of classification, to Class 2 -"Forest", but are characterized by low values of NDVI and high values on granulometric maps.Consequently, in the second step, it can be assigned to Class 5 "Built-up area".For this type of surface, a condition relating to the value of the quotient indicator NIR / B was introduced.This surface receives higher values than actual built-up areas.Pixels classified as built-up area (Class 5), but having high values of this ratio and exceeding the limit valuedefined as B1will be assigned to Class 2. The condition of changes from Class 5 to Class 2 can be determined as follows: CLASSES1 = 5 AND NIR/B > B1
The proposed methodology of image processing has been subjected to thorough testing in many variants and its accuracy has also been assessed.As a result, the elaborated model prepared basing on a specific set of satellite data for test area No. 1 might be tested and calibrated in the image from the same sensor acquired at other times for test area No. 2. In total, 15 variants of the model with different parameter settings of V and G were tested for test area No. 1, and 30 variants for test area No. 2. Area 2 is much more diverse within a single class than area No. 1 (figure 1) and required many additional analyses.
Assessment of the accuracy was performed according to the methodology proposed by Congalton and Green (Congalton and Green, 2009) based on the estimation of the thematic accuracy for the classification results.Assessments were made on the basis of the test sites prepared, which constitute the statistical sample representative for each class.Test fields for each image (classification results) were created on the basis of an independent photointerpretation of multispectral color composition guided by the fundamental principles of sample determination.For test area No. 1, 29 test fields (polygons) were prepared with a total area of over 61 hectares, and for test area No. 2, two variants were prepared: 21 and 23 polygons depending on the date (in total over 200 ha).
For test area No. 1, very good results for all the defined classes of land cover were achieved.Depending on the model variant, and on the set of input parameters, the overall classification accuracy varies border omega from 90% to 97% and kappa from 0.86 to 0.96.The highest producer's accuracy was achieved for the "water" class and it amounts to 99.9%.The weakest is classified in "built-up areas", and in this case, the producer's accuracy ranges from 43% to 82% depending on the model variant employed.The user's accuracy of "water" is 100%, but for other classes it is in the range of 63% to 100%.The lowest results obtained with successive variants were used for the calibration of the model until a satisfactory quality level was achieved.The best results were obtained for the parameters G1 ϵ (40; 50) and G2 ϵ (20; 25).
For test area No. 2, measures of overall assessment of the accuracy indicate relatively good quality classifications: omega from 84% to 98% and kappa from 0.74 to 0.97.The high efficiency classification of "water" has been shown at 97% of the producer's accuracy and 100% of the user's accuracy.For the "built-up area" class, the producer's accuracy at 85% and the user's accuracy level of 95% has been achieved, depending on the model variant, these values are slightly different.The class associated with vegetation obtained results significantly exceeding 90%, both in terms of producer's and user's accuracy.The problematic class is the "bare soil" class, which is often mixed with the class "meadows, pastures and arable land covered with vegetation."It is a result of the fact that both classes are similar in texture and are difficult to distinguish clearly in the image at the end of March, when the growing season has just started and may be more nuanced in moisture in the case of arable land, meadows and pastures.The best results were obtained for the parameters G1=50 and G2 ϵ (40; 50).
The best classified land cover class is "water" irrespective of the area or spatial resolution images.The only obstacle in the classification of water bodies are reflections occurring on the surface forming a very irregular bright patch, which makes it difficult to correct separation.However, this is relatively rare.Other phenomena that can affect the quality of water classification is the freezing of the water surface in early spring or overgrowing by the water plants during the summer.
The weakest classified class is that of "built-up areas and communication lines", due to its large variety of functionality (single-family houses, factories, roads and other infrastructure elements), of structures (compact or distributed buildings, linear objects) and of the spectral properties of the materials from which it is composed.For classes "forests and wooded land", "the meadows, pastures, arable land with vegetation" and "bare soil" classification results are good, but this depends on the appropriate selection of the parameters.Hence, the parameters' setting by the operator has been proposed.
The studies determined that the following parameters: V1=-0,10; V2=0,25; G1=37; G2=22; B1=2,0 (table 2) allow to achieve very good results of automatic classification according to the model presented above.2. Optimal values for V1, V2, G1, G2 parameters These values have proved to be optimal for two of the three used test images (T1 and T2), but also for the third test image (T3).They have allowed to obtain high classification accuracy, and therefore can be represented as optimal, at least on the basis of the data test used.Nonetheless, it is worth paying attention to the possibility of increasing the potential effectiveness of classification by introducing the possibility to modify these values for images with different characteristics.This model will be implemented into the software which allows to evaluate the potential cost related to the purchase of certain form of land cover for industrial investment purposes.

CONCLUSIONS
The research allows to determine the optimal model form and indices values.Using the same parameters, we have obtained very good accuracy of extraction of 5 land cover/use classes: water, low vegetation, bare soil, wooded area and built-up area in all test images (kappa from 0.87 to 0.96).Importantly, even significant changes in parameter values did not cause the significant declination of classification accuracy, which is a clear sign of high robustness of the proposed method.Due to the simplicity of the model and the fact that it relies on singular indices values, it can be easily equalized (adjusted) in order to obtain optimal accuracy of classification.
The developed model of delimitation of particular land cover classes takes into account the need to define the parameters − values that allow to distinguish the different classes.The fully automated process, in the strict sense of the word, and therefore not allowing any operator intervention, should be required to definitively determine those values.In the present case, this kind of test can also be taken.It should be noted, however, that despite work for far-reaching unification of potential source images, they may still differ from each other in either spectral or textural characteristics, due to factors such as the acquisition date, the height of the sun, weather conditions, the state of the atmosphere, the state of vegetation, etc.As a result of this diversity, acceptance of universal parameters for all the images elaborated using this model cannot fully exploit the potential of the model, and consequently the data used.It should be emphasized that the form of the model is designed in such a way that suppression of the objective differences in the processed images discussed above can be achieved solely by changing the individual parameters.

For
research work, test areas were chosen in the territory of Poland and characterized by the typical rural and suburban landscape (figure 1).Test area No. 1 is located in the south of Warsaw, covering the southern part of the Wilanów district and the northern part of the municipality Konstancin; test area No. 2 is located in the province of Western Pomerania.1. 2.

Figure 1 .
Figure 1.Test areas: No. 1suburban area of Warsaw, No. 2rural landscape of Western Pomerania.The research on the formulas and values of certain indices consisted in analyzing spectral (for NDVI) and textural (granulometric maps) data regarding different examples of 5 land cover types being the subject of the classification.Different types of granulometric analysis based on different types of opening and closing operations (such as simple operations and operations with a multiple structuring element) were tested regarding separability of different classes on specific granulometric maps.

Figure 2
Figure 2 presents samples of land cover classes visualized on a spectral colour composition (column A) and a colour composition consisting of different granulometric maps based on operations with multiple structuring elementcolumn B. A B

Figure 3 .
Figure 3. Block diagram of the classification model.Number of classes according to the list presented in section 2.2.

Table 1 .
Test images