Examining Land Use/Land Cover Change and Its Prediction Based on a Multilayer Perceptron Markov Approach in the Luki Biosphere Reserve, Democratic Republic of Congo

: Villages within the Luki Biosphere Reserve and the surrounding cities have undergone rapid demographic growth and urbanization that have impacted the reserve’s natural landscape. However, no study has focused on the spatiotemporal analysis of its land use/land cover. The present research aims at providing a comprehensive analysis of land use/land cover change in the Luki Biosphere Reserve from the year 1987 to 2020, and to predict its future change for the year 2038. Landsat images were classiﬁed in order to provide land use/land cover maps for the years 1987, 2002, 2017 and 2020. Based on these maps, change detection, gradient direction, and landscape metric analyses were performed. In addition, land use/land cover change prediction was carried out using the Multilayer Perceptron Markov model. The results revealed signiﬁcant land use/land cover changes in the Luki Biosphere Reserve during the study period. Indeed, tremendous changes in the primary forest, which lost around 17.8% of its total area, were noted. Other classes, notably savannah, secondary forest, built-up area, fallow land and ﬁelds had gained 79.35, 1150.36, 67.63, 3852.12 hectares, respectively. Based on the landscape metric analysis, it was revealed that built-up areas and fallow land and ﬁelds experienced an aggregation trend, while other classes showed disaggregation and fragmentation trends. Analysis further revealed that village expansion has signiﬁcantly affected the process of land use/land cover change in the Luki Biosphere Reserve. However, the prediction results revealed that the primary forest will continue to increase while built-up area, fallow land and ﬁelds will follow a trend similar to a previous one. As for secondary forest and savannah, the forecast revealed a decrease of the extent during the period extending from 2020 to 2038. The present ﬁndings will beneﬁt the decision makers, particularly in the sustainable natural resources management of the Luki


Introduction
Since the International Geosphere and Biosphere Program (IGBP) and the International Human Dimensions Program (IHDP) on global change initiated a research program on land use/land cover change, numerous studies in this field have been conducted at a global scale due to current issues in climate and global change [1]. The concept of land use/land cover is comprised of two dimensions: land cover, which describes the physical state or the surface of the soil (the type of vegetation, water, and soil), and land use, which illustrates the way human beings utilize the land. Land use relates to different activities that are undertaken on lands, including agriculture, pasture, or settlement [2][3][4]. A change to provide accurate information on LUC change to facilitate the policy-making process in elaborating good strategies in land use planning. On the other hand, the National Research Council [31] has shown that each of these models has its strengths and limitations in predicting LUC change. Therefore, there is need for an integrated approach by these models to overcome the corresponding limitations [27]. Recent studies have integrated a number of prediction methods to simulate LUC change [32,33]. In the present study, the MLPNN-based MCM was used to simulate the future LUC change. The MLPNN-MCM integrates both multiple layer perceptron neural network and Markovian chain model and has been successfully used in various studies [34,35]. Stochastic Markov and multiple layer perceptron neural network techniques appear to provide better prediction results on LUC change compared to the regression models [29,36,37] due to the ability of the MLPNN to predict several transitions simultaneously. MCM can automatically generate numerous parameter values that take a small number of data for training and reduce the calibration time of the model [38]. In addition, it can combine various parameters that can impact the numerous LUC transitions [27,29]. In addition, the MCM also provides the transition potentials of each LUC class to simulate the future LUC changes based on previous changes [39,40].
The goals of the present study are to: (1) examine and quantify land use/land cover change of the Luki Biosphere Reserve, and for all typical zones, including Core zone, Buffer zone, and transition zone; (2) depict spatiotemporal characteristics of village expansion in the reserve and its impacts on land use/land cover change; (3) analyze landscape pattern change of these villages; (4) predict future trends of land use/land cover change in the reserve. From a management standpoint, the present study will provide insights on previous LUC change and its future trends in the Luki Biosphere Reserve. The resulting information can help decision makers to solve issues that may arise, including whether forest losses will increase due to the increase of inhabitants in the region. This can help to elaborate effective strategies of land use management for sustainable development of the region.

Study Area
The present research was conducted in the Luki Biosphere Reserve, one of the protected areas in the Democratic Republic of Congo ( Figure 1). Luki extends over an area of about 32,000 hectares and is located in the southwest part of the country, in the Kongo central province. The annual average rainfall is 1095.66 mm (Figure 2), mostly falling between October and May, with a marked dry season between May and September. The annual average temperature is 24.65 • C. The main vegetation types are tropical forests. The relatively high precipitation in the region allows the local communities to resort to year-round agricultural activities to support their families.
A large stretch of the national road crosses the Luki Reserve from Manterne to Kinzao-Mvuete and is a compulsory passage for vehicles. This partly explains the cause of the rapid growth of the villages on this route. It is worth noting that many illegal operators within the reserve live in these villages, compromising the sustainability of its natural resources.

Image Classification and Accuracy Assessment
Remotely sensed data (Landsat images and a digital elevation model with a spatial resolution of 30 m) were freely downloaded from http://glovis.usgs.gov/ and https: //code.earthengine.google.com websites. From Google Earth Engine, Landsat images of the years 1987, 2002, 2017 and 2020 were downloaded. Also, the digital elevation model for the year 2000 was downloaded. In 2021, field data (500 observations per LUC category) were collected using a GPS receiver to provide a better understanding of different LUC categories in the study area.   Table 1 provides the description of the Landsat images used in the present research. Landsat images of the years 1987, 2002, 2017, and 2020 were radiometrically and atmospherically corrected using the FLAASH method. The landscape of the study area was then classified into five LUC types, including primary forest, savannah, fallow land and fields, secondary forest, and built-up area. The training and testing data of each year were created based on the field data collected for each LUC type in 2020 and the Google Earth images. Supervised classification was performed on the four images using the maximum likelihood algorithm. The maximum likelihood offers better classification results in the field of LUC study [42]. The accuracy assessment for the 1987, 2002, 2017, and 2020 maps was achieved using the Kappa statistic [43][44][45]. For all image processing, ENVI 5.3, QGIS, and ArcGIS were used. In recent decades, a number of studies have revealed the capability of the Artificial Neural Network algorithm in calibrating the LUC change modeling to predict its future trends [46][47][48][49]. The Artificial Neural Network algorithm was created by imitating the nervous system of the human being [50]. The present study used the Artificial Neural Network called multilayer perceptron neural network [51]. In the MLPNN algorithm, neurons or nodes are interconnected to create a complex neural network that is made up of three layers including an input layer, one or several hidden layers, and one output layer [52], as shown in Figure 3. Moreover, this model is made up of two types of input variables: the LUC images as dependent variables and the drivers of LUC change as independent variables. Each neuron in the input layer receives one of the input variables to produce an output value to the next layer. In that way, the input variables behave as signals that go through the layers in the MLPNN and finally produce output values. The mathematical formulation of the weighting process of the nodes in the input layers is presented below: o ij represents the weights between node i and j. o i represents the output resulting from the node i. The mathematical computation of the output from a given node j is presented as follows [52]: In the above formula, the nonlinear sigmoid function (f ) is supposed to be applied on the weighted sum before the signal goes to the next layer. Nevertheless, in order to perform the correction of connection weights, the following formula is used [53]: η represents the learning rate parameter; δ j represents an index of the rate of change of the error; and α represents the momentum parameter.
Meanwhile, the performance of the neural network is affected by some other factors [52], including the number of nodes, number of training samples, number of iterations and the learning rate. Generally, the training accuracy is influenced by the sample size and iteration number. In fact, a small number of training samples may not represent the pattern of each class, while a great number of training samples may result in an overlap. Additionally, very few or too many iterations can lead to poor generalization of the model [53]. Indeed, the model runs with a number of iterations until the minimization of network errors [52]. The objective of performing the network training is to provide appropriate weights in order to perform the classification of unknown pixels [53]. The error rate is computed following the formula of the Root Mean Square (RMS) Error [54]: N represents the number of elements; i represents the index for elements; e i represents the error of the ith element; t i represents the target value (measured) for the ith element; and a i represents the calculated value for the ith element.
Finally, to better represent the training data, the neural network requires a more significant number of nodes in the hidden layer [52]. The mathematical expression to compute the number of nodes in the hidden layer is presented as follows [53]:

Markov Chain Analysis
Markov chain techniques have been successfully applied in predicting LUC change. This model essentially examines the trend of transition probability of LUC classes in a given landscape over the specified number of time units. In this model, the LUC change can be described as a set of alternative states over time. This assumes that only the system state at time t can be used to determine its state at a given time (t + 1). In the present study, the Markov chain analysis was performed for the period between 1987 and 2002 to produce the transition probability matrix [56]. These matrices provide the transition probabilities for all LUC categories over the specified time and were used to predict the LUC map of 2017. The Markov chain analysis can be described using the following mathematical expression [57,58]:  [55].
N h represents the number of hidden nodes; N i represents the number of input nodes; N o represents the number of output nodes.
S (t) and S (t+1) represent the state of the system at time t and t + 1, respectively. P ij represents the transition probability matrix standing for the probability of the LUC class to transition from the state i to the state j, from the time t to the time t + 1. P ij is mathematically computed using the following expression [59,60]: P ij = P 11 P 12 P 13 . . . P 1n P 21 P 22 P 23 . . . P 2n P 31 P 32 P 33 . . . P 3n . . . · · · · · · . . . · · · P n1 P n2 P n3 . . . P nn (7) (0 ≤ P ij ≤ 1 and n ∑ j=1 P ij=1,i,j=1, 2,...,n ). P n represents the state probability of any time; n is the LUC type number; S is LUC status, t; t + 1 is the time point. Thus, any class with a low transition value will be expected to have a probability close to (0), while a high transition will have probabilities close to (1) [60].

Implementation of the MLPNN_Markov Modeling
During the transition potential modeling, the major LUC changes that occurred in the Luki Biosphere Reserve from 1987 to 2002 were identified and examined. This previous analysis was performed for the purpose of defining the transition categories. Consequently, the transition categories were well-defined, and the performance of the predictors was evaluated by means of the Cramer test [61]. In the present research, predictors with Cramer values >0.15 were used for predicting LUC change, as suggested by Hamdy et al. [62]. In fact, the range values of the Cramer test were between 0 and 1. Predictors with values close to 1 are reported to be more associated with the defined transition categories.
Thus, their association with the defined transition categories are expected to be large. With respect to measuring the strength of association between predictors, the Cramer test is a nonparametric statistical analysis whose application does not relate to the variable distribution [63][64][65]. The Cramer's V test is computed as follows: V = Cramer's index, n = sample size, q = smallest value in the rows and columns of the LUC image, O = observed frequency for a class, and E = expected frequency in the corresponding class.
Predictors such as slope, altitude, distance from roads, distance from villages, and economic indicators have often been used for the purpose of predicting LUC change [1,30,[65][66][67]. Thus, on the basis of the LUC change that occurred in the Luki Biosphere Reserve and the existing literature, five predictors were selected for the MLPNN_Markov modeling. These included the following variables: Distance from roads, distance from villages, distance from farmland, slope, and altitude ( Figure 4). A digital elevation model was downloaded from the USGS website, and from this raster the slope raster was computed. Distance from roads and villages was computed by using the proximity tool in QGIS software. The vector data of roads and villages were converted to a raster file, then we calculated the distance of each pixel based on the road and village pixel. As for distance from existing farmlands, pixels of farmland for the corresponding period of prediction were extracted. The distance raster from farmland was then computed using the proximity tool in QGIS software.
After performing the Cramer's V test, the MLPNN model was trained to compute the probabilities for each cell to transition from one type of LUC class to another from the year 1987 to 2002. Therefore, to improve the model performance, it is usually required to customize some parameters. In this case, the number of iterations was set to 10,000. The maximum sample size was customized based on the minimum number of pixels that was expected to transition from the initial date (1987) to the final date (2002). Afterwards, the learning rate, momentum, and hidden layers were customized by their default values. A given weight was then attributed to each major transition. The training process then started by dividing the sample data into validation and training, and then analyses on the accuracy of validation and training data. Finally, the training process was completed when the best accuracy was reached.
After the MLPNN model was finished, each transition contained the appropriate weight. These weights were incorporated in the Markov transition matrix for future prediction. Finally, the 2017 and 2038 maps were predicted using the transition potential maps derived from the MLPNN model and the Markov chain analysis. The workflow of the training and simulation process is illustrated in Figure 5.

Variables Selection
Based on the Cramer's V test, the selection process of variables was performed. Table 2 shows the results related to the test of the driver variables, provided by the Cramer's V test for the prediction of LUC change in the Luki Biosphere Reserve. Since their Cramer values were greater than 0.15, all the variables were used in the present study. These results revealed that all the explanatory variables had an influence on the LUC change process in the Luki Reserve, even though the elevation, distance from roads, and distance from existing fallow land and fields were revealed to be the most important. In general, areas close to existing fallow land and roads and with a high elevation were expected to experience more significant changes. Additionally, most areas close to villages were expected to experience a change in the region.

Model Validation
In this study, the validation process used the Pontius et al. [68] method. Here, the actual map of time 1 and time 2, and the predicted map of time 2 were used to generate the validation result, notably the agreement and disagreement component [69]. The agreement component was composed of two parts, i.e., correct rejection (persistence simulated correctly) and hits (change simulated correctly). The disagreement component was made up of three parts including misses (change predicted as persistence), false alarms (persistence predicted as change), and wrong hits (change predicted as a change to the wrong category). Thus, the reference map of 2002, the reference map of 2017, and the predicted maps of 2017 were used for that purpose.

Gradient Direction Analysis
Upscaling methods are widely used in the field of landscape ecology and gradient direction analysis and have been found to be very useful in examining the trend direction of land use/land cover change over time [70]. Firstly, the method starts by drawing several concentric rings on urban areas following a well-defined interval. The urban area and all land use/land cover classes that are less impacted by human activity (forest, water body, etc.) should be completely covered by the extent of the largest concentric ring. The next step involves drawing 16 fans through the rays extended from the urban center at defined intervals of 22.5 • . Thus, a set of segment zones is generated at the level of the intersection between the concentric rings and the fans. Finally, the land use/land cover class with the largest proportion in each segment is selected to represent that segment.

Landscape Metrics Analysis
Landscape metrics analysis has been widely applied in the field of landscape ecology to quantify changes that may have occurred over time. Landscape structure analysis based on numerous metrics has been commonly used to examine spatiotemporal change in various land use/land cover classes. These changes can be examined and quantified following numerous landscape metrics related to the number, size, shape, and other parameters of the landscape that can be obtained from remote sensing data [71,72].
This study used four landscape metrics (Table 3), notably the Landscape Division Index (DIVISION), Patch Cohesion Index (COHESION), Aggregation index (AI), and Contagion index (CONTAG), to examine the spatial structure of the landscape [70]. Generally, the DIVISION and COHESION indices are used to measure the connectivity of the landscape habitat. The COHESION index ranges between 0 and 100. The null value is observed for the highly subdivided landscape and is physically less connected. The DIVISION varies between 0 and 1. The zero value is observed for the homogeneous landscape (landscape with a single patch), while the value of 1 is observed when the focal patch type consists of a single small patch, one cell in area. The aggregation index provides a measure of the specific class and is ranged between 0 and 100. The value of the index reaches 100 for a landscape whose focal patch type is maximally aggregated into a single compact patch. When the different types of land use are aggregated and geographically self-connected, they are characterized by high values of COHESION and AI indices. On the other hand, landscapes with land use types disconnected from each other are characterized by high values of the division index. CONTAG measures both the interspersion and the dispersion of patches at the landscape level. CONTAG values range from 0 to 100. The null values are observed in the landscape where patch types are characterized by minimum disaggregation, whereas it reaches 100 in the case of maximum aggregation. The FRAGSTATS software (version 4.2) as well as the landscape ecology plugin from QGIS [73], were used to carry out these quantitative analyses. Table 3. Landscape metrics used in this study.

Definition
Description Reference a ij A 2 a ij = area (m 2 ) of patch ij. A = total landscape area (m 2 ). [74-76] p ij * = perimeter of patch ij in terms of number of cell surfaces. a ij * = area of patch ij in terms of number of cells. Z = total number of cells in the landscape.
[77] AI = g ii max→g ii (100) g ii = number of like adjacencies (joins) between pixels of patch type (class) I based on the single-count method. max→g ii = maximum number of like adjacencies (joins) between pixels of patch type (class) i based on the single-count method. [78]   P i = proportion of the landscape occupied by patch type (class) i. g ik = number of adjacencies (joins) between pixels of patch types (classes) i and k based on the double-count method. m = number of patch types (classes) present in the landscape, including the landscape border if present. [79]

Land Use/Land Cover Change of the Luki Biosphere Reserve
The past patterns (1987-2020) of the land use/land cover change were examined using Landsat images and the change detection technique. The overall classification accuracy was more than 80% for 1987, 2002, 2017, and 2020 classified images ( Table 4). The landscape of the Luki Biosphere Reserve was classified into different land use/land cover classes including primary forest, secondary forest, fallow land and fields, savannah and built-up area (Table 5 and Figure 6). Primary forest remained the dominant class of the landscape since it was the major class. The class declined from approximately 87.59, 78.29 to 72.02% of the total area of Luki in 1987, 2002 and 2020, respectively. During the last four decades, primary forest lost 4811.93 hectares.
However, other land use/land cover classes experienced positive change from 1987 to 2020, except for savannah, which underwent negative change from 2002 to 2020. Fallow land and fields, secondary forest, savannah, and built-up area increased by 3852.12, 1150.36, 79.35, and 67.63 hectares, respectively. Secondary forest, built-up, and fallow land and fields quickly increased from 1981 to 2021, although most of the newly increased land use/land cover class was fallow land and fields.

Transition among Land Use/Land Cover Types from 1987 to 2038
The transition matrix (Table 6) provides the comprehensive direction change of different land use/land cover classes of the Luki Biosphere Reserve over the study period. From 1987 to 2020, the greatest loss of area was observed in the primary forest class. Additionally, the largest area of primary forest was converted into fallow land and fields, followed by secondary forest, savannah, and built-up area. Indeed, more than 13.26% of primary forests was converted into fallow land and fields due to the expansion of agricultural activities in the reserve. Further, the fallow land and field class recorded the largest areal gain, mainly converted from primary forest, secondary forests (7.8%), savannah (1.5%), and built-up area (0.8%). However, of the primary forest converted into fallow land and fields, 4% naturally evolved into secondary forest in the absence of human activities. Moreover, the overexploitation of primary forest areas over a long period led to their conversion into savannah. From 1987 to 2038, the highest amount of change is still expected to be observed converting primary forest to fallow land and fields; it will be even greater (16%) than that of the period between 1987 and 2020. Over this period, 3.1, 0.2, and 0.09% of the primary forest area will be converted into secondary forest, savannah, and built-up area, respectively. However, fallow land and fields will continue to increase in area coverage, 23.9, 29.7, and 7.1% of the increase are attributed to conversion from secondary forests, savannah, and built-up areas, respectively.

Land Use/Land Cover Change at Different Zonal Levels of the Luki Biosphere Zone
The landscape change in different zones of the Luki Biosphere Reserve is presented in Tables 7-9. The spatial pattern of change in each zone is also illustrated in Figure 7.  At the Core zone level, from 1987 to 2020 the primary forests lost 287.5 hectares while the fallow land and fields and secondary forest gained 194.5 and 93 hectares, respectively. This could be due to the illegal activities within this zone, including agricultural practices by the local communities, cutting timber for handcrafting, and charcoal for commercial use. With respect to the Buffer zone, a loss of 818.93 hectares of primary forests and 104.54 hectares of secondary forests was revealed. This could have resulted from the slashand-burn agriculture practiced in this area. Thus, the fallow land and field class gained 704.73 hectares. Savannah and built-up areas also experienced an increase of 6.92 and 2.79 hectares, respectively. The Transition zone underwent dramatic changes in turn. In fact, from 1987 to 2020, the primary forest experienced a reduction of 4040.55 hectares. The secondary forests in turn gained 948.7 hectares. Several hectares of forests converted into agricultural areas were found to have evolved into secondary forests. Agricultural areas increased to around 2954.63 hectares. The savannah class recorded a gain of 72.39 hectares from 1987 to 2020. On the other hand, from 2002 to 2020, a decrease of 62.83 hectares was observed due to the practices of defending savannas in certain areas of the reserve.

Landscape Metrics Analysis at the Luki Biosphere Reserve
In the present study, four landscape metrics, including patch cohesion index, contagion index, aggregation index and landscape division, were computed to examine landscape structure and change of different land use/land cover classes (Figure 8). The primary forest was the dominant land use/land cover class in the Luki Biosphere Reserve landscape, and the change induced in this class was also closely related to eco-environmental issues such as deforestation, land degradation, biodiversity lost, wildlife habitats destruction, and wildlife migration. From 1987 to 2020, there was an increase in the landscape division metric within the primary forest (starting from 0.27 to 0.63), while the patch cohesion index decreased from 99.9 to 93.8. In addition, the aggregation index decreased from 97.6 to 94.5. This contributed to the fragmentation and the reduction of the self-connectivity inside the primary forest. A similar decreasing trend in the patch cohesion index was observed in the secondary forest from 2002 to 2020, with a slight increase of the division metrics. The fallow land and fields and the built-up area experienced an opposite trend for both division and cohesion metrics, when compared with the primary forest. In fact, these two land use/land cover classes experienced an increasing trend in the cohesion index and a decreasing trend in the division index, indicating that they gradually expanded outwards from 1987 to 2020. As for the savannah class, a slight increase in division index was observed. With respect to the aggregation and cohesion indices, the savannah class experienced an increasing trend from 1987 to 2002, while a decreasing trend was observed from 2002 to 2020.

Impact of Village Expansion on Land Use/Land Cover Change
The gradient direction analysis was used for the purposes of visualizing and understanding the direction of land use/land cover change over the time. This was in response to population growth and village expansion (Figures 9 and 10) from the center of the four enclaves (four villages legally established in the reserve), including Tsumba-Kituti, Kiobo, Kisavu and Kibuya. In 1987, land use/land cover classes of Tsumba-Kituti comprised primary forest, secondary forest, savannah, fallow land and fields, and built-up area. The region consisted of several villages whose expansion dramatically contributed to observed changes in the landscape. The center of the region was dominated by the settlements, while the southern part showed a pronounced decrease in the primary forest. The secondary forest experienced a slight decrease across all the regions, while the southern region recorded a high decrease compared to other regions. Fallow land and fields increased in all directions. As for Kisavu village, in 1987 the landscape was dominated by primary forest. However, the expansion of the village induced major changes in the landscape. Over the past 40 years, a huge amount of primary forest was converted into fallow land and fields and secondary forest. Most of these changes were observed in the southern and eastern part of the villages where local communities established farm lands for crop cultivation. For the Kiobo village, primary forest was converted into fallow land and fields, particularly in the southern part of the village. Additionally, some areas of primary forest in the northern part of the village were slightly converted into secondary forest due to the artisanal wood logging that often leads to the degradation of primary forest. Finally, the Kibuya region experienced a slight increase in fallow land and fields in the eastern part due to the expansion of croplands.

Land Use Change of the Luki Biosphere Reserve
Natural ecosystems are of crucial importance for the survival of humanity as they provide ecosystem services, including air purification, climate regulation, food provision, regulation of environmental problems, and conservation of water resources. Population growth and expansion of urban areas are the bases for the transformation of natural landscapes through the conversion of natural ecosystems (forest, savannah, water, wetlands, etc.) into built-up areas, crop and degraded areas. Several studies have shown that urban expansion and population growth are at the root of the problems associated with heat islands [80][81][82], the increase in carbon emissions, and the loss of carbon [83]. In the Democratic Republic of Congo, several studies have revealed that population growth and urban expansion increases the proliferation of cultivation areas as well as the degradation of natural ecosystems [23,[84][85][86][87]. Indeed, many natural land use types have been converted into crop, plantation, and degraded areas due to economic benefits. The loss of natural ecosystems and the intensification of cultivation areas or urban expansion have been the bases of several environmental problems [80,88,89].
In the Luki Biosphere Reserve, most of the villages were developed along the first national road, which connects the two main cities in the region to the city province of Kinshasa, and whose demand for natural resources as well as crop products has been increasing. This continues to increase the risks to the natural resources within the reserve, including deforestation, degradation, and disappearance of species of fauna.
The fragmentation of natural landscapes is another problem caused by urbanization and population growth [75,90,91]. According to Gkyer [92], landscape metrics analysis can help natural resources managers to better comprehend changes that have occurred in the landscape. Numerous landscape metrics, particularly the indices of division, aggregation and cohesion, can reveal the presence of fragmentation of natural ecosystems caused by human activities [92]. Indeed, the landscape division metric can be evaluated over time, and its increase can inform on the fragmentation of the landscape. The patch cohesion index and the aggregation index can also inform on the fragmentation of the landscape when their values decrease. Thus, in the Luki Biosphere Reserve, from 1987 to 2020 the observed increase in the landscape division metric within the primary forest (from 0.27 to 0.63), and the decreases in the patch cohesion index (from 99.9 to 93.8) and the aggregation index (97.6 to 94.5) indicated the fragmentation of the primary forest class leading to the reduction of its self-connectivity. As for the secondary forest class, we observed a decrease trend in the patch cohesion index from 2002 to 2020, while the division metrics slightly increased. In fact, due to human activities, primary and secondary forest have been prone to unsustainable exploitation leading to their fragmentation over the years. An opposite trend was observed in the fallow land and fields and the built-up area, for both division and cohesion metrics, when compared with the primary forest. Indeed, these two classes have experienced an expansion over the years leading to the increase in the cohesion metrics and the decrease in the division metrics.
As shown in Figure 11, the CONTAG index demonstrated the fragmentation of land not only in the entire reserve but also at the zonal level. A high CONTAG value means that the dominant land use patches are well connected within a landscape, while low values indicate that the landscape is more fragmented. CONTAG values decreased from the central area to the transition area. However, landscape fragmentation was observed in all zones of the reserve. The Core zone was characterized by the highest CONTAG values, but they steadily decreased from 1987 to 2020. This means that the landscape of the Core zone was less fragmented compared to the other zones. With the most complex land use/land cover types, the Transition zone's landscape was the most fragmented. From 1987 to 2020, there was an increasing trend of land fragmentation in all zones, as CONTAG continuously decreased in all zones. Fragmentation of land can lead to problems related to land and wildlife management. Fragmentation reduces ecological corridors, endangering animal biodiversity. In Bangkok, the fragmentation of green spaces was found to reduce the density of pollinators, even when floral resources were abundant [93]. Similarly, Chaiyarat et al. [94] reported a decline in bird diversity following the fragmentation of green spaces.

Prediction of Land Use/Land Cover Change
The present research examined LUC change in the Luki Biosphere Reserve over the past 33 years (from 1987 to 2020) with an aim of predicting future LUC changes for the year 2038. Five explanatory variables were used for this purpose, including distance from roads, distance from villages, distance from existing farmland, slope, and elevation [93]. Prediction of the future LUC was obtained by means of the land change modeler plugin implemented in IDRISI software. Thus, based on the LUC change between 1987 and 2002, the transition potential maps were created using the MLPNN algorithm. Afterwards, a simulated map of 2017 was produced using the Markovian chain model. The CROSSTAB module in IDRISI was then used for model validation. Noteworthy was the validation process used in the Pontius et al. [68] method. The results obtained revealed a low percentage of disagreement components (3.6%) and a high percentage of agreement components (96.4%) ( Table 10). Since the validation results provided a high accuracy, the transition potential maps between 2002 and 2017 were then used to predict the 2038 land use/land cover map ( Figure 12) following the same approach. The time interval used to predict the 2038 map was 15 years.  Overall, in 2038 there will still be a decrease in primary forest and savannah of around 899.43 and 223.12 hectares, respectively. The secondary forest will also experience a decrease, losing around 518.88 hectares from 2020 to 2038 due to its proximity to roads and villages, which will lead to its conversion to fallow land and fields. However, built-up area and fallow land and fields will still experience a great expansion by gaining 65.07 and 1576.41 hectares, respectively. The decrease in primary forest can be attributed to the development of unsustainable agriculture, wood energy, expansion of built-up area, and artisanal wood logging. Although savannah decrease is attributed to the expansion of farmland and built-up area, it was noticed that the effort of protecting savannah led to the conversion of this land use/land cover class to forest land.

Conclusions
Since 1987, the Luki Biosphere Reserve and all its zones, including the Core zone, Buffer zone and Transition zone, have undergone dramatic land use/land cover change. The primary forest land decreased throughout the reserve and across all the zones. Secondary forest, savannah, built-up area, and fallow land and fields increased from 1987 to 2020. The Transition zone had the highest land use/land cover change intensity, followed by the Buffer zone and Core zone. The prediction results indicated that there will still be a great decrease in primary forest, while other land use/land cover classes will still face a great expansion, except the secondary forest, which will experience a slight decrease from 2020 to 2038. With respect to landscape analysis, most villages showed clear spatial orientations in the process of village expansion. A strengthening of the connectivity of built-up areas and fallow land was observed throughout the reserve, as well as a phenomenon of fragmentation and disaggregation of other classes of land use outside the built-up area (such as primary forest, secondary forest, and savannah). Over the past few decades, most natural land use types were transformed and replaced by newly grown built-up areas, fallow land and fields. Environmental issues such as environmental degradation, loss of biodiversity, and soil erosion can arise from the loss of different natural land use/land cover classes in the Luki Biosphere Reserve. In addition to this, the fragmentation of the landscape can negatively affect the future of biodiversity in the reserve. Therefore, future studies on the environmental impacts driven by land use/land cover change in the Luki Biosphere Reserve should be initiated to provide reliable information for decision makers for sustainable forest management. Additionally, studies on the dependency of local communities on forest resources should be undertaken in different villages and their surroundings.