Agricultural Land Abandonment in Bulgaria: A Long-Term Remote Sensing Perspective, 1950–1980

: Agricultural land abandonment is a globally signiﬁcant threat to the sustenance of economic, ecological, and social balance. Although the driving forces behind it can be multifold and versatile, rural depopulation and urbanization are signiﬁcant contributors to agricultural land abandonment. In our chosen case study, focusing on two locations, Ruen and Stamboliyski, within the Plovdiv region of Bulgaria, we use aerial photographs and satellite imagery dating from the 1950s until 1980, in connection with ofﬁcial population census data, to assess the magnitude of agricultural abandonment for the ﬁrst time from a remote sensing perspective. We use multi-modal data obtained from historical aerial and satellite images to accurately identify Land Use Land Cover changes. We suggest using the rubber sheeting method for the geometric correction of multi-modal data obtained from aerial photos and Key Hole missions. Our approach helps with precise sub-pixel alignment of related datasets. We implemented an iterative object-based classiﬁcation approach to accurately map LULC distribution and quantify spatio-temporal changes from historical panchromatic images, which could be applied to similar images of different geographical regions.


Introduction
Agricultural land abandonment has been a significant land use change, attracting increased global attention since the 2010s [1]. It has been a concern for the European Union, as it has been forecasted that 11% of agricultural land in the EU is at a high risk of abandonment by 2030 [2]. The Common Agricultural Policy (CAP) has further accelerated the dual process of intensification of agriculture and agricultural land abandonment in the EU [3]. This political context and ecological momentum have facilitated large-scale EU-wide agricultural land abandonment examinations [4][5][6], including predictions and modelings for the near future [5,7,8]. Earth observation (EO) data, collected from various satellite platforms (starting in the 1970s with the Landsat program) and processed with remote sensing methods, has been a crucial and ever-developing aspect of these research efforts. In the last decades, for several local scale case studies, earlier monochromatic visual material (starting from the 1940s aerial photographs and the 1960s satellite imagery) has also been incorporated into Land Use Land Cover Changes (LULCC) examinations for more extended periods [9][10][11][12]. There have also been attempts to include historical cadastral maps as even earlier datasets, and in combination with aerial photographs [13][14][15][16][17].
We aim to quantitively analyze the land abandonment problem using very highresolution (VHR) historical aerial and satellite images with image processing techniques.
The satellite image dataset prepared for the Stamboliyski and Ruen study regions includes three images from the collections of declassified U.S. military intelligence missions from the KH optical reconnaissance satellites. We used images from Corona and Hexagon satellites operated between 1960 and 1972, and between 1971 and 1986, respectively. The archive of Corona's declassified satellite images was released and offered for public use in February 1995 [36]. The historical data obtained from the declassified archive of US satellite imageries have been useful in several remote sensing applications, including historical LULCC [37][38][39][40][41][42][43][44][45][46][47]; however, not for Southeast Europe and not in combination with aerial photographs, except for a recent example of combining 1945 aerial photographs with 1975 Hexagon imagery for deltaic areas in Greece [48]. The follow-up Hexagon reconnaissance mission's imagery was declassified in 2002 and 2011. Hexagon was the most sophisticated and the last film-based orbiting photo-reconnaissance satellite [49]. Although Hexagon imagery has high spatial resolution and advanced photogrammetric qualities, such as offering larger footprints without complex panoramic distortions [50], it has not been as widely used as its predecessors, such as Corona. Almost all the Hexagon project imagery taken by KH-9 system cameras with close to global coverage was declassified in 2011 (Declass-3). However, the imagery was only available at the National Archives and Records Administration (NARA) for personal inquiry until 2019, when the third phase of 14,000 declassified rolls was transferred to the USGS, EROS Center. Since then, they are gradually becoming available online for scanning and downloading (https://www.usgs.gov/centers/eros/ science/usgs-eros-archive-declassified-data-declassified-satellite-imagery-3 accessed on 1 October 2022). Previously declassified and online available Hexagon imagery covered mainly parts of North America and polar regions [51].
Historical aerial photographs have been used for LULCC in several studies [63,64], focusing on relatively small areas due to the coverage limitations of aerial photographs [65]. For over ten years, there have also been attempts to create long-term datasets by using aerial photographs and satellite imagery to examine historical LULCC [66], including our recent work, which combined historical maps with aerial photographs and modern satellite imagery [17,67]. Recently, there have also been attempts to automatize object-based image analysis of aerial photographs to conduct LULCC studies [68]. Nevertheless, as we already mentioned, using historical aerial photographs starting as early as the 1940s and the 1950s in combination with early satellite imagery to examine long-term LULCC is a rare undertaking. Furthermore, adding matching demographic data from the censuses is a novel approach to focus on rural depopulation and urban sprawl.

Study Areas and Datasets
We wanted to select two populated regions emphasizing two drivers of LULCC: rural depopulation and rapid population growth. Bulgaria has one of the fastest shrinking populations worldwide, and this trend is expected to continue at least until 2050 [69]. Bulgaria is in the middle of a demographic crisis on the national level. Furthermore, due to the hundreds of partially depopulated villages, it is possible to use the term demographic catastrophe at the rural level [70]. Moreover, the Plovdiv district has the country's highest share of the rural population [71]. Therefore, we decided to focus on the Plovdiv district. Thanks to the online availability of the National Register of Populated Places of the National Statistical Institute of Bulgaria, it is possible to detect populated places with drastic demographic changes from 1934 to 2020 [72]. Every populated place is given a unique Unified Classification of Administrative-Territorial and Territorial Units (UCATTU, EKATTE in Bulgarian) number; 63238 and 51980 for Ruen and Stamboliyski, respectively. After studying the long-term demographic dynamics of populated places in the Plovdiv district, we opted for Stamboliyski and Ruen as illustrative examples of drastic demographic changes in Table 1, possibly triggering dramatic LULLC. Stamboliyski experienced rapid urbanization between 1934 and 1965, had a stagnant population between 1975 and 1992, and then the population decreased. On the other hand, Ruen is an excellent example of  We looked for and found more candidates within the Plovdiv district with suitable demographic profiles, which could be examined for urban sprawl or rural depopulation as forms of historical LULCC. The availability of high-quality aerial and satellite imagery was our second constraint in choosing our populated regions. After successfully locating an adequate number of images of sufficient quality and at suitable intervals for Stamboliyski and Ruen, we finalized our case studies.
Our study areas are shown in Figure 1. The industrialized town of Stamboliyski, located 15 km west of Plovdiv city in South-Central Bulgaria, with an estimated population of 10,738 inhabitants at the end of 2020, represents a highly dynamic landscape over the last 70 years. The village of Ruen, located 17 km south of the city of Plovdiv, with an estimated population of 268 inhabitants in 2020, represents a relatively stable area with slow and low historical land use and land cover (LULC) change. The sizes of the Stamboliyski and Ruen study areas are about 14 km 2 and 6 km 2 , respectively. We acquired two panchromatic Corona images from the Key Hole 4B (KH-4B) mission from 1968, and two panchromatic Hexagon images from the Key Hole 9 (KH-9) mission from 1975 and 1980. We selected images extending over study areas (e.g., Stamboliyski and Ruen) with minimal cloud, fog, and snow cover. We purchased filmstrips of  We acquired two panchromatic Corona images from the Key Hole 4B (KH-4B) mission from 1968, and two panchromatic Hexagon images from the Key Hole 9 (KH-9) mission from 1975 and 1980. We selected images extending over study areas (e.g., Stamboliyski and Ruen) with minimal cloud, fog, and snow cover. We purchased filmstrips of chosen Corona and Hexagon images covering our study sites through the EarthExplorer graphical user interface of the U.S. Geological Survey portal. Each Corona image file consisted of four image frames (Figure 1a-d), with 7 to 10 subsets of each Hexagon image file. We list the exact dates, file names, and satellite missions of used Corona and Hexagon images, and the total number of frames and frames of each image we used covering the Ruen and Stamboliyski study areas, in Table 2. We present the footprints of the chosen images and the boundaries of Stamboliyski city and Ruen village in  Other than satellite images from declassified reconnaissance archives, we obtained monochromatic aerial photographs covering our study regions, which were acquired in 1945, 1952, 1965, and 1975. Except for the photograph covering Ruen in 1945, these aerial photographs can be obtained from the Geoportal of the Ministry of Defense of the Republic of Bulgaria (https://gis.armf.bg/ accessed on 1 October 2022). We acquired the 1945 Ruen aerial photo from the Luftbilddatenbank (https://www.luftbilddatenbank.de/ accessed on 1 October 2022), a company specializing in locating unexploded ordnance. We prepared aerial photo mosaics for two study sites for four dates 1945, 1952, 1965, and 1975, utilizing a total number of 21 aerial photographs. Detailed information about the Other than satellite images from declassified reconnaissance archives, we obtained monochromatic aerial photographs covering our study regions, which were acquired in 1945, 1952, 1965, and 1975. Except for the photograph covering Ruen in 1945, these aerial photographs can be obtained from the Geoportal of the Ministry of Defense of the Republic of Bulgaria (https://gis.armf.bg/ accessed on 1 October 2022). We acquired the 1945 Ruen aerial photo from the Luftbilddatenbank (https://www.luftbilddatenbank.de/ accessed on 1 October 2022), a company specializing in locating unexploded ordnance. We prepared aerial photo mosaics for two study sites for four dates 1945, 1952, 1965, and 1975, utilizing a total number of 21 aerial photographs. Detailed information about the number of aerial photographs combined to cover the study areas for the chosen years is shown in Table 3. To establish the correct map coordinate information for each cell in the aerial photographs, we used the spline (triangulation) transformation, a rubber sheeting method, to achieve high local accuracy [73]. Compared to polynomial functions, the spline transformation offers a very high goodness of fit to the reference data, but requires a far higher number of ground control points (GCPs) [17,74,75].
The rubber sheeting technique uses a fifth-order polynomial transformation to estimate optimal values using a mathematical function that minimizes overall surface curvature. The rubber sheeting method avoids sharp bending and provides a smooth surface that passes through the GCPs [76]. Spline-based transformations' more comprehensive image-warping options offer better accuracy for the geometric correction of historical topographic data with positional inaccuracies and geometric distortions [77,78].
We used the georeferencing tool in ArcGIS to conduct geometric correction of multitemporal and multi-source images. Georeferencing guarantees that satellite images and aerial photographs are correctly positioned in space, and that geometric distortion is avoided to the possible extent. Among the available datasets, we chose the one with the highest spatial resolution and widest coverage extent to be used as a reference for the further stages after conducting geometric correction using ground control points (GCPs) on the images and the Google Earth images. We chose the aerial photograph from 1975 as the reference data for the Ruen study location, since it encompasses the entire village in one frame and has a high spatial resolution. Then, we implemented a spline (triangulation) transformation to georeference the image, which provides higher accuracy for the geometric rectification of historical aerial photographs. Due to the 47 years of terrain changes between the aerial photograph from 1975 and the Google Earth image from 2021, we tried to select invariant GCPs during the period. We ensured the homogenous distribution of GCPs across the aerial photograph. Figure 3 shows the GCPs of the "1975 Ruen aerial photo". The zoomed views of three random GCPs are also exemplified in Figure 3.
For the geometric correction of the remaining aerial photographs encompassing the Ruen village, we used the georeferenced 1975 Ruen aerial photograph defined in WGS 84 Datum and UTM Zone 35N projection as the reference. Figure 4 shows GCP samples collected from aerial photographs in 1952 and 1965, and their counterparts on a 1975 Ruen aerial photograph georeferenced at the first stage.
We created the 1952 Ruen complete aerial photograph by mosaicking the two aerial photographs covering the Ruen village. Similarly, we made the mosaicked 1965 Ruen aerial photograph. We utilized Erdas Imagine image processing program for mosaicking. Figure 5 shows aerial photographs captured in 1952 and 1965, two for each year, and the result of mosaicking these photographs into the final mosaicked images, one for each year (1952 and 1965) covering most of Ruen village boundaries. one frame and has a high spatial resolution. Then, we implemented a spline (triangulation) transformation to georeference the image, which provides higher accuracy for the geometric rectification of historical aerial photographs. Due to the 47 years of terrain changes between the aerial photograph from 1975 and the Google Earth image from 2021, we tried to select invariant GCPs during the period. We ensured the homogenous distribution of GCPs across the aerial photograph. Figure 3 shows the GCPs of the "1975 Ruen aerial photo." The zoomed views of three random GCPs are also exemplified in Figure 3.   We created the 1952 Ruen complete aerial photograph by mosaicking the two aerial photographs covering the Ruen village. Similarly, we made the mosaicked 1965 Ruen aerial photograph. We utilized Erdas Imagine image processing program for mosaicking. Figure 5 shows aerial photographs captured in 1952 and 1965, two for each year, and the result of mosaicking these photographs into the final mosaicked images, one for each year (1952 and 1965) covering most of Ruen village boundaries. For the Stamboliyski study region, we geometrically corrected two aerial photographs of Stamboliyski from 1975 using the GCPs collected from the Google Earth image at the first stage. Each of the aerial photographs covers an area of approximately 27 km 2 . Figure 6 shows the distribution of GCPs utilized for geometric correction. It presents detailed views of some random GCPs on aerial photographs from the Stamboliyski region and their counterparts on Google Earth images. For the Stamboliyski study region, we geometrically corrected two aerial photographs of Stamboliyski from 1975 using the GCPs collected from the Google Earth image at the first stage. Each of the aerial photographs covers an area of approximately 27 km 2 . Figure 6 shows the distribution of GCPs utilized for geometric correction. It presents detailed views of some random GCPs on aerial photographs from the Stamboliyski region and their counterparts on Google Earth images.
We mosaicked the geometrically-corrected aerial photographs of Stamboliyski from 1975 and used them as a reference for rectifying other aerial photographs from 1952 and 1965. Figure 7 shows the two aerial photographs of Stamboliyski from 1975 and the mosaicked aerial photograph. The red color illustrates the boundaries of the Stamboliyski town.
We then georeferenced the aerial photographs encompassing Stamboliyski from 1952 and 1965, according to the geometrically corrected and mosaicked 1975 Stamboliyski aerial photograph. We produced the 1952 Stamboliyski aerial photograph by joining five aerial photographs, each covering an area of approximately 3 km 2 , while three aerial photographs from 1965, each covering an area of about 12 km 2 , were used to mosaic the 1965 Stamboliyski aerial photograph. Figures 8 and 9 show the individual aerial photographs from 1965 and 1952, respectively. Additionally, the mosaicked 1952 and 1965 aerial photographs are represented in Figures 8 and 9, respectively. We mosaicked the geometrically-corrected aerial photographs of Stamboliyski from 1975 and used them as a reference for rectifying other aerial photographs from 1952 and 1965. Figure 7 shows the two aerial photographs of Stamboliyski from 1975 and the mosaicked aerial photograph. The red color illustrates the boundaries of the Stamboliyski town. We then georeferenced the aerial photographs encompassing Stamboliyski from 1952 and 1965, according to the geometrically corrected and mosaicked 1975 Stamboliyski aerial photograph. We produced the 1952 Stamboliyski aerial photograph by joining five aerial photographs, each covering an area of approximately 3 km 2 , while three aerial photographs from 1965, each covering an area of about 12 km 2 , were used to mosaic the 1965 and 1965, according to the geometrically corrected and mosaicked 1975 Stamboliyski ial photograph. We produced the 1952 Stamboliyski aerial photograph by joining five ial photographs, each covering an area of approximately 3 km 2 , while three aerial ph graphs from 1965, each covering an area of about 12 km 2 , were used to mosaic the 1 Stamboliyski aerial photograph. Figures 8 and 9 show the individual aerial photogra from 1965 and 1952, respectively. Additionally, the mosaicked 1952 and 1965 aerial p tographs are represented in Figure 8 and Figure 9, respectively.

Geometric Correction of Satellite Images
After rectifying and mosaicking aerial photographs, we clipped the region covering Ruen and Stamboliyski from Corona and Hexagon satellite image frames, and georeferenced the clipped images according to the geometrically corrected and mosaicked 1975 Ruen and Stamboliyski aerial photograph as reference images. Figures 10 and 11 present the rectified declassified images of 1968, 1975, and 1980 that encompass the Ruen and Stamboliyski borders, respectively.

Geometric Correction of Satellite Images
After rectifying and mosaicking aerial photographs, we clipped the region covering Ruen and Stamboliyski from Corona and Hexagon satellite image frames, and georeferenced the clipped images according to the geometrically corrected and mosaicked 1975 Ruen and Stamboliyski aerial photograph as reference images. Figures 10 and 11 present the rectified declassified images of 1968, 1975, and 1980 that encompass the Ruen and Stamboliyski borders, respectively.

Geometric Correction of Satellite Images
After rectifying and mosaicking aerial photographs, we clipped the region covering Ruen and Stamboliyski from Corona and Hexagon satellite image frames, and georeferenced the clipped images according to the geometrically corrected and mosaicked 1975 Ruen and Stamboliyski aerial photograph as reference images. Figures 10 and 11 present the rectified declassified images of 1968, 1975, and 1980 that encompass the Ruen and Stamboliyski borders, respectively.    We list the spatial resolutions of the rectified aerial photographs and declassified satellite images of Ruen and Stamboliyski in Table 4.  Figure 12 shows spatial resolution comparisons between aerial and satellite images of Ruen on different dates. The aerial photograph from 1945 has the highest spatial reso- We list the spatial resolutions of the rectified aerial photographs and declassified satellite images of Ruen and Stamboliyski in Table 4. Figure 12 shows spatial resolution comparisons between aerial and satellite images of Ruen on different dates. The aerial photograph from 1945 has the highest spatial resolution, whereas the Corona image from 1968 has the lowest spatial resolution. As shown in Figure 11, buildings, roads, and even trees are detectable in satellite images and aerial photos of higher spatial resolution. At the same time, it is hard to detect buildings and borders of roads in low spatial resolution images.  Figure 12 shows spatial resolution comparisons between aerial and satellite images of Ruen on different dates. The aerial photograph from 1945 has the highest spatial resolution, whereas the Corona image from 1968 has the lowest spatial resolution. As shown in Figure 11, buildings, roads, and even trees are detectable in satellite images and aerial photos of higher spatial resolution. At the same time, it is hard to detect buildings and borders of roads in low spatial resolution images.

Classification of Aerial and Satellite Images
We classified the geometrically-corrected aerial photographs and satellite imageries to generate the LULC maps of the two study areas. We produced eleven LULC maps, six for Ruen and five for Stamboliyski. We used aerial photos to generate the LULC maps of 1945, 1952, 1965, and 1975, and the Corona and Hexagon satellite imageries were used to form the LULC maps of 1968 and 1980.
The object-based classification approach produced better results than pixel-based approaches, and multiple studies showed that object-based classifications have significantly superior accuracy [79,80]. The OBIA technique's segmentation step divides an image into useful image objects and meaningful entities composed of many nearby pixels. The most widely-used segmentation technique is multi-resolution segmentation (MRS) [81]. In the MRS method, a scale parameter determines the size of the resulting objects. It depends on the study area's characteristics and the image's spatial resolution. When a higher scale parameter value is selected, more merging is allowed, and the object is produced at a larger size.
In contrast, fine-scale values enable the formation of smaller objects. The ideal volume of a user-defined scale parameter for a study area should be set by testing the overall effect of segmentation results on object-based classification at various scales [82]. In the MRS algorithm, each pixel is considered an individual object, and pairs of image objects are combined throughout the pairwise clustering process in numerous subsequent steps in an interactive region-growing approach to generate more significant segments [83]. Merging the adjacent pixels is achieved by measuring the nearby object heterogeneity and the internal homogeneity of an object. The objects are merged until the assumption for internal homogeneity, and the threshold for surrounding object heterogeneity, are maintained. A combination of color (spectral values) and shape features make up the homogeneity criterion [84].
The objects are labeled in the second stage of the OBIA technique, known as the classification phase. The objects can be classified through the object-oriented supervised classification method, or the knowledge-based classification approach [85,86].
In the knowledge-based classification technique, the classification rules-also referred to as rule bases or knowledge models-are established and refined to identify the objects. Object identification is performed by utilizing a set of features, including texture (homogeneity, dissimilarity, entropy, contrast), shape (compactness, rectangular fit), position (distance, coordinates), extent (length/width), layer values (mean, standard deviation, mode), etc. An analyst's knowledge of an object's coverage, shape, texture, location, and size is utilized to determine the best threshold for the feature properties offered by e-cognition tools.

Classification Schema
The LULC map in this study provides eight classes and is defined using the land cover classification system adopted in European Space Agency (ESA) WorldCover project. During the ESA WorldCover project, detailed and high-resolution information on worldwide LULC cover was provided by the classification of Sentinel-1 and 2 images (ESA-WorldCover, 2020). The land cover classes in our work are derived from the ESA land cover classification scheme including cropland, built-up, bare/sparse vegetation, grassland, shrubland, water bodies, and tree cover. The explanation of each land cover class can be found in "WorldCover_PUM_v1.0" (https://worldcover2020.esa.int/, accessed on 1 October 2022), and we also provide examples from aerial photographs and satellite images for each land cover class in Figure 13. In addition to these land cover classes, we added a land cover class to our classification scheme to classify permanent cropland. In the ESA classification schema, perennial woody plants and plantations (e.g., fruit trees, nut trees, oil palm, olive trees) are classified as "tree cover". In this study, we separated permanent cropland from forest tree cover to examine changes in the agricultural region of the study sites.

Methodology
We first chose one photograph for each study region with high spatial resolution, encompassing the large extent of the study areas. Then, we implemented a geographic object-based image analysis (GeOBIA) using the eCognition ® program (Trimble Germany GmbH, Munich, Germany). The segmentation and classification phases of the OBIA technique are used to locate objects in an image and categorize the objects, respectively. For Ruen, the 1965 aerial photograph was the preferred image. Although the Ruen aerial photograph from 1945 has the highest resolution, the Ruen village is covered in >99% of the 1965 photograph, as opposed to only 93% in the 1945 photograph. The 1965 aerial photograph is the recommended image for Stamboliyski. Although the 1952 Stamboliyski aerial photograph has the highest resolution, the city of Stamboliyski is covered by more than 98% of the 1965 photograph, compared to only 74% in the 1952 photograph. It was also advantageous to pick the same year for both cases as the base year. Furthermore, between 1956 and 1965, Stamboliyski experiences the steepest rise and Ruen the most significant In addition to these land cover classes, we added a land cover class to our classification scheme to classify permanent cropland. In the ESA classification schema, perennial woody plants and plantations (e.g., fruit trees, nut trees, oil palm, olive trees) are classified as "tree cover". In this study, we separated permanent cropland from forest tree cover to examine changes in the agricultural region of the study sites.

Methodology
We first chose one photograph for each study region with high spatial resolution, encompassing the large extent of the study areas. Then, we implemented a geographic object-based image analysis (GeOBIA) using the eCognition ® program (Trimble Germany GmbH, Munich, Germany). The segmentation and classification phases of the OBIA technique are used to locate objects in an image and categorize the objects, respectively. For Ruen, the 1965 aerial photograph was the preferred image. Although the Ruen aerial photograph from 1945 has the highest resolution, the Ruen village is covered in >99% of the 1965 photograph, as opposed to only 93% in the 1945 photograph. The 1965 aerial photograph is the recommended image for Stamboliyski. Although the 1952 Stamboliyski aerial photograph has the highest resolution, the city of Stamboliyski is covered by more than 98% of the 1965 photograph, compared to only 74% in the 1952 photograph. It was also advantageous to pick the same year for both cases as the base year. Furthermore, between 1956 and 1965, Stamboliyski experiences the steepest rise and Ruen the most significant fall in their populations, according to census data (Table 1).
We segmented the 1965 aerial photographs into several objects using the multiresolution segmentation algorithm. In this segmentation method, a parameter called scale determines the size of resulting objects, and the shape and compactness parameters determine the boundaries of objects. The segmentation process of the aerial photographs was implemented with 0.2 shape, and 0.8 compactness parameter values. To minimize the number of objects facilitating the manual classification of the aerial photograph, the optimal values for scale parameters were applied according to each land cover class. The photograph was subdivided into the largest possible meaningful objects of a land cover class, with the scale appropriately defined. We applied different values for scale parameters, including 800, 600, 400, 200, 100, and 50, at different stages. At each step of the multi-scale segmentation approach, objects assigned to one of the land cover classes were classified and labeled with the related class name. However, others not applicable to classification due to encompassing more than one land cover type were segmented further using the smaller scale. The scale parameter of 800 was suitable for segmenting objects for most of the "tree cover" and "cropland" land covers. As seen in Figure 14, there was no need to use smaller scales to segment lands covered by trees (i.e., tree cover), since the objects resulting from the scale of 800 meaningfully encompass tree cover and separate it from neighboring objects of other land cover classes. We then manually assigned the category of "tree cover" to these objects. use smaller scales to segment lands covered by trees (i.e., tree cover), since the objects resulting from the scale of 800 meaningfully encompass tree cover and separate it from neighboring objects of other land cover classes. We then manually assigned the category of "tree cover" to these objects. Subsequently, the 600, 400, and 200 scale parameters were utilized to segment and classify image objects. Most of the shrubland, grassland, water bodies, bare/sparse vegetation, and permanent crops were classified at these stages. We used the smaller objects resulting from scales of 100 and 50 in the multi-resolution segmentation configuration to classify the remaining regions. Mainly, the objects relating to the built-up land cover category were classified at this stage. While land covers with simple outlines, such as croplands, water bodies, and forests, are better suited to being segmented and classified at larger scales, smaller objects associated with built-up structures with complex boundaries necessitate the creation of smaller image segments to distinguish built-up structures from neighboring regions. Figure 15 shows the objects resulting from the last step of the multi-scale segmentation procedure, which is suitable for segmenting and classifying parts such as roads, buildings, grasslands within the forest area, and shrubland within agricultural land. Subsequently, the 600, 400, and 200 scale parameters were utilized to segment and classify image objects. Most of the shrubland, grassland, water bodies, bare/sparse vegetation, and permanent crops were classified at these stages. We used the smaller objects resulting from scales of 100 and 50 in the multi-resolution segmentation configuration to classify the remaining regions. Mainly, the objects relating to the built-up land cover category were classified at this stage. While land covers with simple outlines, such as croplands, water bodies, and forests, are better suited to being segmented and classified at larger scales, smaller objects associated with built-up structures with complex boundaries necessitate the creation of smaller image segments to distinguish built-up structures from neighboring regions. Figure 15 shows the objects resulting from the last step of the multi-scale segmentation procedure, which is suitable for segmenting and classifying parts such as roads, buildings, grasslands within the forest area, and shrubland within agricultural land. at larger scales, smaller objects associated with built-up structures with complex bou ries necessitate the creation of smaller image segments to distinguish built-up struc from neighboring regions. Figure 15 shows the objects resulting from the last step o multi-scale segmentation procedure, which is suitable for segmenting and classi parts such as roads, buildings, grasslands within the forest area, and shrubland w agricultural land. At the end of the first phase of image classification, we had the land cover map of the Ruen and Stamboliyski from the 1965 aerial photographs. The additional aerial photographs and satellite images were then segmented and classified using the land cover map of the 1965 aerial photographs. We implemented the segmentation and classification of the images in multiple stages. First, an image was segmented using the scale parameter set to 80, and the shape parameter and the compactness were set as 0.7 and 0.3, respectively. The classes given to objects in the 1965 land cover map were then assigned to the segmented objects in stages, one class at a time. We later checked the classified objects manually, and we edited the segments which were misclassified. Misclassified objects can include segments that have altered over time, or items that appear different in a current image due to differences in quality or spatial resolution. Figure 16 shows an example of a land cover class assignment based on the 1965 map, and subsequent revisions based on the 1952 aerial photograph from Ruen. It depicts the tree cover class allocated to the 1952 aerial shot based on the 1965 land cover map, and the altered pieces of this land cover class that were later assigned to different land cover classes.
At each phase, we assigned one land cover class to the unclassified objects based on the 1965 land cover map of Ruen and Stamboliyski, and modified images in the process until all the objects were classified. As shown in the legend of Figure 16, we had a class called "segment again" that we utilized for objects requiring further segmentation due to spanning regions with various land covers. We set scales of 50 and 25 to segment objects in the "segment again" class. While applying a region's prepared land cover map to other photographs from other dates saves time and effort, assigning a class at a time narrates the number of misclassified objects to edit and reduces human errors. manually, and we edited the segments which were misclassified. Misclassified objects can include segments that have altered over time, or items that appear different in a current image due to differences in quality or spatial resolution. Figure 16 shows an example of a land cover class assignment based on the 1965 map, and subsequent revisions based on the 1952 aerial photograph from Ruen. It depicts the tree cover class allocated to the 1952 aerial shot based on the 1965 land cover map, and the altered pieces of this land cover class that were later assigned to different land cover classes. At each phase, we assigned one land cover class to the unclassified objects based on the 1965 land cover map of Ruen and Stamboliyski, and modified images in the process until all the objects were classified. As shown in the legend of Figure 16, we had a class called "segment again" that we utilized for objects requiring further segmentation due to spanning regions with various land covers. We set scales of 50 and 25 to segment objects in the "segment again" class. While applying a region's prepared land cover map to other photographs from other dates saves time and effort, assigning a class at a time narrates the number of misclassified objects to edit and reduces human errors.

Accuracy Assessments-OA and Kappa
Satisfactory overall accuracies (>85%) for the classification results are prerequisites for accurate land cover mapping [75]. We assessed the accuracy of our classification results with a set of 20 reference point samples for each LULC class. We targeted well and homogenously distributed sample points within the entire mapping area. We used the original aerial photographs and satellite images as the reference data. Tables 5 and 6 present the overall accuracy and Kappa metrics of land cover maps of Ruen and Stamboliyski. Tables 7 and 8 show the per-class accuracy results for the aerial photographs and satellite images of Ruen and Stamboliyski. Of the nine LULC classes, other than tree cover, shrubland, and grassland land classes in which marginal mapping confusion occurred, other LULC classes had high per-class accuracy values (>85%).   Figure 17 presents the LULC map of the Ruen site in 1945, 1952, 1965, 1968, 1975, and 1980, the classification results of the aerial photographs (1945-1952-1965-1975), Corona satellite image (1968), and Hexagon imagery (1980). Table 9 presents the statistical analysis of the land cover classes in each of the land cover maps.  Table 1) resulted in a decrease of 1/3 in cropland for Ruen in the last decade under investigation. This development is in accordance with expectations of rural depopulation's retarded impact on agricultural land abandonment. Figure 18 presents the LULC map of the Stamboliyski site in 1952, 1965, 1968, 1975, and 1980, the classification results of the aerial photographs  (1952-1965-1975), Corona satellite image (1968), and Hexagon imagery (1980). Table 10 presents the statistical analysis of the land cover classes in each of the land cover maps.  Figure 17 presents the LULC map of the Ruen site in 1945, 1952, 1965, 1968, 1975, and 1980, the classification results of the aerial photographs (1945-1952-1965-1975), Corona satellite image (1968), and Hexagon imagery (1980). Table 9 presents the statistical analysis of the land cover classes in each of the land cover maps.     Table 1) resulted in a decrease of 1/3 in cropland for Ruen in the last decade under investigation. This development is in accordance with expectations of rural depopulation's retarded impact on agricultural land abandonment. Figure 18 presents the LULC map of the Stamboliyski site in 1952, 1965, 1968, 1975, and 1980, the classification results of the aerial photographs (1952-1965-1975), Corona satellite image (1968), and Hexagon imagery (1980). Table 10 presents the statistical analysis of the land cover classes in each of the land cover maps.    Table 1) and then leveling of population growth (from 13,313 in 1975, to 13,449 in 1985; see Table 1) can be seen reflected in a rapid and then gradual increase in the built-up area. The mirror image of this dynamic can be observed in the fast fall but then leveling of cropland.

Discussion
Accurate identification of the spatio-temporal distribution of LULC classes is essential to understanding past human activities and their impact on the environment, such as agricultural land abandonment, urbanization, and deforestation. Although there are complicated algorithms to process recent high-resolution aerial photographs and satellite images, this is a challenging task for historical photographs and satellite platforms, such as Corona and Hexagon. The problem is mainly due to their limited spectral and radiometric resolutions and the availability of various geometric distortions. With the recent public availability of Hexagon images, reliable methodologies for the geometric correction and classification of these images and integrated usage of these images with another historical geospatial dataset became important. Although there are several studies with Corona KH-4 images, studies with Hexagon images are fairly limited [43,48,87,88]. First studies with Hexagon KH-9 images were generally conducted for glacier monitoring, archaeology, and digital surface model generation, similar to Corona images [48,[50][51][52]54,57,[89][90][91]. The number of LULC mapping studies with Key Hole satellites is minimal, and a few studies were conducted with Corona images [41,43]. Utilizing visual interpretation methods, Fekete [45] used all Key Hole satellite series, including KH-9 Hexagon, to determine natural hazard risk and urban growth. Rashid and Aneaus used high-resolution declassified satellite data from 1965 (from Corona mission) and 1980 (Hexagon mission) to quantify land system changes. They employed on-screen digitization, with cognitive inputs from analysts [92]. Prokop utilized a Corona image from October 1965 to detect long-term land-use changes and they applied a visual interpretation technique for the LULC mapping [93]. Our study is unique in terms of the proposed methodology and the usage of Hexagon data in conjuction with other historical geospatial data for the LULC mapping perspective. We applied a feasible geometric correction approach by implementing rubber sheeting and collecting reference points from the Google Earth system that can be transferred to any other geographic region if similar datasets are available. Most of the studies in the literature use manual vectorization techniques, which are time-consuming and subjective. We propose an iterative GeOBIA approach to accurately generate multi-temporal LULC maps from multi-modal data, which researchers can use for various purposes in the case of similar remotely sensed data. Specifically, scale parameters specific to different land categories could be applied to various geographical areas in the case of using similar geospatial data sources. Historical census data is also an essential information source to understand and interpret historical LULC changes.
Implementing the ESA WoldCover project's LULC classes with a modification in tree classes to determine permanent agriculture classes is also unique to our study, which is essential to monitor agricultural changes apart from forests. After revising current ESA WorldCover data for tree cover classes and quality controlling its accuracy, a comparison could also be made to detect LULC changes from the 1940s to today, which will be considered in our future work. However, it is important to emphasize that ESA WorldCover product resolution is 10 m, which is lower than our very-high resolution historical geospatial data; therefore, generalization might be needed, resulting in the loss of some spatial details. Another alternative is to use VHR satellite or UAV images to find the long-term changes from the past to today, something that we will consider in our future work to find up-to-date VHR data.

Conclusions
Agricultural land abandonment is a massive threat to the population geography of Bulgaria, which has one of the fastest shrinking populations on a national scale worldwide. Combined with an ongoing rural depopulation, agricultural land abandonment will probably cause a rural population collapse in Bulgaria in the coming decades. Therefore, it is crucial to understand the timing and magnitude of these decisive demographic and economic changes. Our study roughly examines the first three decades of communist rule in Bulgaria, during which drastic social engineering practices, including radical collectivizations in the agricultural sector, have occurred. In both of our chosen locations, we can see a decrease in cropland in the period. This finding is striking at first glance, since we opted for Ruen as a representative example of stark rural depopulation and Stamboliyski for urban development. The fact that we observe a fall in Ruen and stagnancy in Stamboliyski, in the areas used for agricultural production in the period, hints at a systematic issue, most probably politically caused, and long in the making of agricultural land abandonment.

Conflicts of Interest:
The authors declare no conflict of interest.