The new concept of natural numerical networks is introduced in this paper. It is used as a novel tool for automated classification of protected habitats using satellite images. The natural numerical network is based on the numerical solution of the nonlinear forward-backward diffusion (FBD) equation on complete graph. The network, represented by the numerical
discretization of the FBD equation, can be classified as a new deep learning method. Usually, a deep learning method is formed by an artificial neural network with many hidden layers. In our discretization scheme, we use a sequence of time steps resolving the dynamics of the diffusion equation, and one hidden layer corresponds to one time step of the numerical scheme for solving the FBD equation. The proposed deep learning method does not use artificial neural network principles in its construction, see e.g.1990JGCD , BCC or Goodfellow-et-al-2016 . In addition, diffusion equations are widely used in modelling phenomena in natural sciences, such as biology, physics or chemistry; thus, we call the proposed network natural. Indeed, the method seems to truly be a natural procedure for the clustering and supervised classification of data, so we call the proposed method the natural numerical network or the natural network for short. In building the natural network, we are inspired by the work HR in which Eldad Haber and Lars Ruthotto showed the relation between a successful deep learning model, the so-called Residual Neural Network (ResNet) RN ; ReversArch
, and the numerical solution of the system of ordinary differential equations using the forward Euler method. Subsequently, they designed parabolic and hyperbolic networks for deep learning based on the appropriate partial differential equationsHR . Our natural network uses another type of PDEs, the nonlinear forward-backward diffusion equations. The forward diffusion causes the movement of points in a feature space toward each other. The opposite effect, when the points are kept away from each other, is caused by backward diffusion. This yields the desired classification. The natural network based on FBD equations contains just a few parameters that are optimized in the learning phase of the method. After learning the parameters and optimizing the topology of the graph of the network, the classification is performed. The relevancy map for each habitat is created as a novel tool for validating the classification, studying the relation of the habitat classification with the species composition and finding new Natura 2000 habitat appearances. It shows abilities of developed method in distinguishing between next-standing mixed deciduous forests with similar species composition and between two types of riparian forests as shown in the alluvium of the Danube river, and also potential in automated finding new localities of protected habitat areas as shown on softwood floodplain forest newly discovered in Slovakia by the proposed method.
The use of satellite data has become one of the essential methods for effectively and directly acquiring information on the Earth’s surface Liu2015 ; Randin_etal.2020 . Together with standardized botanical records (plots) and regular in situ measurements, remote sensing is a powerful monitoring instrument Corbane_etal.2015 ; Lausch_etal.2020 playing an irreplaceable role in acquiring data essential for evaluating and implementing environmental policy by data analysis Liu2015 ; Ullerud_etal.2018 ; Chi_etal.2016 . Remote sensing is also one of the most important tools in ecology and nature conservation for achieving the effective monitoring of ecosystems in space and time Rocchini ; thus, using satellite images to monitor habitats and biome dynamics has been highlighted in many types of research activities. Ecosystems representing Natura 2000 habitats are complex plant communities including tree, shrub, and herb layers together with typical fauna EEANatura2000 ; Natura2000 for which it was impossible to reach an automated identification and classification with existing methodologies based on satellite data Bock ; Vanden so far. In this paper we develop method which performs the first step to fill this gap. Although the method is generally designed to work with any type of optical data monitoring the Natura 2000 habitats, we use the optical information from spectral bands of the Sentinel-2 satellite ESASentinel freely available on European Space Agency (ESA) servers.
In summary, the aims of presented study are 1) to give the complete mathematical definition of the natural numerical network based on forward-backward diffusion equations, 2) to present the numerical scheme corresponding to the natural numerical network in all phases of the classification algorithm, 3) to optimize the natural numerical network on training dataset of Natura 2000 habitats, and 4) to apply the trained natural numerical network to classification of Natura 2000 habitats, construction of the habitat relevancy maps for results validation and for finding new appearances of protected habitats.
2.1 Mathematical model
Let us define a graph as an ordered pair, where is a finite set of vertices, and is a set of the two-element subsets of representing the edges of the graph teoriaGrafovKnor . We denote the number of vertices of the graph by . Let us consider that the graph is a complete graph that means that each vertex is connected to each other vertex by an edge. Let us suppose that graph is undirected and thus the edges do not have an orientation.
Let us consider the function representing the spatial coordinates of the vertex in time . In our case, is a dimension of feature space . A diffusion of on the graph is formulated as a partial differential equation (PDE)
where represents a so-called diffusion coefficient, see also calculusOnGraphs . We consider equation (1) together with an initial condition , . The boundary conditions are not necessary to prescribe because diffusion occurs between all vertices of the complete undirected graph .
Let us define the distance between two vertices as the Euclidean distance between two points and in the feature space and denote it . Since every pair of vertices of , and , forms one edge, , we can simplify the notation for the distance and denote it as a length of the edge , which will be used throughout the following text.
We will design the diffusion coefficient depending on the length of the edges of the graph . It will give a nonlinear diffusion model on the graph representing a generalization of the Perona-Malik model from the image processing Perona-Malik . We consider equation (1) with the diffusion coefficient in the form
The value of in the diffusion coefficient depends on the type of diffusion which is applied on every single edge. If we need to apply forward diffusion, we choose as a positive constant. The forward diffusion, represented by a positive diffusion coefficient, averages the values of a diffused quantity. The averaging reflects the smoothing property of the standard diffusion equation. It is called forward because it describes diffusion towards the future and in our application, it causes a moving and thus clustering of points together. On the other hand, the backward diffusion is represented by a negative diffusion coefficient, is a small negative value, and can be understood as returning to the past in a diffusion process. It is an inverse process of the smoothing (averaging) values of a diffused quantity and in our application, it gives a repulsion of the points belonging to different clusters. If we only used the backward diffusion model, the points would be moving away from each other, and the whole system became unstable. But by a suitable combination of the forward and backward diffusion, when we choose a small negative coefficient
for backward diffusion, we do not observe any calculation instability. Such a model is a suitable and natural tool for supervised learning - the points inside the given clusters are kept together while points of different clusters are kept away from each other. Such proper behaviour is realized by the model (1)-(2).
In Fig. 1 we illustrate behaviour of the model (1)-(2) on three given clusters where by blue lines we plot some of the links of forward diffusion and by red lines some of the links of backward diffusion. This figure depicts the basic features and behaviour of the natural network. The points inside a given cluster are attracted by the forward diffusion while there is a repulsion of points of different clusters by the backward diffusion.
Additionally, in Fig. 2 we illustrate the situation arising in supervised learning and application phases of the classification method when a new observation is added into the network. Only the forward diffusion is applied to all links of the vertex representing the new observation. It is depicted by the blue lines connecting the new observation (black square) with every other point. Thus, this new observation is attracted by a certain diffusion speed to all existing clusters which themselves are subject to the forward-backward diffusion as described before. The dynamics of the network decides about the cluster membership of the new observation.
The model (1)-(2) contains, together with forward-backward diffusion switch , also a weighting coefficient . The constant controls how the length of the edge affects the diffusion of the vertices and over time. If is large, the diffusion coefficient is close to 0 which means that the diffusion process will be slow, and the points are not diffusing (moving to each other) by an averaging. If is small, the diffusion coefficient is close to 1, the diffusion process is faster and points are moving to each other fast by the diffusion. Since the coefficient is multiplying the squared length of the edge, the distant points in the feature space are averaging (moving) slower than the close points.
The model (1)-(2), after discretization, represents a basic natural numerical network for supervised deep learning classification and, with just a positive coefficient , it is also a proper model for unsupervised clustering (which is, however, not discussed in this paper). This basic model allows useful modifications which will be used also in our final classification algorithm. First of all, we slightly modify the diffusion coefficient into the form
Now, the parameters , , represent weights for each coordinate ,
, of the vector
By this modification, we can control the diffusion speed in each direction of the -dimensional feature space and achieve more accurate classification results.
Next modification of the basic model (1)-(2) controls a forward diffusion coefficient on the edges of the new observation points, see Fig 2 and 3. We can reduce the forward diffusion influence on the vertex by using the diffusion coefficient in the form
on all edges of , where is a parameter of the "diffusion neighborhood" size. In classification of the new observation, the aforementioned modification causes that only the points in a "-diffusion neighborhood", i.e. for which the diffusion coefficient is large than , are attracting new observation point . This modification is illustrated in Fig. 3.
2.2 Numerical discretization - natural network construction
Let us denote by any of the coordinates of . To discretize the equation (1), we use i) the balance of diffusion fluxes (inflows and outflows) in each vertex and ii) the approximation of the diffusion flux to the vertex along its edge .
First, let us define the diffusion flux approximation, which depends on the difference of values of function at the vertices and , as
for each edge , where represents the diffusion coefficient on the edge . If , it represents the diffusion inflow of the quantity into the vertex . On the other hand, if , it represents the diffusion outflow of the quantity from the vertex . Then the balance of diffusion fluxes in vertex is expressed by the equation
that means, the time derivative of in the vertex is positive - the value of increases in time, if the overall inflow to the vertex is greater than the overall outflow from the vertex . Vice versa, the time derivative of in the vertex is negative if the sum of inflows and outflows in the vertex is negative, i.e. outflows from the vertex are greater than inflows. When we substitute the approximation of the diffusion flux (6) into the balance equation (7), we obtain
The right hand side of the equation (8) in the graph theory represents the so-called "graph-Laplacian" (see equation (12) in calculusOnGraphs or equation (2.5) in edge-basedLaplacian ) which is given for a weighted complete undirected graph by relation
where represents a measure of the vertex and represents a "weight" of the edge in the weighted graph. In numerical mathematics, we would understand the "Laplacian" defined in this way as an averaged Laplace operator on a "finite volume" with a measure (area/volume) . Our numerical discretization (8) of the diffusion equation on the complete undirected graph corresponds to the choice , which is the standard choice for the vertex measure also in the graph theory, see calculusOnGraphs .
For the time discretization, we use the semi-implicit approach, see e.g. semiImpl . The time interval is divided uniformly into time steps and let denote the size of the time step. For the approximation of time derivative we use the finite difference method and obtain
Since the diffusion coefficient at the edge can depend on the unknown quantity , see (2)-(5), and thus can change over time, we take its value from the previous time step. We denote it by and obtain the semi-implicit scheme in the form
The semi-implicit scheme (11) can by rewritten in each time step into the system of linear equations
This system of equations is represented by a full matrix and as we have said before, for a complete undirected graph it is not necessary to define any boundary condition.
In the case of classification of the data from -dimensional feature space, our diffusing variables are the Euclidean coordinates of the vertices of the graph . In general, we get in each time step systems of linear equations
which are interconnected by the diffusion coefficient , which depends on all , , and can be written in the form
Let us denote the -th cluster by , , where is the number of clusters. Let us have the vertices and , where , . The value in the diffusion coefficient (14) is given by the following values
For a reader convenience in Fig. 4 we present dynamics of the network by (13) - (15) in order to show how the points inside the given clusters are moving together and clusters themselves are keeping away. In practice, in the learning phase and also in the application phase, we add new observation to the network and run the dynamics of the network. In the application phase we add to the network completely new observation while in the learning phase the new observation is taken away from the learning dataset. In both cases, the dynamics is modified in such a way that all other points are moving by (13) - (15) but for the new observation , , , the diffusion coefficient is set to
are given constants.
A stopping criterion is applied for dynamics of the network. We call it the histogram stopping criterion because it calculates the number of occurrences (frequency) of evolving points in prescribed spatial cells in every time step. That means the D grid with cells given by a specific spacing is created, see the grey grid in Fig. 6 for 2D case illustration. In numerical experiments presented below, we always use the spacing . We also determine which of the given clusters is the smallest one, and let the variable represent the number of points in this smallest cluster. Then, in every time step the histogram development is monitored, and whenever the number of points (denoted as the frequency in Fig. 6 left up corner) inside a grid cell is greater or equal to , the cell is marked, see the red cells in Fig. 6. At the same time, a specific neighborhood of every marked cell is examined, see the yellow subdomains in Fig. 6. The examined neigborhood is given by the interior of the concentric squares with Chebyshev radius and , respectively, in Fig. 6 we illustrate the situation where and , and such parameters are also used in the computations presented below. If there are only zero values of the frequency in all the cells inside the examined neighborhood of the marked cell, we claim that a cluster was formed in the marked cell. The dynamics of the network is stopped when the number of clusters formed is equal to the number of clusters given.
For quantifying the relevancy of classification of any new observation , , the relevancy coefficient is defined. We define it by using the information to which cluster the new observation is classified combined with the information about distances of the point representing the new observation and the centroids of all given clusters.
The new observation is classified into the cluster , if after stopping the network dynamics its closest point in the network is a point from the cluster and, at the same time, their distance is less than
. Otherwise, the point is not classified to any cluster, it is called the outlier and its relevancy will be equal to 0 in all clusters.
Now, let us assume that the new observation is classified into the cluster . The new observation has given its initial position and we calculate the centroids of the given clusters at the final time step by
where is the number of points in the cluster . Then we calculate the distance between the new observation and the centroid of the cluster to which it is assigned by the network dynamics,
and the average distance of the new observation to all other cluster centroids,
The above distances are used to define the quantity
which is in the range and which is the basis for definition of the relevancy coefficient. When the position of the new observation is close to the centroid of the cluster to which it is assigned then is close to 1, see Fig. 7 left. In this case, the relevancy of the resulting classification should be high. The quantity is close to or less if the distance of the new observation to the centroid of the cluster to which it is assigned is similar or greater than its distance to other clusters centroid, see Fig. 7 right. In this case, the resulting relevancy of the classification should be significantly reduced. We use the logistic function
which after linear rescaling from the interval to the interval give the final definition of the relevancy coefficient for any new observation ,
While is linearly decreasing from 1 to 0, depending on the distances ratio , the final relevancy coefficient has the nonlinear character, see Fig. 8. It sets the relevancy values to be close to 1 for all new observations belonging to a neighborhood (size of which depends on ) of the centroid of the cluster to which it is assigned. In all numerical experiments presented below, we use . The relevancy coefficient defined by (22) is used in the definition of relevancy maps which give a useful information on Sentinel-2 image pixels and image subareas membership in respective clusters, see Tabs. 3 - 6.
2.3 Ground-based vegetation data sampling - Application to habitat identification and prediction
One of the possible applications of the natural numerical networks described above is ecology and nature conservation. The identification and classification of the Natura 2000 habitats by using remote sensing data still has strong limitations due to habitats’ natural character and variability. In our application, we use a natural numerical network for the classification of the protected Natura 2000 forest habitats in the territory of Slovakia. There were four Natura 2000 forest habitats dominant in Western Slovakia chosen for the classification, as shown in Table 1. Therefore, we have four clusters , where , and . All habitat areas borders were semi-automatically segmented in NaturaSat software segSemiAuto ; Ambroz and checked in the field by botany experts during vegetation seasons of 2019 and 2020.
|Code||Name of habitat||Color|
|91E0||Alluvial forests with Alnus glutinosa and Fraxinus excelsior||red|
|91G0||Pannonic woods with Quercus petraea and Carpinus betulus||yellow|
|9110||Luzulo-Fagetum beech forests||purple|
There were 125 areas segmented in red, green, blue and near-infrared channels of Sentinel-2 data ESASentinel . We denote the segmented areas as , where , and . This input from field experts contains 30 segmented areas of 91E0 habitat, which consists of mixed ash-alder alluvial forests in temperate and boreal Europe (Alno-padion, Alnion incanae, and Salicion albae); 29 91F0 areas, which consist of riparian mixed forests of Quercus robur, Ulmus laevis and Ulmus minor, Fraxinus excelsio or Fraxinus angustifolia, along the great rivers of the Atlantic and Middle European provinces; 32 91G0 areas, which consist of Pannonic woods with Quercus petrea and Carpinus betulus; and 34 9110 areas, which consist of Luzulo-Fagetum beech forests habitats. The Sentinel-2 data from September 10, 2018 covering Western Slovakia were used, as shown in the left square in Fig. 9. More specifically, we studied the Natura 2000 habitats in the Podunajská nížina lowland along the Danube river (see Fig. 10 left up), the Záhorská nížina lowland along the Morava river (see Fig. 10 right up) and the habitats in the Malé Karpaty Mts. (see Fig. 10 bottom row).
2.4 Obtaining of multispectral data and habitat spectral characteristics
Sentinel-2 is a satellite of the European Space Agency (ESA) designed for the Copernicus European Union’s Earth observation program focused on the observation of the atmosphere, land, seas and climate on Earth Copernicus . All data acquired by the Sentinel-2 satellite are systematically processed, and only the products of this process are available for users ESASentinel . We use the Level-2A product that provides Bottom Of Atmosphere (BOA) reflectance images and offers 17 channels. In addition to these 17 channels, we calculate one more, the normalized difference vegetation index (NDVI) EOS
giving further useful information on habitat status. Thus, for the classification and the feature space construction, we use 18 channels; and for every channel, we compute the statistical characteristics - the mean, the standard deviation, the minimum value and the maximum value in a prescribed image subarea. Therefore, the feature space is the 72-dimensional Euclidean space; i.e., its dimensions are .
For classification purposes, first, in each segmented area , where , we created a square with a randomly chosen center in a pixel and Chebyshev radius . For large segmented areas, can be chosen; and for small areas, can be chosen smaller, as shown in Fig. 11. Such squares are used for building the learning datasets and also for the construction of the so-called relevancy maps defined below. The statistical characteristics of the above mentioned 18 channels are computed for every square , and they form the coordinates of points in the 72-dimensional feature space. The initial network graph is constructed such that every vertex of the graph is given by one such point corresponding to one square (or to one pixel or one segmented area , as we may say).
2.5 Relevancy map
The relevancy map is a grayscale image with the same size as the images from Sentinel-2 optical channels. After finding the optimal parameters of the natural network, the square is created in every image pixel (and not only inside the segmented areas as described above). For every , the statistical characteristics of the square are computed and considered as new observation . We note that to create required squares along the Sentinel-2 image boundary and compute the statistical characteristics there we use the reflection of corresponding values from the image interior to the exterior. The new observation is added to the graph as a new vertex. Every new observation is classified by the natural network and its relevancy coefficient is computed by (22). Finally, depending on the Chebyshev radius of the square , the relevancy map , , is defined for every cluster in every pixel as follows
2.6 Learning phase and network graph topology optimization
First, let us consider the graph having
vertices as described above and call it the (initial) learning dataset LDS125. All data from the learning dataset are labelled by the number of clusters to which they belong. As we stated above, we apply PCA to the original 72-dimensional points representing the data. PCA finds the coordinate system (basis) in which the highest variance (variability) of the data is in the first coordinate and it decreases subsequently with further coordinates. The change of the basis is represented by the linear transformation (matrix) that is applied to every point of the dataset, and we obtain coordinates of every point in the new coordinate system. Then, we are able to reduce the dimension of the feature space to, considering only the first two coordinates of every point, as shown in the top left of Fig. 4. As we observed experimentally, further coordinates do not help differentiate clusters and can be abandoned. After PCA matrix transformation, we also scale the coordinates of the points into the range which helps in the model parameter tuning process. In fact, for any data, we can then use the same ranges of the natural network model parameters depending on the point distances.
The main goals of the learning phase and the network graph optimization are to tune the parameters of the model (13)-(16) and optimize the structure of the graph itself to achieve the highest possible classification accuracy for observations from the learning dataset. To achieve that goal, we subsequently remove the cluster label from each vertex of graph , set it as the new observation and classify it using the model, as shown in Fig. 5. We vary the model parameters , and and choose the combination of parameters that results in the greatest number () of correctly classified observations from the learning dataset. If is the number of all observations, our goal is to achieve a success rate close to 1.
In the tuning process, the range for the diffusion coefficient parameters is with the step size , where ; and the range for parameter is the interval with the step . As we tune three parameters only, we are able to go through all discrete parameters combinations for every new observation (a type of brute force approach) and find the combination that gives the highest possible classification accuracy . For any parameter combination, we apply at most time steps of the natural network dynamics. We have chosen , but in most cases, the histogram stopping criterion is fulfilled for fewer time steps and classification is fast. For the numerical solution of the linear systems of equations (13), we use the SOR (Successive Overrelaxation) iterative method. The first row of Table 2 shows the results. We achieved the best success rate by using the model parameters , and .
The achieved the best success rate of 0.84 should be improved, so we performed further steps in the learning phase. It is quite clear that the initial random choice of the representative squares inside the segmented areas , where , is not the optimal approach. The statistical characteristics of those squares do not necessarily correspond optimally to the statistical optical characteristics of the corresponding habitat in the Sentinel-2 image data. A strategy to solve this problem is based on spatial adjustment (shift) of representative squares inside the segmented areas such that we obtain a higher classification success rate with new (updated) vertices of the graph . We can understand this step as a learning dataset adjustment. Since the representative squares have various radii (depending on the size of the segmented area), we construct the relevancy maps , as shown in section 2.5, for . Since the relevancy map gives the relevancy of classification for every image pixel, we will be able to compare the relevancies of the pixels inside the segmented areas and choose the pixel with the highest relevancy. The relevancy maps are constructed by using the graph corresponding to the initial learning dataset LDS125 and by using the optimal network parameters , , and . Then, we check whether there is an such that we are able to find a new pixel in which (ideally close to 1), where is the relevancy map of the cluster to which the segmented area belongs. If we find such a pixel in the relevancy map , the new square is constructed. Every for which it was possible to find the new representative square, with higher , is replaced by the new one and the adjusted learning dataset LDS125adj is created. The network parameter tuning process is conduct again, but it now uses the LDS125adj learning dataset. The result is shown in the second row of Tab. 2; and the best success rate was , which is higher than the previous one. It was achieved by using the parameters , , and . We can conclude that the adjustment of representative squares significantly increased the classification success rate.
However, in the previous adjustment step, it was not always possible to find an adjusted representative square for all , where . This was caused by the fact that for some segmented areas , only zero values of relevancy were computed in all interior pixels in all relevancy maps , where . We checked whether these only zero values in all relevancy maps were changed in those segmented areas when using LDS125adj dataset with its optimal parameters. Thus, we again constructed the relevancy maps. Checking the results, we found that there are still seven areas for which no inner pixel can be found such that it has relevancy greater than zero. It is clear that those areas cannot contribute in any way to increase the classification success rate. We call the area fulfilling the condition
unclassifiable and we remove all squares corresponding to the unclassifiable areas from the learning dataset LDS125adj and create the final learning dataset LDS118 containing only labeled observations. Using this adjustment, we changed the topology of the graph itself, and we call this step the network graph topology optimization. The model parameters were tuned again and the best success rate of classification was , as shown in the third row of Tab. 2, for the optimal parameters , , and . This success rate is high and allow us to use such optimally tuned (trained) natural network in practical applications presented in the next subsections.
3 Results and discussion
3.1 Relevancy maps for Western Slovakia
In the Methods section, we discussed the learning process that led to the successfully optimized and trained natural network. Finally, the graph of the trained natural network contains vertices and the optimal parameters are , , and . This subsection is devoted to the construction of the final relevancy map for each explored Natura 2000 habitat. It was clear from the learning phase that the radius of the chosen square areas should vary between 3 (for small segmented areas) and 5 (for large segmented areas) to get the optimal results. For this reason, we compute three relevancy maps for the squares with radii , and , as shown in Fig. 12. The final relevancy map , where , is obtained by taking the maximum of those three relevancies in every pixel , i.e., we define it by
The qualitative (visual) comparisons of the Natura 2000 habitat segmented areas with the final relevancy maps for habitats 91E0 and 91F0 are plotted in Fig. 13. As the figure shows, the interior of the segmented areas on each relevancy map contains bright colors. This reflects the correct high relevancy of classification inside the segmented areas and thus the correct assignment of the image pixels to the Natura 2000 habitat.
For the quantitative evaluation of the classification accuracy, we calculate the mean relevancy inside each segmented area. It is clear that the relevancy values inside the inner narrow band of width are irrelevant because the squares , for from the narrow band, also partially cover the pixels outside the area, cf. Fig. 12. Thus, we shrink the boundary curve by the distance , and the mean value of the relevancy is computed inside the smaller curve. In Tables 3 - 6, we show the mean relevancy inside the segmented areas of the 91E0, 91F0, 91G0 and 9110 Natura 2000 habitats in each of the final relevancy maps.
The studied Natura 2000 habitat segmented areas obtain nonzero values in not only one final relevancy map. This occurs because the final relevancy map is constructed by the combination of the relevancy maps with different s and by the fact that pixels inside the segmented areas can be classified differently by the natural network. We notice that most segmented areas obtain the highest mean relevancy for the habitat to which they belong. This is a feature supporting the usability of the final relevancy maps in the classification of new areas and finding new appearances of habitats. There are also some cases where the segmented area has the highest mean relevancy in a different cluster than that to which it should belong. This is, however, in a large majority of cases explainable. E.g., it may be caused by a similar species composition of habitats or by the very often difficult decision of experts in the field when classifying the habitat. Thus, we can conclude that the final relevancy maps give very useful information for an expert decision on habitat classification. This is further discussed below.
The 91E0 softwood floodplain forests are the best detectable habitat among the studied habitats, as shown in Tab. 3, where all areas have the highest mean relevancy in the correct 91E0 cluster. The only disputable result is the habitat area 91E0_Sihot, which was classified correctly as 91E0 but the relevancy values are high for both 91E0 and 91F0 habitats. It can be caused by the high cover of Ulmus minor, the hardwood tree species at this locality, which is also typical for the 91F0 habitat, but in combination with different trees. The floodplain forests at the Sihoť Island are unique due to the specific co-occurrence of softwood floodplain forest species Populus sp. and Salix sp. together with Ulmus minor; therefore, the obtained relevancy values in Tab. 3 are meaningful. It seems that when hardwood tree species are mixed in the 91E0 forest habitat, the mean relevancy of such an area also has a higher value for the 91F0 habitat cluster. This fact indicates recognition of transitional floodplain forest type that is common in natural conditions and also partly occurs in the localities in Číčov (Cicov_1, Cicov_5) and Medveďov (Medvedov_2).
Within the 91F0 habitat, three areas are misclassified, and some others have close relevancies in both 91F0 and 91E0 habitats, as shown in Tab. 4. The segmented areas 91F0_Kopac1, 91F0_Vysoka_pri_Morave_6 and 91F0_Vysoka_pri_Morave_8 were classified contrary to expert opinions. These forests represent a transitional type between 91E0 and 91F0 habitats, and experts classified them as 91F0. In many cases, it is hard to classify transitional floodplain forests fully objectively. The categorization is dependent on the subjective decision of the field expert. The relevancy shows possible classification within the 91E0 habitat of these three areas, which should be discussed among botany experts. In the alluvial country with natural characteristics, the softwood (91E0) and the hardwood (91F0) floodplain forests are often interconnected. They form transitional zones with hardwood tree species in the softwood habitat and, conversely, with softwood species in the hardwood floodplain forest. Our results show that the relevancy maps are sensitive to these transitional forests and can give information about such mixed composition of the floodplain forest and possible classification within the next-standing habitat.
A few misclassified areas are found in the 91G0 habitat, as shown in Tab. 5. They are misclassified mainly within 91F0 instead of the correct 91G0. In addition, via a more in-depth look, we see that all misclassified areas contain Quercus petraea species together with Carpinus betulus. For the 91F0 habitat, Quercus robur is typical, but almost no difference between the two Quercus at the species level is identified. However, since the 91F0 habitat occurs in river alluvia while the 91G0 habitat occupies a lower hill, adding the information from a digital elevation model to the feature space may solve this problem.
In Table 6
, we present the results for the 9110 habitat areas. In some cases, they form transitional zones with 91G0 habitat because they share some tree species, Carpinus betulus, Tilia cordata, Prunus avium, Acer platanoides, etc. The beech forests in the Carpathians’ upper parts are monodominant and thus are very easily detectable while forests in lower altitudes host more 91G0 species. Thus, the information from the digital elevation model can again improve the classification. One of the misclassified areas (9110_Raca_2) represents a very old-growth forest with a lower canopy cover, and probably some undergrowth species with different optical characteristics influenced the result.
|Area code||91E0 habitat||91F0 habitat||91G0 habitat||9110 habitat|
|Area code||91E0 habitat||91F0 habitat||91G0 habitat||9110 habitat|
|Area code||91E0 habitat||91F0 habitat||91G0 habitat||9110 habitat|
|Area code||91E0 habitat||91F0 habitat||91G0 habitat||9110 habitat|
The natural vegetation has a continuous character, and Natura 2000 habitats stated in the habitats directive describe the essential, most characteristic forms of habitats. In real nature, these types are often interconnected and hard to classify, even by experts in the field. The presented results show the high classification success; moreover, it gives us information about the admixture of some tree species typical for next-standing habitats. This finding opens new opportunities in phytosociological and ecological research, especially in the fields of ecotone zones and vegetation gradients.
3.2 Exploration of other regions in Slovakia and finding new appearances of protected habitats
After obtaining the successfully trained natural network, we can explore further areas of the occurrence of protected habitats in territories neighboring Western Slovakia, as shown in Fig. 9. We explored the areas of Central and Eastern Slovakia and obtained very promising results. By using the relevancy maps, we confirmed the appearance of the 91F0 habitat around the Latorica River in the south of Eastern Slovakia and found the new appearance of 91E0 protected habitat along the Rimava River near the village of Dubovec in the southern part of Central Slovakia. The main goal of the relevancy map construction, to find the new appearances of protected habitats in an automated way, was thus achieved.
During the Natura 2000 mapping campaigns, the 91F0 floodplain forests around the Latorica River were sampled using vegetation plots. We use the trained natural network to create the final relevancy maps for the Eastern Slovakia region, and the details of the Sentinel-2 image and the 91F0 final relevancy map can be seen in Fig. 14. The figure shows 11 points (yellow color) representing the vegetation plots - phytosociological relevés, which mark the appearance of the habitat. Usually, they are not far from the boundary of the habitat. First, we see the correspondence of these points and bright colors in the relevancy map of the 91F0 habitat, which indicates the correct classification of the pixels along the Latorica River. Furthermore, we automatically segmented the habitat area segAuto starting the segmentation from 11 habitat marking points. The evolution of the segmentation curves undergoing topological changes is presented in Fig. 14 with the final segmentation result plotted in the bottom right picture both on the Sentinel-2 image and relevancy map. The automatic segmentation identified the compact habitat area and also included a thin river branch, a forest road and small areas with younger floodplain forests, which is accepted because these objects are common parts of the habitat. Furthermore, the relevancy map was able to detect even these small parts, and high relevancy values were assigned only to the best representative area of the habitat. We computed also the mean value of the relevancy in the 91F0 relevancy map inside the segmented area; and it is , which indicates the correct classification of the area along the Latorica River.
During the exploration of the relevancy maps created for the southern area of the Central Slovakia region, we found that the natural network assigns high 91E0 relevancy to the area around the Rimava River close to the village of Dubovec (Fig. 15 right). The discovered area has never been identified as a target habitat, and no databases contain such information (the Slovak vegetation database or the database of the State Nature Conservancy of the Slovak Republic where all currently known areas of Natura 2000 habitats are collected). We segmented the area by applying the automatic segmentation method segAuto , as shown at the bottom right in Fig. 15, to obtain the final segmentation result. The computed mean relevancy inside the segmented area is equal to , and it again indicates a possible new appearance of the 91E0 habitat. Finally, the botanists went into the field and confirmed that the newly discovered area is classified as 91E0 Natura 2000 habitat, as shown in photos of the area in Fig. 16.
3.3 Exploration of the alluvial forests along the Danube River in the Central and South Europe
During the field exploration of the Danube River floodplains, the appearance and positions of the habitats were sampled using vegetation plots. We use the natural network to explore the alluvial forests along the Danube River and we made a qualitative and quantitative comparison of the vegetation plots and the relevancy map.
The floodplain forests occur on the banks of the Danube River from Central to South Europe. A large area of these forests is situated in Upper Austria around Linz city, where softwood floodplain forests (91E0 habitat) form a mosaic with monodominant plantations of Canadian poplars or maples. Vegetation plots were sampled in this area, and they were used to verify the relevancy map. All of the vegetation plots fit inside the regions identified by the relevancy map. Fig. 17 depicts two of them, situated inside the areas, which the natural network denotes as 91E0 habitat areas. Automatic segmentation was applied to these areas, and the result is also depicted in Fig. 17. The final curves surround the areas with high 91E0 relevancy. The 91E0 mean relevancy of the areas was equal to for the smaller area and for the larger area.
The relevancy map was created, and vegetation plots were also sampled in the Danube River alluvia in Vojvodina (North Serbia). The identification of 91E0 habitat areas was slightly less successful in this area. Two plots were not recognized due to the presence of the moistest forest types dominated by sparse willows and absenting poplars that are not so widespread in the upper parts of the Danube alluvia where the natural network was trained. One habitat area was not found because the plot occurs along a very thin forest line, which is explainable by the Sentinel-2’s resolution and the mechanism of creating the relevancy map. All other plots were identified correctly, and two examples of habitat areas in the relevancy map are shown in Fig. 18 and Fig. 19, respectively. The first example is situated on the left bank of the Danube River near the village of Vajska. This vegetation plot was used as a starting point for the automatic segmentation algorithm, as shown in Fig. 18. This thin habitat area was segmented with a 91E0 mean relevancy equal to .
A complex habitat system from the Karadordevo nature reservation, identified in the relevancy map and segmented by automatic segmentation, is presented in Fig. 19. The starting point for the segmentation was the vegetation plot located on the north border of this area. The entire area of the 91E0 habitat was successfully segmented while the mosaic of wet meadows and young clear-cuts inside the habitat area was excluded. The mean relevancy of the segmented area is , which means correct classification within the 91E0 habitat.
Another studied area was the delta of the Danube River in Romania. It is difficult to explore the region due to the complicated accessibility of most of the area. The character of the landscape means that localities can only be reached by small boats; thus, the importance of creating relevancy maps of such areas is significant. Fig. 20 plots the Saint George Branch of the Danube River, and the relevancy map corresponds very well with the vegetation sampled by plots. In the detail of the Danube delta’s branch relevancy map, we can clearly see bright-colored 91E0 habitat and black-colored surrounding swamps and other ecosystems, which indicates correct classification.
The presented results give new perspectives for diverse science disciplines - mathematical modelling, vegetation science and ecology, remote sensing, nature conservation and mapping, and subsequent delivery of ecosystem services. The identification of plant communities on the scale of, e. g., Natura 2000 habitats, using remote sensing has been an open challenge for field ecologists over the last decades. Rapidly developing remote sensing techniques and different data-mining approaches were implemented to monitor land surface types such as agricultural land, water bodies, abandonment land, natural and plantation forests of different types, meadows with various management practices, built-up areas, etc, see e.g. Franklin_etal.2001 ; Fagan_etal.2015 ; Laurin_etal.2016 ; Noi ; Erinjery_etal.2018 ; Cheng&Wang2019 ; Wasnievski_etal.2020 . However, due to the complicated character of target nature phenomena, it was not possible to reach the detailed scale of Natura 2000 habitats by using existing methodologies and satellite data. Thus, we developed the novel method - the natural numerical network - for that purpose. We introduced its definition in the form of the discretized forward-backward diffusion equation and applied it to the classification of Natura 2000 forest habitats. We also introduced the new concept of relevancy maps, which allow the automatic recognition of Natura 2000 habitat areas in remote sensing data provided by Sentinel-2 optical information and finding new appearances of protected habitats via the satellite data. We are not aware of any other method that would identify and explore the Natura 2000 habitats or similarly detailed plant communities on such an exact level and accuracy by using the remote sensing data.
This work was supported by the grants APVV-16-0431, APVV-19-0460 and ESA contracts 4000122575/17/NL/SC and 4000133101/20/NL/SC.
- (1) S. E. Dreyfus, Artificial neural networks, back propagation, and the Kelley-Bryson gradient procedure, Journal of Guidance Control Dynamics 13 (5) (1990) 926–928. doi:10.2514/3.25422.
- (2) B. C. Csaji, H. Eikelder, Approximation with artificial neural networks, Master’s thesis, Eötvös Lorand University (2001).
- (3) I. Goodfellow, Y. Bengio, A. Courville, Deep Learning, MIT Press, 2016, http://www.deeplearningbook.org.
- (4) E. Haber, L. Ruthotto, Stable architectures for deep neural networks, Inverse Problems 34 (1) (2018). doi:10.1088/1361-6420/aa9a90.
- (5) K. He, X. Zhang, S. Ren, J. Sun, Deep residual learning for image recognition, 2016, pp. 770–778. doi:10.1109/CVPR.2016.90.
B. Chang, L. Meng, E. Haber, L. Ruthotto, D. Begert, E. Holtham,
architectures for arbitrarily deep residual neural networks
, in: 32nd AAAI Conference on Artificial Intelligence, AAAI 2018, Vol. 32, 2018, pp. 2811–2818.
- (7) P. Liu, A survey of remote-sensing big data, Frontiers in Environmental Science 3 (2015) 45. doi:10.3389/fenvs.2015.00045.
- (8) C. Randin, M. Ashcroft, J. Bolliger, J. Cavender-Bares, N. Coops, S. Dullinger, T. Dirnböck, S. Eckert, E. Ellis, N. Fernández, G. Giuliani, A. Guisan, W. Jetz, S. Joost, D. N. Karger, J. Lembrechts, J. Lenoir, M. Luoto, X. Morin, D. Payne, Monitoring biodiversity in the anthropocene using remote sensing in species distribution models, Remote Sensing of Environment 239 (2020) 111626. doi:10.1016/j.rse.2019.111626.
- (9) C. Corbane, S. Lang, K. Pipkins, S. Alleaume, M. Deshayes, V. Millán, T. Strasser, J. Vanden Borre, T. Spanhove, M. Förster, Remote sensing for mapping natural habitats and their conservation status - new opportunities and challenges, Intenational Journal of Applied Earth Observation and Geoinformation 37 (2015) 7–16.
- (10) A. Lausch, M. Heurich, P. Magdon, D. Rocchini, S. Karsten, J. Bumberger, D. King, A Range of Earth Observation Techniques for Assessing Plant Diversity, Remote Sensing of Plant Biodiversity, 2020, pp. 309–348. doi:10.1007/978-3-030-33157-3_13.
- (11) H. Ullerud, A. Bryn, R. Halvorsen, L. Hemsing, Consistency in land-cover mapping: Influence of field workers, spatial scale and classification system, Applied Vegetation Science 21 (2) (2018) 278–288. doi:10.1111/avsc.12368.
- (12) M. Chi, A. Plaza, J. Benediktsson, Z. Sun, J. Shen, Y. Zhu, Big data for remote sensing: Challenges and opportunities, Proceedings of the IEEE 104 (11) (2016) 2207–2219. doi:10.1109/JPROC.2016.2598228.
D. Rocchini, V. Petras, A. Petrasova, N. Horning, L. Furtkevicova, M. Neteler, B. Leutner, M. Wegmann, Open data and open source for remote sensing training in ecology, Ecological Informatics 40 (2017) 57–61.doi:10.1016/j.ecoinf.2017.05.004.
- (14) European Environmental Agency, The natura 2000 protected areas network, https://www.eea.europa.eu/themes/biodiversity/natura-2000/the-natura-2000-protected-areas-network (2020).
- (15) State nature conservation SR, Natura 2000, http://www.sopsr.sk/natura/index1.php?p=3&lang=sk (2020).
- (16) M. Bock, P. Xofis, J. Mitchley, G. Rossner, M. Wissen, Object-oriented methods for habitat mapping at multiple scales – case studies from northern germany and wye downs, uk, Journal for Nature Conservation 13 (2-3) (2005) 75–89. doi:10.1016/j.jnc.2004.12.002.
- (17) J. Vanden Borre, D. Paelinckx, S. Mucher, L. Kooistra, B. Haest, G. De Blust, A. M. Schmidt, Integrating remote sensing in natura 2000 habitat monitoring: Prospects on the way forward, Journal for Nature Conservation 19 (2) (2011) 116–125. doi:10.1016/j.jnc.2010.07.003.
- (18) European Space Agency, Sentinel 2, https://sentinel.esa.int/web/sentinel/missions/sentinel-2/data-products (2020).
- (19) A. Bondy, U. S. R. Murty, Graph Theory, 1st Edition, Springer-Verlag London, 2008.
- (20) J. Friedman, T. Jean-Pierre, Calculus on graphs, CoRR cs.DM/0408028 (2004).
- (21) P. Perona, J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE Transactions on Pattern Analysis and Machine Intelligence 12 (7) (1990) 629–639. doi:10.1109/34.56205.
- (22) J. Friedman, J.-P. Tillich, Wave equations for graphs and the edge-based laplacian, Pacific Journal of Mathematics 216 (2) (2004) 229–266. doi:10.2140/pjm.2004.216.229.
- (23) K. Mikula, N. Ramarosy, Semi-implicit finite volume scheme for solving nonlinear diffusion equations in image processing, Numerische Mathematik 89 (3) (2001) 561–590. doi:10.1007/PL00005479.
- (24) K. Mikula, J. Urbán, M. Kollár, M. Ambroz, I. Jarolímek, J. Šibik, M. Šibíková, Semi-automatic segmentation of natura 2000 habitats in sentinel-2 satellite images by evolving open curves, Discrete and Continuous Dynamical Systems - Series S 14 (3) (2021) 1033–1046. doi:10.3934/dcdss.2020231.
- (25) M. Ambroz, M. Kollar, K. Mikula, Semi-implicit scheme for semi-automatic segmentation in naturasat software, in: ALGORITMY 2020 - 21st Conference on Scientific Computing, Vysoké Tatry - Podbanské, Slovakia, September 10 - 15, 2020. Proceedings of contributed papers, Vydavateľstvo SPEKTRUM STU, 2020, ISBN: 978-80-227-5032-5, 2020, pp. 171–180.
- (26) Copernicus, Copernicus: Europe’s eyes on earth, https://www.copernicus.eu/en/about-copernicus/copernicus-detail (2020).
- (27) Earth Observing System, Normalized difference vegetation index, https://eos.com/ndvi/ (2020).
- (28) I. T. Jolliffe, Principal Component Analysis, 2nd Edition, Springer-Verlag New York, 2002. doi:10.1007/b98835.
A. Rencher, Methods of Multivariate Analysis, Wiley-Interscience, 2002.doi:10.1002/0471271357.
- (30) K. Mikula, J. Urbán, M. Kollár, M. Ambroz, I. Jarolímek, J. Šibik, M. Šibíková, An automated segmentation of natura 2000 habitats from sentinel-2 optical data, Discrete and Continuous Dynamical Systems - Series S 14 (3) (2021) 1017–1032. doi:10.3934/dcdss.2020348.
- (31) S. E. Franklin, Remote Sensing for Sustainable Forest Management, 1st Edition, Lewis Publishers/CRC Press, 2001.
- (32) M. Fagan, R. Defries, S. Sesnie, P. Arroyo, C. Soto-Castro, A. Singh, P. Townsend, R. Chazdon, Mapping species composition of forests and tree plantations in northeastern costa rica with an integration of hyperspectral and multitemporal landsat imagery, Remote Sensing 7 (5) (2015) 5660–5696. doi:10.3390/rs70505660.
- (33) G. Vaglio Laurin, N. Puletti, W. Hawthorne, V. Liesenberg, P. Corona, D. Papale, Q. Chen, R. Valentini, Discrimination of tropical forest types, dominant species, and mapping of functional guilds by hyperspectral and simulated multispectral sentinel-2 data, Remote Sensing of Environment 176 (2016) 163–176. doi:10.1016/j.rse.2016.01.017.
P.-T. Noi, M. Kappas, Comparison of random forest, k-nearest neighbor, and support vector machine classifiers for land cover classification using sentinel-2 imagery, Sensors (Switzerland) 18 (1) (2018) 18.doi:10.3390/s18010018.
- (35) J. J. Erinjery, M. Singh, R. Kent, Mapping and assessment of vegetation types in the tropical rainforests of the western ghats using multispectral sentinel-2 and sar sentinel-1 satellite imagery, Remote Sensing of Environment 216 (2018) 345–354. doi:10.1016/j.rse.2018.07.006.
- (36) K. Cheng, J. Wang, Forest type classification based on integrated spectral-spatial-temporal features and random forest algorithm—a case study in the qinling mountains, Forests 10 (7) (2019) 559. doi:10.3390/f10070559.
- (37) A. Waśniewski, A. Hoscilo, B. Zagajewski, D. Mouketou-Tarazewicz, Assessment of sentinel-2 satellite images and random forest classifier for rainforest mapping in gabon, Forests 11 (9) (2020) 941. doi:10.3390/f11090941.