Introduction

Today, due to modern satellites, it is possible to obtain high-quality images of the Earth's surface. High spatial resolution images are of special practical value [1]. These remote sensing data are used for monitoring the Earth's surface, mapping, tracking deforestation, assessing the consequences of natural disasters, and in many other areas of human activity [2]. Such modern satellites as the WorldView-2 and the WorldView-3 are capable of producing images with a spatial resolution of up to 1.24 m per pixel in the panchromatic channel and up to 0.3 m per pixel in the multispectral channels. The WorldView-3 satellite is capable of covering up to 680,000 square kilometers of the Earth's surface per day, and WorldView-2 up to one million square kilometers. The huge number of multi-channel satellite images arriving to the Earth every day requires fast, high-quality processing to obtain useful information in a timely manner [3]. There is a need to develop methods of automated computer processing of data. Satellite images are actively used in agriculture. With the help of images obtained in the near-infrared parts of the spectrum, it is possible to assess the state and stages of crop growth, to determine the type of vegetation. Recently, aerial images have also been actively used in agriculture.

One of the urgent tasks of remote sensing data processing is the counting of trees in plantations. Cultivation of fruit trees, almonds, walnuts and hazelnuts, oil palms is an important part of the agricultural industry and a significant source of income in such countries as the USA, Malaysia, Brazil, etc [4]. In order to effectively organize work and forecast the harvest, it is necessary to know the number of trees in the cultivated area. Until now, human labor was used to count and assess the condition of individual trees on plantations, which was very time-consuming and time-consuming. Orchards and tree plantations can cover huge areas and have thousands of trees, so the task of automating harvest counting is very important. Thanks to satellite images, it is possible to create an effective automatic system of counting trees in plantations [5].

To solve this problem, it is necessary to choose effective methods of processing and analyzing data obtained by means of remote sensing. The task of identifying and counting trees can be divided into two subtasks: 1) identification of areas of the Earth's surface where there are agricultural plantations, in particular fruit tree plantations; 2) identification of individual objects in the selected area of the image for the purpose of counting them. In this work, these problems are solved using the methods of texture analysis of digital images based on the adjacency matrix of gray levels, image segmentation based on textural features, and selection of individual objects based on their spectral features.

Related works

Attempts have been made repeatedly to create an automatic tree identification and counting system [6]. Remote sensing technology is extensively utilized for monitoring and quantifying canopy growth, detecting potential plant diseases, and tracking changes within forest structures. Timely analysis of this data is critical for optimizing yields and evaluating the response of forests to climate anomalies [7,8].

The cultivation of oil palm is a vital contributor to agricultural productivity in numerous developing countries across the tropics. As such, conducting research to investigate and accurately quantify oil palm cultivation is both valuable and meaningful. In papers, CNN and RCNN were applied for tree detection and counting [9,10,11].

In paper [12], an algorithm for searching for trees in the image using a geometric-optical tree model was proposed, assuming that the center of the tree is at the point of maximum similarity between the model and the sample. The dome-shaped crown of the tree was taken as the basis of the geometricoptical model of the tree. Information about the lighting in the image was used to create it. The sun height, azimuth, and tree width parameters were determined automatically for each individual image, but the result of automatic parameter determination was no better than the user's visual assessment, and needs improvement.

LiDAR data were used in [13,14,15]. Using the LiDAR-based canopy height model and segmentation methods the tree crowns were delineated.

In paper [16], a method of counting fruit trees was also proposed. To find the crowns of individual trees, the LoG filter was used, which gives the strongest response to round objects. Since the LoG filter detects all objects that have a shape similar to the crown of a tree, the authors use vegetation indices and take into account only those objects that are green vegetation. Binarization of the red channel of the image and calculation of the NDVI vegetation index were used to find the vegetation in the image.

Vegetation indices were used in paper [17] to identify individual trees in a plantation. The most common vegetation indices were investigated, and those that give the maximum difference between the spectral characteristics of trees and the background were selected. The borders of the plantations were outlined manually on the pictures, the crowns of the trees were located as local maxima of the values of the vegetation indices.

In [18], an algorithm for searching for young palm trees in a satellite image using Haar features was proposed. Seven Haar features describing the shape of the tree were taken. In order to avoid incorrect identification of objects that are not trees, the obtained results were classified by the support vectors (SVM) method.

Materials and methods

Satellite images obtained by the WorldView-2 and the WorldView-3 satellites were used in this work. The WorldView-2 is the first commercial instrument with an eight-channel high-resolution spectrometer that includes traditional spectral channels as well as four additional ones. They provide higher accuracy in detailed analysis of the state of vegetation, selection of objects, analysis of the coastline and coastal water area. In its characteristics, the WorldView-2 meets the highest requirements. The data received from this satellite have a root mean square error (RMS) of no worse than 4 m without ground reference points [19].

The WorldView-3 satellite is designed for shooting in panchromatic and multispectral modes. Camera equipment is completely similar to that installed on the WorldView-2 satellite. The accuracy of the geo-positioning in the plan is 6.5 m or 4 m SKP without additional correction of the plan coordinates by ground reference points. The WorldView-3 conducts shooting in the following modes: VNIR ( Visible and Near Infrared )multispectral visible and near-infrared range, only 8 channels; SWIR ( Shortwave Infrared )mid-infrared range, allows shooting through haze, fog, smog, dust, smoke and clouds, only 8 channels; CAVIS ( clouds , aerosols , vapors , ice , snow )allows for atmospheric correction of only 12 channels with a spatial resolution of 30 m at nadir and wavelengths from 0.4 μm to 2.2 μm .

One of the methods that ensure the selection of useful information from remote sensing data is image segmentation. Segmentation allows you to identify homogeneous areas on the Earth's surface that correspond to certain natural objects [20].

Analysis of a number of satellite images and aerial photographs has shown that each image contains several types of land cover (forest, grass, soil, etc.), which can be characterized by their spectral textural properties. Fig. 1 shows a fragment of a satellite image. Texture fragments of vegetation can be divided into two main types: visually different in terms of spectral properties and structure; visually close in terms of spectral properties and structure. Such fragments belong to the same class, for example "forest", "soil", "grass", etc. Texture fragments within the same class are close in their characteristics [21].

Figure 1: Satellite image

The complex structure of satellite and aerospace images does not allow us to solve these problems using only the spectral properties of the images. The result of segmentation of the satellite image (Fig. 2) based on spectral characteristics showed that it is impossible to distinguish areas with a uniform texture using this method. Spectral properties of objects on the Earth's surface do not always provide complete information about the objects, as they depend on many factors, such as relief, soil type, climate, geographical location of the area. Additional a priori information such as image acquisition geometry and image context information must be used to improve the reliability of feature class decisions. To identify fruit tree plantations, textural properties of images were used in this work.

Calculation of the texture properties of an image based on the gray-level co-occurrence matrix

Texture analysis methods are widely used in the segmentation, analysis and deciphering of remote sensing imagery [22]. Various approaches to the detection and description of textures are known: statistical, geometric, structural, spectral, and model [23]. Texture analysis methods based on a onedimensional frequency histogram do not take into account the relative position of image pixels. They allow one to take into account only group properties, pixels belonging to one object on an aerial image.

The adjacency matrix of gray levels, or gray-level co-occurrence matrix (GLCM) [24] allows one to take into account the relative position of the pixels of the image and thus analyze the texture with pronounced spatial regularity. GLCM is calculated in one of the image channels and has dimensions L×L, where L is the number of gray levels in the image channel. It shows how often pixels with value i border pixels with values j in the horizontal (0 o ), vertical (90 o ), or diagonal (45 o and 135 o ) direction. We denote the GLCM as P:

𝑃 𝑟,𝜃 (𝑖, 𝑗) = |{(𝑘, 𝑠), (𝑡, 𝑣): 𝐼(𝑘, 𝑠) = 𝑖, 𝐼(𝑡, 𝑣) = 𝑗}|, (1) where i,j are matrix brightness levels 𝑃(𝑖, 𝑗 = 1, 𝐿); 𝐼(𝑘, 𝑠)and 𝐼(𝑡, 𝑣) are the values of the elements of the matrix P with the coordinates (𝑘, 𝑠)and (𝑡, 𝑣); r is the distance between the elements 𝐼(𝑘, 𝑠)and 𝐼(𝑡, 𝑣); θ is the angle between the elements 𝐼(𝑘, 𝑠)and 𝐼(𝑡, 𝑣)relative to the horizontal axis.

Based on the calculated adjacency matrices, the following indicators of textural features are calculated:

Average value:

𝜇 𝑖 = 𝜇 𝑗 = ∑ [𝑖 ∑ 𝑃(𝑖, 𝑗) 𝐿−1 𝑗=0 ] . 𝐿−1 𝑖=0(2)

Energy:

𝑇 𝑒 = ∑ ∑[𝑃(𝑖, 𝑗)] 2 𝐿−1 𝑗=0 . 𝐿−1 𝑖=0(3)

Variation:

𝜎 𝑖 2 = ∑ [(𝑖 − 𝜇 2 ) 2 ∑ 𝑃(𝑖, 𝑗) 𝐿−1 𝑗=0 ] . 𝐿−1 𝑖=0(4)

Homogeneity:

𝑇 ℎ = ∑ ∑ 𝑃(𝑖, 𝑗) 𝐿−1 𝑗=0 (1 + |𝑖 − 𝑗|) ⁄ ; 𝐿−1 𝑖=0(5)

where 𝑃(𝑖, 𝑗)is the frequency of occurence in the window of two pixels with brightness i and j at an angle α at a distance d; σ 2 is the root mean square deviation of pixel values in the window. By utilizing statistical points 1-4, it is possible to generate texture features that consider the relative positions of neighboring pixels within a given window. As a result, these features are particularly effective at describing textures that exhibit significant spatial regularity. The following assessments also apply:

The second angular moment

𝑇 2 = ∑ ∑(𝑃 (𝑖, 𝑗) 𝑀 ⁄ ) 2 , 𝐿 𝑗=1 𝐿 𝑖=1 (6)

where M is the total number of pairs of adjacent elements, this number is a measure of image homogeneity. For d =1, α=0, M =2 L y ( Lx -1). Contrast

𝑇 𝐶 = ∑ 𝑛 2 𝐿−1 𝑛 [∑ ∑(𝑃 (𝑖, 𝑗) 𝑀 ⁄ ) 𝐿 𝑗=1 𝐿 𝑖=1 ] , |𝑖 − 𝑗| = 𝑛,(7)

is determined by the magnitude of local variations of pixel values: the larger it is, the higher the contrast is.

Correlation coefficient

𝑇 𝑐𝑐 = 𝜎 𝑥 −1 𝜎 𝑦 −1 ∑ ∑ [𝑖𝑗 ( 𝑃(𝑖, 𝑗) 𝑀 ) − 𝑚 𝑥 𝑚 𝑦 ] , 𝐿 𝑗=1 𝐿 𝑖=1(8)

where m x, m y, σ x, σ y are mean values and root mean square deviations for 𝑝 𝑥 (𝑖) = ∑ 𝑃 (𝑖, 𝑗) 𝑀 ⁄ 𝐿 𝑗=1

and 𝑝 𝑦 (𝑗) = ∑ 𝑃 (𝑖, 𝑗) 𝑀 ⁄ 𝐿 𝑖=1

, respectively. The correlation coefficient is a measure of the linearity of the regression dependence of brightness in the image.

Variance

𝑇 𝐷 = ∑ ∑(𝑖 − 𝑚) 2 (𝑃 (𝑖, 𝑗) 𝑀 ⁄ ) 𝐿 𝑗=1 , 𝐿 𝑖=1(9)

determines the brightness variations from the average value.

characterizes the uneven distribution of brightness of image elements. The above characteristics were calculated for test satellite images.

Texture image segmentation

Image segmentation is a crucial stage in image processing that involves partitioning an image into homogenous areas based on shared pixel characteristics. This process has a significant impact on subsequent calculations of object properties and classification results [25]. To achieve the most accurate outcomes, it's essential to carefully select the optimal segmentation method and its parameters that are best suited for the specific problem at hand. Segmentation methods can be divided into automatic ones and interactive ones, the latter requiring the participation of the user. Automatic methods are also divided into two classes: 1) selection of image areas with certain properties specific to a specific subject area (marker methods, binarization); 2) dividing the image into homogeneous regions. The methods that divide the image into homogeneous areas are the most versatile, since they are not focused on a specific subject area or specific analysis tasks. Such algorithms are the most widespread in the computer vision, they include methods of water separation, the method of boundary selection, and methods based on multidimensional histogram clustering.

Assessing the quality of segmentation methods is not a straightforward process, as there is no universally accepted objective criterion for doing so. The optimal choice of method ultimately depends on the specific problem that needs to be addressed. To facilitate comparisons between segmentation techniques, reference image databases with known "ground-truth" segmentations can be utilized to evaluate the quality of each method's performance. One of the problems solved in the section is determining which segmentation methods and what values of their parameters are most suitable for use in the classification of multidimensional photogrammetric images with high spatial resolution.

The multi-scale segmentation of an image based on the calculated texture properties was used in this work to select areas of the image corresponding to a specific texture [26]. The multi-scale segmentation method (Multiresolution segmentation) is based on the technique of sequential merging of adjacent image elements. This is an optimization procedure that minimizes the average heterogeneity of image objects. The multi-scale segmentation method was applied to the satellite image shown in Fig. 1 to identify areas with similar textural characteristics. The calculated textural characteristics are taken as input data: the mean value, the variance, the homogeneity, the second angular moment, the correlation, the contrast, the divergence, and entropy. Fig. 3 shows the result of segmentation of the satellite image (Fig. 1) based on its textural features. The segments belonging to the texture corresponding to a plantation of trees are separated using the supervised classification method based on the parallelepiped algorithm. For controlled classification, reference areas are used, which are chosen by the operator in accordance with their belonging to a certain information class. The spectral features representing one class of pixels in the image are determined for each reference area. Each pixel is put to one of the classes by successively comparing it with all reference features. In controlled classification, information classes and their number are determined first followed by the determination of their corresponding spectral classes.

Parallelepiped algorithm. This classification algorithm is based on the Boolean logic and statistical indicators of the training sample in n spectral ranges. First, for each class c and range k, the average brightness value in the training sample is calculated 𝜇 ck and 𝜎 ck . Than the following rule is applied to classify the pixels of the image. A pixel belongs to a class if and only if its brightness BV ijk satisfies the following condition: 𝜇 ck − 𝜎 ck ≤ BV ijk ≤ 𝜇 ck − 𝜎 ck , (12) where c= 1,2,3... m is the class, and k= 1,2,3... m is the spectral range.

If we denote the lower and upper bounds of the inequality as:

𝐿 ck = 𝜇 ck − 𝜎 ck , 𝐻 ck = 𝜇 ck + 𝜎 ck , (13) the inequality will take the form:

𝐿 ck ≤ BV ijk ≤ 𝐻 ck . (14)

The set of points corresponding to inequality ( 14) forms a parallelepiped in the n -dimensional space of spectral features. If the brightness values of the pixels belong to the parallelepiped, the pixel belongs to this class. In this way, the segments of the image corresponding to the tree plantation are determined (Fig. 4).

Counting objects in satellite images

When the region of an image corresponding to a certain texture is found, it is possible to count the structural elements that make it up. This paper proposes an algorithm for counting objects in satellite images. The algorithm consists in selecting individual objects that create a texture based on their spectral characteristics.

The proposed algorithm consists of the following steps:

1. Selection and calculation of the spectral characteristics by which the search is carried out.

2. Binarization of the image according to the obtained spectral characteristics.

3. Finding connected components of the binarized image. 4. Selection of connected components corresponding to the parameters. 5. Calculation of the connected components. At the first step, the spectral characteristics according to which objects are searched for are selected. When searching for such objects as trees, it is possible to use vegetation indices as spectral features: NDVI, SAVI, ARVI, EVI. The NDVI index is commonly used [27].

The normalized vegetation index NDVI (Normalized Difference Vegetation Index) is used to solve problems using quantitative estimates of vegetation cover. The physiological state of the plant cover is largely determined by the content of chlorophyll and the level of moisture. It is advisable to use relative indicators of the state of vegetation, in particular forests, obtained taking into account spectral indices closely correlated with the level of plant chlorophyll and moisture [28]. The relative vegetation index NDVI is an indicator of the amount of photosynthetically active biomass and is calculated according to the formula:

𝑁𝐷𝑉𝐼 = 𝑁𝐼𝑅 − 𝑅𝐸𝐷 𝑁𝐼𝑅 + 𝑅𝐸𝐷 ,(15)

where NIR is the reflection of light in the near-infrared region of the spectrum, RED is the reflection in the red region of the spectrum (8 and 5 channels of WorldView-2 and WorldView-3 images, respectively). The calculation of NDVI is based on the two most stable areas of the spectral reflectance curve of vascular plants. The maximum absorption of solar radiation by chlorophyll takes place in the red region of the spectrum (0.6-0.7 μm); the IR region (0.7-1.0 μm) is the region of maximum reflection of leaf cellular structures.

A high photosynthetic activity, usually associated with dense vegetation, reduces the reflection in the red region of the spectrum and increses it in the IR one. Comparing these indicators with each other allows one to clearly distinguish plant objects from other environmental objects. The use of the normalized difference between the minimum and maximum reflectance allows one to reduce the effect of the image illumination and the radiation absorption by the atmosphere. Natural objects not associated with vegetation have a fixed NDVI value, which allows using this parameter for their identification. Depending on the objects on the Earth's surface, the NDVI index takes the following values: − negative values (when calculating in the range from -1 to 1) for water bodies; The main advantage of NDVI is the ease of calculation because no additional data are required except for the remote sensing data and the survey parameters. The NDVI index can be calculated on the basis of any images that have spectral channels in the red and near-IR ranges. The NDVI index has many modifications: SAVI, ARVI, EVI, etc. They are designed to reduce the impact of interfering factors.

Several high spatial resolution images with plantations of different types of trees were investigated (Fig. 5 and 7). The trees of a plantation usually have the same size, shape and spacing. This greatly facilitates the task of counting in comparison with counting trees in forest areas where the trees are located randomly. Three classes of objects are mainly present in the images used in this work: trees, soil, and tree shadows. Considering the fact that trees in gardens and plantations are arranged in an orderly manner, they have the same size. Due to their spectral characteristics, they are easily distinguished from the rest of the image pixels. Fig. 5 shows an image fragment and the NDVI index calculated on its basis (Fig. 5 (b)). Fig. 6 shows the histogram of the NDVI index. Using the method of automatic histogram splitting, it is possible to distinguish objects with a high index value, which corresponds to trees. The calculation of the vegetation index showed that in some images (Fig. 7) its value is too high in the entire area corresponding to the tree plantation. In this case, the histogram of the image is shifted from the center, as shown in Fig. 8. This leads to the impossibility of distinguishing the objects that correspond to trees using the histogram thresholding. In this case, it is possible to use other spectral characteristics of the image, in particular the NSVDI shadow identification index [29]:

𝑁𝑆𝑉𝐷𝐼 = 𝑆 − 𝑉 𝑆 + 𝑉 , (16)

where S is the image saturation and V is the brightness. To obtain the S and V components of the image, it is transformed from the RGB color model to the HSV color model. NSVDI takes values from -1 to 1, shadow areas have high index values. This work uses the NSVDI index to identify and count trees. For the image shown in Fig. 1, the proposed object counting algorithm was applied:

The first step of the algorithm is to find shadows in the image. This work used the image whose fragment is shown in Fig. 1. The input image was transformed into the HSV color model. The obtained image was then used to find the NSVDI index (Fig. 9). The pixels of the image that belong to the shadow areas take higher values, while the pixels that take smaller values belong to the non-shadow areas.

The second step of the algorithm is image binarization. An algorithm for automatically finding the optimal histogram threshold was applied to the image of the NSVDI index (Fig. 10). The image of the NSVDI index is divided into two classes according to the optimal threshold [30]. The third step of the algorithm is to find connected components of the obtained binary image. In this way, we obtained numbered segments of the image corresponding to the shadows (Fig. 11). The resulting image contains a large number of small segments that may correspond to noise, texture inhomogeneity, etc., while the large segments correspond to the shadows from individual trees.

The fourth step is to distinguish the connected components obtained at the second step by size. We discard the segments whose size is smaller than the given threshold, thus leaving only the segments that corresponding to the tree shadows. Further, it is not difficult to count the obtained segments. Fig. 11 shows the result of distinguishing individual objects using the connected components. The OpenCV computer vision algorithm library has been used for image processing and processing results visualization. As a result of the test image processing, 773 image elements corresponding to trees were obtained (Fig. 11).

Conclusions and plans for the future

The paper proposes a method of processing satellite images of high spatial resolution based on textural characteristics to solve the applied problem of identifying areas of the Earth's surface with agricultural plantations, in particular fruit tree plantations. An algorithm for the identification of individual objects in the selected area of an image for the purpose of their counting is proposed.

Texture segmentation was used to solve the problem of identifying areas of the Earth's surface where agricultural plantations are located. The calculation of textural properties of the image was performed on the basis of the gray-level co-occurrence matrix. The multiresolution segmentation method was applied to a satellite image to identify areas with similar textural characteristics. The calculated textural characteristics were taken as input data for segmentation: the mean value, the variance, the homogeneity, the second angular moment, the correlation, the contrast, the divergence, and the entropy.

The algorithm for selecting individual objects is based on the calculation of the spectral properties of images. Various spectral characteristics of images were analyzed. As a result, it was found that it is advisable to use vegetation indices and shade identification indices to identify trees in large agricultural plantations. On the basis of the spectral properties of the images and the result of their binarization, the objects that create the texture were distinguished. These objects were counted by means of software.

The actual number of trees was calculated on the studied samples and compared with the result of the algorithm. The obtained results showed a sufficiently high accuracy of the calculation. Further work will be devoted to improving texture segmentation, developing methods for automatically determining the minimum size of segments corresponding to the shadows from individual trees.