Analysing LULC transformations using remote sensing data: insights from a multilayer perceptron neural network approach

ABSTRACT

The study examines the complex dynamics of changes in LULC over three decades, focused on the years 1992, 2002, 2012, and 2022. The research highlights the significance of comprehending these alterations within the framework of environmental and socio-economic consequences. The changes in land use and land cover (LULC) have significant and far-reaching effects on ecosystems, biodiversity, and human livelihoods. This study offers useful information for politicians, conservationists, and urban planners by examining historical patterns and forecasting future changes. The study utilized a Multilayer Perceptron Neural Network (MLP-NN), a well-known machine learning technique that excels at collecting intricate patterns. This model’s design had three layers: input, hidden, and output. The model underwent 10,000 iterations during its training process, and a thorough statistical analysis was conducted to assess the impact of each driving component. The MLP-NN model demonstrated impressive performance, with a skill measure of 0.8724 and an accuracy rate of 89.08%. The accuracy of the LULC estimates for 2022 was verified by comparing them with observed data, ensuring the model’s reliability. Moreover, the presence of evidence likely was found to be a significant factor that had a substantial impact on the accuracy of the model. The study highlights the effectiveness of the MLP-NN model in accurately predicting changes in LULC. The model’s exceptional accuracy and proficiency make it a powerful tool for future LULC forecasts. Identifying the primary causes of model performance and understanding their implications may help to enhance land management strategies, encourage spatial planning, guide accurate decision-making, and facilitate the development of policies that align with sustainable growth and development.

KEYWORDS:

• LULC
• MLP-NN
• Rawalpindi mainland
• land transition
• future simulation

1. Introduction

The study of LULC dynamics, which refers to changes in the Earth’s land cover, has gained significant attention in environmental change. Understanding the dynamics of Land Use and Land Cover (LULC) is vital for gaining insights into environmental shifts, particularly in the context of swift urbanization and ecological changes. Despite considerable research in this area, a substantial knowledge gap exists in comprehending the complex relationship between human activities and LULC variations, especially in developing regions such as Pakistan. This research aims to bridge this gap by proposing an enhanced predictive model that leverages the advantages of current techniques while surpassing their limitations. LULC dynamics across various temporal and spatial scales is vital to modern environmental research. Teck et al. (2023) and Wubie et al. (2016) have shed light on this complex interplay arising from the intricate relationship between human activities and ecological reactions. Solomon and Lukas (2022) have extensively explored this theme, highlighting its multifaceted nature. Human activities, such as expanding agriculture, exploiting forest resources, and developing urban areas, have a significant impact on the land use and cover (Mehmood, Anees, Luo, et al. 2024). Deforestation is one of the notable effects of these changes, which has been extensively studied by Hassen and Assen (2018) and Muhammad et al. (2023). The consequences of these phenomena probably involve land degradation, loss of biodiversity, and a deterioration in environmental services (Shukla et al. 2019; Sohail et al. 2023; Yamanoshita 2019). To gain a comprehensive understanding of the complex dynamics of LULC, scholars have increasingly relied on earth observation satellite data and advanced spatial analytics systems. These tools provide a cost-effective and accurate method for conducting investigations, as demonstrated by the research conducted by Astou Sambou et al. (2023).

Previous studies have relied heavily on earth observation satellite data and spatial analytics to monitor changes in LULC. The accuracy and cost-effectiveness of large-scale analysis methods have been highlighted as crucial (Astou Sambou et al., 2023; Kafy et al., 2021; Mehmood, Anees, Rehman, et al. 2024). However, accurately predicting future LULC dynamics can prove challenging, especially in areas experiencing rapid urbanization and environmental changes. The traditional Markov Chain Model (MCM) is commonly used to estimate LULC changes. However, it needs to consider complex spatial patterns and non-linear interactions between different land cover types, as Mishra and Rai (2016) and Mishra et al. (2018) noted. This study presents a solution to these obstacles by fusing the Artificial Neural Network (ANN) and the Cellular Automata-Markov Chain (CA-Markov) model. This hybrid approach elevates predictive precision by capitalizing on the ANNs’ ability to manage intricate datasets and identify complex patterns in LULC changes, a quality that is absent in conventional models (Ke and Huang 2020; Stateczny et al. 2022). Additionally, integrating the CA-Markov model offers a more comprehensive comprehension of spatial dynamics and temporal variations.

Multilayer Perceptron Neural Networks (MLP-NN) is an innovative and essential tool for studying Land Use and Land Cover (LULC) transformations. MLP-NN is highly effective in handling complex datasets and patterns, particularly relevant due to the significant implications of LULC changes on ecosystems, biodiversity, and socio-economic factors (Badshah et al., 2024). Recent studies, on multiplicative long short-term memory with improved Mayfly optimization for LULC Classification, have highlighted the strengths of neural networks in LULC classification, demonstrating the potential of MLP-NN in this domain (Stateczny et al. 2022). The versatility of MLP-NN, as demonstrated in quality prediction for injection moulding by using MLP-NN supports its applicability in diverse fields, including environmental studies (Ke and Huang 2020). The effectiveness of neural network algorithms in LULC mapping has been demonstrated in employing NNs Algorithm for LULC mapping (Abujayyab and Karas 2020), providing a precedent for their use in LULC research. Similarly, studies such as a multilayer perceptron classifier for monitoring network traffic (Guezzaz et al. 2019) and rainfall runoff modelling by MLP-NN for Lui River Catchment have shown the successful application of MLP-NN in complex domains such as predictive maintenance and malware detection, demonstrating its capability in analysing complex phenomena (Basheer et al. 2022).

The significance of this study cannot be overstated, given the rapid urbanization occurring in areas such as Rawalpindi District, Pakistan, where conventional techniques have fallen short in accurately forecasting changes in land use and land cover (LULC). With its more resilient and all-encompassing modelling methodology, this research makes a valuable contribution to the progress of LULC research, helping to shape sustainable urban planning strategies. We aim to improve the comprehension of historical and projected changes in land use and land cover in Rawalpindi District, Pakistan. Additionally, we present a thorough and effective technique for predicting these changes 1992, 2002, 2012, and 2022) by combining the Artificial Neural Network (ANN) with the Cellular Automata-Markov Chain (CA-Markov) model. This investigation marks a noteworthy advancement in predictive modelling of LULC changes by tackling the intricate dynamics of urbanization and ecological effects within a swiftly evolving global context.

2. Methodology

2.1. Study area

The scope of this study is centred on the Rawalpindi district, situated in the northern region of the Punjab province in Pakistan. The district spans an area of around 5,220 square kilometres and is divided into seven administrative territories: Rawalpindi, Taxila, Gujar Khan, Kalar-Syedan, Kahuta, Kotli Sattian, and Murree as shown in the Figure 1. Rawalpindi district is located at a latitude of 33.4095°N and a longitude of 72.9933°E. It showcases a wide range of physical and climatic features (Khan et al. 2019). The district exhibits a diverse topography, with heights ranging from 297–2274 metres above sea level. This variation in elevation contributes to the presence of different climates throughout the different parts of the district. The northern region, characterized by places like Murree, undergoes harsh winters and moderate summers. Conversely, the southern region, exemplified by locations such as Gujar Khan, experiences scorching summers and mild winters. In June, the highest temperatures range from 25.6°C to 39.4°C, while in January, the lowest temperatures vary from 3.2°C to 16.7°C (Khan et al. 2019). The district experiences an average annual precipitation of 1550 mm, with Murree receiving the maximum amount of rainfall. The topography of the territory consists of hills and gentle slopes that slope towards rivers in the eastern and western directions. These rivers merge with the Soan River in the southern part of the region [DOI, 2012; EIA, 2005]. A substantial area of Rawalpindi district is encompassed by forests, which include both Evergreen Coniferous (CF) and Subtropical Deciduous Forests (STF). The region is inhabited by lush green trees such as Silver Fir (Abies alba), Deodar (Cedrus deodara), serai (Quercus incana), Chir (Pinus roxburghii) and kail (Pinus wallichiana), which possess significant commercial and ecological value (Ansari et al. 2022; N. Muhammad et al. 2023). The majority of these forests are situated in the Bhurban, Masot, and Patriata regions of Murree tehsils, with some also found in specific locations of Kotli Sattian. In Rawalpindi district, you can encounter a variety of tree species, including Shisham (Dalbergia sissoo), Keekar (Acacia arabica), Shareen (Albizia lebbeck), Bari (Ziziphus jujuba), and fruit-bearing trees like Guava and Citrus. Notably, some local farmers have also initiated the cultivation of mini forests featuring Eucalyptus and Symbal trees in recent times (rawalpindi.punjab.gov.pk). The dominant occurrence of Subtropical Deciduous woodlands includes species of phulai (Acacia modesta), kau (Olea cuspidata), and sanatha (Dodonaea viscosa) is observed in the tehsils of Rawalpindi, Gujjar Khan, Kallar Syeddan, and Taxila.

Figure 1. Study area’s geographical location.

2.2. Data set

The MLP-CA-MC model utilizes historical LULC data as response variables to simulate future land-use dynamics. To acquire the historical LULC data, satellite imageries with multi-spectral capabilities, free from cloud cover, were obtained for the years 1992, 2002, 2012, and 2022. The images utilized in this study were obtained from the United States Geological Survey Earth Resources Observation and Science (USGS EROS) facility, which can be accessed at https://earthexplorer.usgs.gov. The purposeful choice of acquisition years, specifically was made based on careful consideration and in response to significant historical occurrences. These years were marked by key events in the country’s history, to compare the impact of both urbanization and afforestation programmes initiated in 2014 under the title Billion Tree Tsunami Afforestation Program (BTFP) on the overall conditions of the natural resources (Kharl and Xie 2017; S. Muhammad 2023). Significantly, the aforementioned images were obtained during the same time frame, thereby reducing the potential influence of seasonal fluctuations (Tariq and Shu 2020).

The MLP-NN model, as described by (Subiyanto, Amarrohman, and Rahmah 2021), allows for the incorporation of many driver variables throughout the modelling procedure. Therefore, this study examined various significant factors that contribute to changes in land usage. The drivers considered in this study included elevation, Distance from the river, aspect, proximity to roads, proximity to urbanization, and Evidence likelihood is shown in Figure 3. The aforementioned datasets were employed for visual image interpretation, as well as the identification and classification of land use. Both the Landsat TM and Landsat 8 OLI-TIRS sensors provide a spatial resolution of 30 metres. Table 1 presents a comprehensive summary of the distinctive features of Landsat data. In the course of image selection, emphasis was placed on the prioritization of imagery that was devoid of clouds and shadows. The inclusion of clouds in the imagery has the potential to considerably undermine the precision of land use classification endeavours (Cerbelaud et al. 2021; Ur Rehman et al. 2021). As a result, it was unfeasible to consistently procure imagery from the same month for the entirety of the study duration to uphold the data’s quality and reliability.

Table 1. Overview of the satellite data utilized.

Satellite	Sensor	Path/Row	Spatial Resolution (M)	Date of Acquisition	Sources
Landsat 5	TM	150/37,36	30	16/06/1992	USGS
Landsat 5	TM	150/37,36	30	18/06/2002	USGS
Landsat 5	TM	150/37,36	30	10/06/2012	USGS
Landsat 8	OLI	150/37,36	30	10/06/2022	USGS

2.2.1. Image pre-processing

The field of image processing involves the analysis and manipulation of digital images using various algorithms and techniques (Dougherty 2020). The utilization of image processing, in combination with classification methodology, plays a crucial role in extracting valuable information from satellite imagery. The ArcGIS Pro 3.1 software was utilized in this study for image processing. In the first step, a technique called layer stacking was employed in the software to combine three bands (specifically, bands 5, 4, and 3 for Landsat 8, and bands 4, 3, and 2 for Landsat 5 TM) into a unified composite layer (Jwan, Mansor, and Khuzaimah 2020). The above-mentioned procedure was carried out on each of the four images that were being examined. Following that, the study area located within the Punjab Province of Pakistan was effectively isolated using the sub-set tool, resulting in the generation of a subset image. The image subset was subsequently transformed to the UTM (Universal Transverse Mercator) zone 42 N and underwent resampling to attain a spatial resolution of 30 metres. Given the nature of our research, which involved the analysis of four separate satellite images captured at different points in time, it was crucial to ensure accurate alignment by establishing Ground Control Points (GCPs). The reference image chosen for image registration was the Landsat-8 scene captured in 2022, specifically path 150 and rows 37, and 36. This reference image was used to align and match the images acquired in 1992, 2002, and 2012. To carry out the registration process, a set of Twenty-Five geographically distributed ground control points (GCPs) were utilized. The alignment accuracy was ensured by employing a second-degree polynomial model (Afwani and Danoedoro 2019; Amalisana and Hernina 2017; Raynaldo, Mukhtar, and Novarino 2020).

Atmospheric correction constitutes a crucial stage within the image processing pipeline. The satellite’s signal may be subject to influence by atmospheric constituents, including gases, solid particles, and liquid particles. The radiance that is detected close to the sensor is commonly known as Top of Atmosphere (TOA) radiance. To address the atmospheric distortions, present in the images (Kabir, Leigh, and Helder 2020; Niraj, Kumar Gupta, and Praise Shukla 2022), we implemented a correction technique inspired by the methodology proposed by (Smith et al. 2021). Figure 2 depicts the comprehensive methodological framework utilized in the present investigation.

Figure 2. A comprehensive methodological procedure of the research.

2.3. Satellite image classification and land use/land cover change assessment in multitemporal analysis

The categorization of satellite images plays a fundamental role within change detection frameworks. The frameworks outlined depend on the use of multitemporal datasets to qualitatively evaluate the temporal dynamics of different phenomena while quantifying the observed alterations. Furthermore, change detection is the act of seeing differences in the state of an item or phenomenon through the observation of that object or phenomenon at several points in time (Alawamy et al., 2020). The satellite imagery in our study was subjected to a classification procedure, leading to the assignment of each image into six separate categories. These categories include (i) Coniferous Forest, (ii) Sub Tropical Forest, (iii) Arable Land, (iv) Barren Land (v) Buildup Area, and (vi) Water. The classification was performed using a supervised-based Maximum Likelihood classifier. Table 2 provides a thorough depiction of the LULC classes.

Table 2. Description of LULC classes.

LULC Classes	Description
Coniferous Forest (CF)	Land with coniferous forest and closed canopy forest.
Sub-Tropical Forest (STF)	Land with scrub forest and open canopy forest
Arable Land (AR)	It includes Crops land, grass land, and fallow land
Barren Land (BA)	It covers Barren land and rocks
Buildup Area (BU)	Urban and settlement area, residential, suburbs and villages, industrial, commercial area and transportation area
Water (WA)	It includes Water bodies, canals, rivers, lakes, ponds, and reservoirs

Following that, we proceeded to assess the spatial variations taking place within each LULC category over thirty (30) years. Additionally, the detection of changes in land-use and land-cover (LULC) was accomplished by utilizing a cross-tabulation module inside the ArcGIS platform, as described in detail in reference (S. W. Wang et al. 2020). The results of this methodology produced a comprehensive land use and land cover change matrix,

which is an important analytical tool for extracting relevant information about the nature and spatial patterns of changes in land use (Christensen and Arsanjani 2020). The change matrix functions as a fundamental framework that aids in the identification and measurement of key categories of changes and their corresponding orientations within the designated study region (Ma et al. 2023). Following this, a thorough examination was conducted to analyse the changes in LULC over specific time intervals. The objective was to identify and understand the patterns of gains and losses that occurred during each of the three distinct periods: the initial period spanning from 1992 to 2002, the subsequent period from 2002 to 2012, and the final period from 2012 to 2022. The change matrices associated with these temporal segments were utilized to determine the changes in the range of land-cover categories, including both gains and losses. To determine the gains for each class, the persistence value was deducted from the total of the respective column. Conversely, losses were calculated by subtracting the persistence value from the total of the corresponding row. The careful and thorough methodology employed guarantees a comprehensive comprehension of the changing land use and land cover patterns through time in the studied geographic area.

2.4. Accuracy assessment of digitally classified land cover maps

Accuracy assessment of digitally classified land cover maps is a crucial step in validating image analysis outcomes. This process involves comparing classification results with ground truth data, typically acquired from reliable geographical references (Stehman and Foody 2019). Ground reference data for our study area, including Rawalpindi city, Gujar Khan, and the Murree sub-division’s rugged terrain, was gathered using Google Earth supplemented by a limited number of GPS points (Koubodana et al., 2019; Stehman and Foody 2019). To enhance classification precision, field trips were conducted in collaboration with local forest department personnel. Data collection was facilitated using Garmin Global Positioning System (GPS) devices, particularly the 12-channel Garmin eTrex 30-Summit mode (Rodríguez-Pérez et al., 2007). These points were evenly distributed across each land cover category for the four-year intervals, as depicted in Figure 4.1.

In the 2022 image, a comprehensive set of 300 testing pixels was constructed, with 50 pixels being sampled from each land cover class (Omo-Irabor 2016), where 180 validation points were physically investigated in the field. These testing pixels were randomly distributed over the whole research area. The testing pixels were then compared with the categorized map. The researchers utilized an error matrix to evaluate the correlation between the categorized map and the reference data. This enabled them to compute various metrics such as producer’s accuracy, user’s accuracy, overall accuracy, and the kappa coefficient (Stehman and Foody 2019). The aforementioned indicators had a collective role in assessing the final land cover maps that were developed for the study. The assessment of accuracy for the historical imagery datasets (specifically, those from 1992, 2002, and 2012) was conducted by accessing the Google Earth Pro archives and ESA WorldCover 10 m 2021 v200 (Zanaga et al. 2022) Landcover historical dataset as the primary source of reference data. The assessment of the accuracy of the classified imagery was performed utilizing the subsequent mathematical expressions, denoted as Equations (1)(1) $\mathit{Overall}\mathit{Accuracy}=\frac{\left(\mathit{Total}\mathit{Number}\mathit{of}\mathit{Correctly}\text{break}\mathit{Classified}\mathit{Pixels}\left(\mathit{Diagonal}\right)\right)}{\left(\mathit{Total}\mathit{Number}\mathit{of}\mathit{Reference}\mathit{Pixels}\right)}\times 100$(1) to (4(4) $\mathit{Kappa}\mathit{Coefficient}\left(\mathit{T}\right)=\frac{\left(\left(\mathit{TS}\times \mathit{TCS}\right)-\sum \left(\mathit{Column}\mathit{Total}\text{break}\times \mathit{Row}\mathit{Total}\right)\right)}{\left(\left(\mathit{TS}-\sum \left(\mathit{Column}\mathit{Total}\times \mathit{Row}\mathit{Total}\right)\right)\right)}\times 100$(4) ). Furthermore, the annual rate of alteration within each distinct category was ascertained through the application of mathematical equations as delineated in Equations (5)(5) $\mathit{Percent}\mathit{of}\mathit{change}=\left(\mathrm{Later}\text{year}-\mathrm{Intial}\text{year} \mathrm{Intial}\text{year}\right)\times 100$(5) and (6(6) $\mathit{Rate}\mathit{of}\mathit{change}\left(\mathit{ha}/\mathit{yr}\right)=\left(\frac{\mathrm{Later}\text{year}}{\mathrm{Intial}\text{year} \mathrm{Time} \mathrm{interval}}\right)$(6) ).

(1) $\mathit{Overall}\mathit{Accuracy}=\frac{\left(\mathit{Total}\mathit{Number}\mathit{of}\mathit{Correctly}\text{break}\mathit{Classified}\mathit{Pixels}\left(\mathit{Diagonal}\right)\right)}{\left(\mathit{Total}\mathit{Number}\mathit{of}\mathit{Reference}\mathit{Pixels}\right)}\times 100$(1)

(2) $\mathit{Use}{r}^{\prime }\mathit{s}\mathit{Accuracy}=\frac{\left(\mathit{Number}\mathit{of}\mathit{Correctly}\mathit{Classifed}\text{breakPixels}\mathit{in}\mathit{each}\mathit{Catagory}\right))}{\left(\mathit{Total}\mathit{Number}\mathit{of}\mathit{Classified}\mathit{Pixels}\mathit{in}\mathit{that}\mathit{Catagory}\text{break}\left(\mathit{RowTotal}\right)\right)}\times 100$(2)

(3) $\mathit{Produce}{r}^{\prime }\mathit{s}\mathit{Accuracy}=\frac{\left(\mathit{Number}\mathit{of}\mathit{Correctly}\mathit{Classifed}\text{breakPixels}\mathit{in}\mathit{eachCatagory}\right)}{\left(\mathit{Total}\mathit{Number}\mathit{of}\mathit{Classified}\mathit{Pixels}\mathit{in}\mathit{that}\mathit{Catagory}\text{break}\left(\mathit{Column}\mathit{Total}\right)\right)}\times 100$(3)

(4) $\mathit{Kappa}\mathit{Coefficient}\left(\mathit{T}\right)=\frac{\left(\left(\mathit{TS}\times \mathit{TCS}\right)-\sum \left(\mathit{Column}\mathit{Total}\text{break}\times \mathit{Row}\mathit{Total}\right)\right)}{\left(\left(\mathit{TS}-\sum \left(\mathit{Column}\mathit{Total}\times \mathit{Row}\mathit{Total}\right)\right)\right)}\times 100$(4)

(5) $\mathit{Percent}\mathit{of}\mathit{change}=\left(\mathrm{Later}\text{year}-\mathrm{Intial}\text{year} \mathrm{Intial}\text{year}\right)\times 100$(5)

(6) $\mathit{Rate}\mathit{of}\mathit{change}\left(\mathit{ha}/\mathit{yr}\right)=\left(\frac{\mathrm{Later}\text{year}}{\mathrm{Intial}\text{year} \mathrm{Time} \mathrm{interval}}\right)$(6)

2.5. A robust approach for LULC change prediction: Markov-CA

The Markov-Cellular Automata (Markov-CA) approach has become increasingly prominent in the field of land use change prediction. This is primarily attributed to its strong performance in comparison to other methodologies, as evidenced by previous studies conducted by (Jadawala, Shukla, and Tiwari 2021; Koko et al. 2020; Moradi, Seyed Kaboli, and Lashkarara 2020; Naimur Rahman et al., 2023). The Markov-CA model has been observed to possess the ability to predict bidirectional transitions between LULC classes (Vinayak, Lee, and Gedem 2021; S. W. Wang et al. 2020). The anticipation of future LULC changes using a Markov-Cellular Automata (CA) model involves a systematic procedure consisting of three fundamental stages:

2.5.1. Initial application of Markov chain analysis

The first phase entails employing Markov chain analysis on historical LULC maps from different periods, specifically those from 1992, 2002, 2012, and 2022. The main aim of this study is to generate transition matrices that encompass the probabilistic data regulating LULC changes over a specific period.

2.5.2. Generation of transition potential maps

The computation of transition potential maps for LULC is conducted as an integral component of the second stage. The maps presented here serve as a visual depiction of the spatial likelihood for changes between distinct LULC categories (Aguejdad 2021). The execution of the Cellular Automata Model: The Cellular Automata (CA) model is subsequently employed to analyse the acquired transition matrices and transition potential maps in the concluding phase. The cellular automaton (CA) model incorporates the probabilistic transition data along with the geographical attributes outlined in the transition potential maps to forecast the forthcoming spatial arrangement of LULC (Mansour, Al-Belushi, and Al-Awadhi 2020). The utilization of the Markov-CA framework as a methodological approach provides a robust and empirically validated mechanism for forecasting changes in LULC. This methodology enables a holistic comprehension of the temporal dynamics of landscapes.

2.5.3. Estimation of transition probabilities

Markov model, LULC change prediction framework was applied as indicated in Figure 2. The key element of this methodology is the projection of the distribution of each LULC class, based on the transition probability pij. This probability represents the possibility of moving from one LULC class i to another LULC class j during a specific time interval from t to t + 1 (Lee, Judge, and Takayama 1965; S. Wang and Zheng 2023a). The calculation of p_ij is formally defined as follows:

(7) $P_{\mathit{ij}}=\frac{{\mathit{n}}_{\mathit{ij}}}{n_i}$(7)

It is essential to underscore that, in our calculations, P_ij is bounded within the range [0,1], thereby ensuring that the sum of transition probabilities for all potential transitions from class i equals 1:

(8) $\sum _{\mathit{j}=1}^k{P}_{\mathit{ij}}=1$(8)

2.5.4. Projection of LULC classes at time t + 1

To estimate the distribution of each LULC class at time t + 1 (M_t + 1), we utilized the transition probability matrix (P) in combination with the LULC distribution at the starting time t (Mt). The aforementioned relationship is denoted as:

(9) $\mathit{P}\times M_t=M_{\mathit{t}+1}$(9)

The procedure entails the multiplication of the transition probability matrix (P) with the LULC distribution matrix at time t (Mt), resulting in the projected LULC distribution at time t + 1 (M_{t + 1}). The foundation of our methodology relies on the utilization of land use/land cover maps for the years 1992, 2002, 2012, and 2022, which were developed through the classification procedure. The maps were utilized as the basis for generating transition probability matrices, which represented the changes between LULC classes across the periods of 1992 to 2002, 2002 to 2012, and 2012 to 2022. It is crucial to acknowledge that a first-order Markov model was utilized, which asserts that the probability distribution of LULC classes in the next stage is only determined by the current distribution (Asif et al. 2023; Azubike, Kurkalova, and Mulrooney 2019; S. Wang and Zheng 2023b). The transition matrices were subsequently utilized in the subsequent stages of our land use and land cover change prediction approach.

2.5.5. Spatial suitability analysis and MLP-NN performance

The spatial suitability and assessment of the performance of the MLPNNs are essential processes in comprehending the dependability and efficacy of the model within the framework of simulating LULC change. The MLP-NN utilizes the Backpropagation (BP) learning technique, in which the forward and backward passes are iteratively performed until the network successfully acquires knowledge of the unique attributes associated with each class (Hemeida et al. 2020; P. S. Kumar et al. 2020). The importance of selecting appropriate variables cannot be overstated, as it significantly affects the accuracy of the model’s learning process (Akdeniz, Serdaroglu Sag, and Inam 2023; Debnath et al. 2023). Therefore, it is crucial to perform sensitivity analysis to evaluate how the model responds to changes in these variables before conducting any simulations related to LULC changes. The MLP classifier, which is built into the TerrSet image processing toolset, was utilized to assess the parameters and performance of the model.

The sensitivity analysis entailed a methodical procedure. The model was initially trained using all the explanatory factors, as depicted in Figure 2, along with the 1992 LULC data. Following this, a series of skill tests were undertaken to assess the comparative impact of each particular explanatory variable. The achievement of this task was carried out by deliberately keeping the inputs from particular variables unchanged (Singh et al. 2022). The disparity in expertise, represented as ‘S’ and calculated using Equation (9)(9) $\mathit{P}\times M_t=M_{\mathit{t}+1}$(9) , offers valuable insights into the importance of each variable.

Figure 3. Input layer nodes: 1. PAspect, 2. Proximity to urbanization, 3. Proximity to roads 4. Elevation, 5. Distance from the river, and 6. Evidence likelihood.

Three separate sensitivity analyses were conducted in this study. The first analysis involved keeping a single independent variable constant, while the second analysis involved holding all independent variables constant except for one. The third analysis utilized a backward stepwise constant forcing mechanism. The skill statistic, denoted as ‘S’, exhibits a range of values between −1 and 1. A skill value of 1 represents impeccable predicting ability, while a value of −1 indicates performance that is inferior to random chance (Dolui and Sarkar 2023). Furthermore, a skill value of 0 predicts results that are comparable to those expected by random chance. It is important to acknowledge that achieving a model accuracy of 80% or more is considered desirable to evaluate the training outcomes as adequate (Mansour et al. 2023).

(10) $\mathrm{S}=\mathit{\hspace{0.17em}}\left(\mathrm{A}-\mathrm{E}\left(\mathrm{A}\right)\mathit{\hspace{0.17em}}\right)\mathit{\hspace{0.17em}}\left(\mathit{\hspace{0.17em}}1\mathit{\hspace{0.17em}}-\mathrm{E}\left(\mathrm{A}\right)\mathit{\hspace{0.17em}}\right)$(10)

Where ‘A’ represents the quantified accuracy based on measurements, and ‘E(A)’ denotes the anticipated accuracy, a value derived by considering the number of transitions within the sub-model (‘T’) and the count of persistence classes (‘P’) as defined in Equation (10)(10) $\mathrm{S}=\mathit{\hspace{0.17em}}\left(\mathrm{A}-\mathrm{E}\left(\mathrm{A}\right)\mathit{\hspace{0.17em}}\right)\mathit{\hspace{0.17em}}\left(\mathit{\hspace{0.17em}}1\mathit{\hspace{0.17em}}-\mathrm{E}\left(\mathrm{A}\right)\mathit{\hspace{0.17em}}\right)$(10) :

(11) $\mathrm{E}\left(\mathrm{A}\right)\mathit{\hspace{0.17em}}=\mathit{\hspace{0.17em}}1\mathit{\hspace{0.17em}}\left(\mathrm{T}+\mathrm{P}\right)\mathit{\hspace{0.17em}}$(11)

The performance and prediction efficiency of MLP-NN were evaluated by employing the Root Mean Square Error (RMSE). Root Mean Square Error (RMSE) is a crucial metric used to measure the extent of discrepancy between projected values and their corresponding actual values. It is important to acknowledge that a decrease in the root mean square error (RMSE) indicates an increase in the level of accuracy in predictions (Shahhosseini et al., 2021).

(12) $\mathit{RMSE}=\sqrt{\frac{{\sum }_{\mathit{i}=1}^N\left(\stackrel{\approx }{\text{relbar}\phantom{{\mathit{ }}_{\text{vboxto}.5\mathit{ex}\text{vss}}}}\mathrm{i}-\mathrm{u}\right)}{\mathit{N}}}$(12)

where ‘ûi’ represents the modelled value for sample ‘i’, ‘u i’ corresponds to the observed data for sample ‘i’, and ‘N’ signifies the total number of samples encompassed within the analysis.

2.6. Transition potential modeling and LULC change projection

The Land Change Model’s (LCM) transition sub-model initially incorporated driver factors that had undergone skill testing. Following this, the MLPNNs was employed to produce probable transition maps using the dependent variables, namely the T1 and T2 imageries. In this phase, the transition potential images successfully depicted the influence of the driver variables, consistent with previous studies (Leta, Demissie, and Tränckner 2021), thereby offering valuable insights into the applicability of cell transformation for particular land cover categories (Dolui and Sarkar 2023; Yaghoobi et al. 2022). Subsequently, the generation of both hard and soft predictions for the LULC changes in the year 2021 (referred to as T’3) was carried out following the procedure principles illustrated in Figure 2. The aforementioned predictions were utilized as visual representations to validate the MLP-CA-MC model.

2.7. Validation of MLP-CA-MC model output

The validation process of the MLP-CA-MC model’s output was conducted to evaluate the degree of agreement and disagreement between the observed satellite-derived LULC map of 2022 (T3) and the simulated LULC map (T’3) as shown in Figure S1. The validation method holds significant importance in assessing the accuracy and credibility of the MLP-CA-MC model’s forecasts for the upcoming years 2032, 2042, and 2050, as evidenced by prior research (Deng and Quan 2022; Gaur et al. 2020; Xu, Gao, and Coco 2019). This study utilized two independent validation methodologies, encompassing both hard and soft predictions. The validation process was performed utilizing the VALIDATE and ROC modules that are accessible within the TerrSet software. The Receiver Operating Characteristic (ROC) module was utilized to compute the area under the receiver operating characteristic curve (AUC), often known as ROC statistics. The T’3 soft prediction map was employed as a comparative map in this procedure (Giglioni et al. 2021; Parsons 2020). The AUC values span a range of 0 to 1, with 0 representing a test that is entirely wrong and 1 indicating a completely accurate test (Baig et al. 2022; V. Kumar and Agrawal 2023a). In the meanwhile, the VALIDATE module performed computations on the kappa index statistics, utilizing the T’3 hard prediction map as the reference for comparison.

The kappa statistics obtained in this study encompassed various aspects, including kappa for grid cell-level location (K_location), kappa for cases with missing information (referred to as K_no), kappa for stratum-level location (referred to as K_{locationStrata}), and the standard kappa (referred to as K_standard) (Idrissou et al. 2022; Näschen et al. 2019). The aforementioned validation techniques jointly offer a complete evaluation of the model’s performance and its capacity to reliably forecast LULC changes for the specified future years. Previous studies (Basu et al. 2021; Idrissou et al. 2022; Toma, Belete, and Ulsido 2023) have indicated that an Area Under the Curve (AUC) and Kappa value that is resilient and acceptable is often expected to be 80% or higher. In the framework of predicting future outcomes, the LULC maps for the years 1992, 2012, and 2022 were utilized as dependent variables to model and predict the LULC maps for the years 2032, 2042, and 2052. The projection was derived from a business-as-usual (BAU) scenario, which postulates the continuation of current trends without substantial interventions or policy alterations. It is worth mentioning that a reservoir characteristic was identified as a limiting factor in the modelling procedure, and zoning methodologies were implemented in the phase of allocating changes. This facilitated the multiplication of the transition potentials linked to each land cover transition with the limitations map. The utilization of a constraint-based strategy was crucial in influencing the anticipated changes in LULC.

3. Results

3.1. Assessing changes in LULC extent and magnitude over time

The study examined the magnitude and scope of changes in LULC during a designated period. The total area of coniferous forests in the years 1992, 2002, 2012, and 2022 was determined to be 413.66, 336.40, 276.20, and 178.64 km², respectively. The subtropical forest exhibited a cohabitation of mixed tree species, spanning a total land area of 1213.93 km² in the year 1992. However, this area experienced a reduction and declined to 914.28 km² by the year 2022. The majority of the population in the study area primarily resided in both urban and rural regions and significantly depended on agricultural operations to fulfil their household needs. The agricultural land area in 1992 was recorded as 2939.07 km², whereas in 2022 it expanded to 3515.29 km². The build-up class experienced a significant increase in size, expanding from 82.87 km² to 344.10 km² over 30 years. In contrast, the barren land category exhibited the largest extent of ground coverage, encompassing 536.59 km² in 1992 and then diminishing to 240.47 km² in 2022. Diverse land cover classes showed noticeable changes. The trend of increase and decrease in land cover classes in percentage are clearly illustrated in Figure 5. Coniferous forests (CF), which made up 7.92% of the total area in 1992, predominated the landscape, followed by subtropical forests (STF), which made up 23.25% of the total area, agricultural land (AL), that filled up 56.30% of the total area, barren land (BL), which made up 10.28% of the total area, built-up areas (BU) cover 1.59% of the total area, and water areas (WA) at 0.65%. Agriculture area has grown significantly over time, growing by 40% by 2022. On the other hand, during the same period, coverage of forests (CF and STF) decreased by 11.2% and 10.20%, respectively. Conversely, the built-up area (BU) displayed a contrasting pattern, indicating rapid urbanization. Between 1992 and 2002, the BU expanded at a decadal rate of 2.36%, reflecting the growth of urban areas. This expansion accelerated significantly from 2002 to 2012, with a decadal rate of 11.67%, and continued at a similar pace from 2012 to 2022, growing by 12.08% over the decade. This substantial increase in built-up areas signifies the ongoing urban development and its associated impacts on land cover. The overall accuracy rate of classified imageries for the years 1992, 2002, 2012, and 2022 are 0.81, 0.82, 0.91, and 0.91% with a Kappa coefficient of 0.90, 0.89,0.54, and 0.25 respectively as shown in Table 3. The annual variation in the stated categories was calculated using the widely accepted methodology that considers both spatial and temporal data. The forest-based classifications of coniferous forests and subtropical forests both showed a decreasing trend, with yearly rates of change of 0.58% and 1.35%, respectively, throughout the designated research region. To obtain a thorough understanding of the areas examined and their respective proportions Table S2 and Figure 4 presented the findings of this study.

Figure 4. LULC classes maps and graph of (a) 1992, (b) 2002, (c) 2012, (d) 2022.

Table 3. Classification accuracy assessment results for 1992, 2002, 2012 and 2022.

Year	CF		STF		AL		BL		BU		WA		OA	KA
	PA	UA	PA	UA	PA	UA	PA	UA	PA	UA	PA	UA
1992	0.97	0.90	0.89	0.93	0.83	0.95	0.95	0.90	0.97	0.90	1.00	1.00	0.93	0.92
2002	0.90	0.89	0.83	0.82	0.76	0.97	0.88	0.82	0.98	0.77	0.97	1.00	0.88	0.85
2012	0.94	0.90	0.83	0.86	0.84	0.92	0.93	0.84	0.94	0.88	0.94	1.00	0.90	0.88
2022	0.95	0.95	0.84	0.90	0.84	0.88	0.93	0.95	0.95	0.88	0.97	0.90	0.91	0.90

Figure 5. LULC Classes trend from 1992–2022.

3.2. Analysis of the transition probability matrix

This study presents an analysis of the transition probability matrix within the framework of a Markov chain model. This analysis specifically examines the transition probabilities that were observed throughout three unique periods: 1992–2002, 2002–2012, and 2012–2022, as outlined in Table S3. The matrix presented in Table S3 utilizes a row-column format, where the rows represent past classifications of LULC, while the columns indicate the expected or projected LULC categories. In each of the transition matrices, the diagonal elements represent the probability of a specific land cover class remaining unchanged during the study period (Ntakirutimana and Vansarochana 2021). The utilization of this transition probability matrix is crucial for comprehending the dynamics related to the shifts between different LULC categories. Additionally, it allows for the quantification of the probabilities connected with a cell, belonging to a specific LULC category, changing to any other category within the designated time frame.

3.3. Analysis of the transition area matrix

In this study, the historical LULC data from 1992 to 2022 were utilized to create the transition area matrixes. These parameters are crucial for analysing the future dynamics of LULC. The transition area data offer useful insights into the spatial changes that occur between different LULC categories, namely the shift from Forest Lands to Agriculture Land. The matrix in Table 4 displays the transition area file, which provides a quantitative representation of the total number of pixels or regions within cells that are expected to experience changes in the following period. The transition area matrix was generated for each pair of LULC areas, such as 1992–2002, 2002–2012, and 2012–2022, to illustrate the anticipated changes for each LULC category. Changes in LULC have a substantial impact on the environment and the management of land resources. In this particular study, the transition area files corresponding to the period from 2002 to 2022 were utilized to predict the distribution of LULC for the year 2022. The LULC projection was afterwards utilized as a validation dataset to evaluate the correctness of the Cellular Automata Markov (CA_Markov) model in comparison to the observed 2022 LULC. In contrast, the transition area dataset that corresponds to the time frame of 1992 – 2022 was utilized to forecast the distribution of LULC for the years 2032, 2042, and 2052.

Table 4. Markov transition area file matrix depicting transitions among diverse LULC categories during 1992–2002, 2002–2012, 2012–2022 and 1992–2022.

LULC Class Type		1992/2002
	AL	BL	BU	CF	STF	WA	GT
AL	2461.82	269.39	28.71	17.70	123.96	0.78	2902.36
BL	336.85	148.51	28.97	4.39	52.89	2.50	574.12
BU	15.31	7.08	39.15	0.52	11.66	0.30	74.01
CF	4.48	0.25	0.32	266.24	142.59	0.01	413.89
STF	330.55	59.31	13.62	47.63	766.12	0.56	1217.78
WA	0.91	5.71	0.28	0.02	1.67	25.18	33.77
GT	3149.92	490.25	111.06	336.50	1098.88	29.33	5218.31
	2002/2012
AL	2733.63	108.54	113.68	1.36	192.16	2.72	3152.09
BL	384.48	52.52	29.60	0.04	19.02	4.85	490.49
BU	38.34	4.67	63.20	0.01	4.22	0.41	110.85
CF	30.64	0.03	0.96	178.49	126.35	0.00	336.47
STF	276.06	15.27	27.03	98.49	681.39	1.32	1099.55
WA	0.90	5.54	0.24	0.00	0.07	22.47	29.22
GT	3464.06	186.56	234.70	278.38	1023.22	31.77	5218.69
	2012/2022
AL	2924.36	171.33	195.11	1.24	165.71	4.54	3462.28
BL	88.20	52.15	33.35	0.02	5.18	7.60	186.50
BU	56.30	18.28	152.51	0.11	6.51	0.91	234.62
CF	16.94	1.53	0.18	139.93	119.54	0.30	278.42
STF	340.21	11.02	13.87	37.00	616.83	3.36	1022.30
WA	10.50	4.88	1.96	0.01	1.65	12.83	31.81
GT	3436.50	259.18	396.97	178.31	915.42	29.55	5218.31
	1992/2022
AL	2478.85	141.70	193.51	0.30	84.80	4.20	2903.35
BL	358.34	81.52	96.46	0.90	32.92	4.35	574.50
BU	14.37	1.69	53.57	0.07	3.70	0.67	74.06
CF	52.15	1.53	1.33	163.55	194.85	0.64	414.05
STF	522.89	27.85	50.65	13.51	598.34	5.38	1218.62
WA	11.18	5.09	1.74	0.00	1.53	14.40	33.94
GT	3437.79	259.38	397.26	178.33	916.15	29.62	5220.54

3.4. Fuzzy standardization and suitability map generation

The CA_Markov modelling approach involves a judiciously guided procedure for the spatial allocation of anticipated LULC patterns, leveraging categorical suitability maps as a foundational component (EL-Kawy et al. 2019; Ntakirutimana and Vansarochana 2021). As illustrated in Figure 3, these suitability maps encapsulate both physical and socio-economic factors, each of which contributes to the overall assessment of suitability for specific land use categories. These generated suitability maps serve as essential inputs in the formulation of the transition suitability image, facilitating the integrated analysis process within the CA Markov model.

The significance of the transition potential images as a useful instrument for managing the spatial arrangement of LULC (Hasan et al. 2020). As a result, we applied the fuzzy standardization tool (Ruda, Kolejka, and Batelková 2017), to convert each contributing element in our study into a relative scale. This scale ranged from 0 to 255. In this modified scale, a pixel with a value of 0 represents the lowest level of suitability, indicating unsuitability. Conversely, a pixel with a value of 255 indicates the highest level of appropriateness for a particular LULC category. The socio-economic factors are considered to be the primary drivers of land use and land cover change. Nevertheless, the study focused solely on physical parameters, perhaps neglecting the intricate dynamics that may be involved. This study suggests that the exclusive reliance on physical parameters may prove inadequate in accurately predicting future land use and land cover dynamics inside the CA_Markov model. One possible explanation might be attributed to the impact of socio-economic factors, such as population density, which can play a significant role in driving LULC changes. Regions characterized by a significant concentration of inhabitants can potentially contribute to the depletion of water bodies, forested areas, and shrublands. Consequently, the exclusion of socio-economic factors may result in less precise forecasts of LULC consequences. To improve the accuracy of future LULC simulations, it is recommended to include a blend of physical and socio-economic data, along with other relevant characteristics, as suggested by (de Sousa Oliveira et al. 2022; Taiwo et al. 2023). The implementation of this comprehensive methodology has the potential to enhance the accuracy and reliability of our comprehension and forecasting of LULC dynamics. Furthermore, empirical studies have demonstrated that regions characterized by a lack of access roads tend to exhibit lower levels of human interference and disturbance, mostly due to the absence of easily accessible entry sites.

3.5. Gain and loss

The Land Cover Change Analysis Module was employed to ease the acquisition of gains and losses throughout the selected research period. The results are graphically represented in Figure 6. The most notable expansion in land cover was mostly found in the categories of arable land and urban land. Simultaneously, there was a significant decrease in the categories related to forests across this research.

Figure 6. Gains and losses graph between 1992–2002, 2022–2012 and 2012–2022.

In the years 1992 and 2002, there was a notable increase in the agricultural domain, with an expansion of 722.46 km². Additionally, urban land experienced a growth of 68.04 km² during the same period (Figure 7). The expansion observed between 2002 and 2012 can be attributed mostly to the land types of urban land, arable land, and barren land, as depicted in Figure 7. Between the years 2002 and 2012, there was an observed increase in arable land by 742.01 km², while urban land expanded by 163.04 km². This expansion was mostly at the cost of coniferous and sub-tropical forest areas. In contrast, an area of 108.40 square kilometres underwent a shift from being arable to becoming barren land, while concurrently obtaining an additional 306.7 square kilometres from categories associated with forests. It is important to notice that the extension of barren land was significantly influenced by arable land and build-up. The period from 2002 to 2012, illustrates a significant conversion of 159.4 km² of coniferous forest and 432.5 km² of subtropical forest into different land use categories.

Figure 7. Contributions to net change (Sq. km) in coniferous forest (CF), subtropical forest (STF), arable land (AL), barren land (BL), buildup (BU), and water (WA) between 1992–2002, 2002–2012, 2012–2022.

Between the years 2012 and 2022, there was a significant shift in land use patterns seen within the geographic area under investigation. The transformative process resulted in agricultural land experiencing a significant increase of 504.4 square kilometres, making it the primary benefactor. The expansion described is of particular significance due to its utilization of pre-existing land areas that were formerly classified as subtropical forest (340.2 km²), coniferous forest (16.9 km²), built-up areas (56.3 km²), and barren land (88.2 km²), as depicted in Figure 7. Concurrently, the observed change in land use was accompanied by a proportional decrease in forest-covered land, amounting to a total of 534.26 square kilometres. The decline in forested coverage can mostly be attributable to the increasing agricultural activities that occurred throughout the specified period. Furthermore, it is worth noting that the growth of Urbanization throughout the designated timeframe resulted in a net increase of 209 square kilometres. The process of expansion resulted in a reduction of the coniferous forest by a minimum area of 0.18 square kilometres, subtropical forest by a significant area of 13.87 square kilometres, arable land by a notable expanse of 195.1 square kilometres, and barren land by a combined area of 33.3 square kilometres. It is important to acknowledge that the significant expansion of the build-up was primarily driven by the transformation of arable and barren land parcels into different land use categories.

Over the entire period of investigation, the build-up class displayed noticeable transformations primarily attributed to modifications in land cover classifications, including barren land, and coniferous and subtropical forests. The aforementioned transitions can be concisely summarized in Table 4 and visually represented in Figure 9. There was a significant increase in the extent of agricultural land, which expanded from 2939.07 km² in 1992 to 3515.3 km² in 2022. In contrast, the land devoid of vegetation experienced a steady decrease in size, decreasing from 536.5 km² in 1992 to 476.2 km² in 2002, and continuing to reduce from 2012 to 2022. The drop observed throughout the period from 1992 to 2012 can be attributed to the transformation of unproductive land into cultivable land and subsequent development. Furthermore, there was an observed increase in the built-up area, arable land, and water bodies. In contrast, a significant decrease was seen in the areas of coniferous forest, declining from 413.6 to 336.4 km², and sub-tropical forest, declining from 1213.9 to 1098.8 km², with a net change of 77.4 and 120.7 km² respectively between the years 1992 and 2002. This loss persisted from 2002 to 2022 as well, as depicted in Figure 9.

The net changes in land cover during the period from 1992 to 2002 exhibited significant variations, characterized by a positive net expansion of 269 square kilometres in arable land and an additional 32.07 square kilometres in built-up areas. A comparable trend was observed during the next decade spanning from 2002 to 2012, wherein the net alteration equated to 319.81 square kilometres for arable land and 116.88 square kilometres for built-up areas. During the last decade of investigation, a significant increase was found in the BU land cover categories, as depicted in Figure 8 and elaborated upon in Table S2. The transformation maps depicted in Figure 9 were created to clarify the alterations that occurred during the periods of 1992–2002, 2002–2012, and 2012–2022. These maps were made by employing image difference operations within the software ArcGIS Pro. The maps shown in this study illustrate the distribution of unaltered values, which represent the percentage of pixels that have remained consistent over time. Additionally, the maps also depict the decreasing class, which indicates changes towards water bodies, and coniferous and subtropical forests. Moreover, a rise in class signifies the extension of agricultural activities, the utilization of unproductive land, and the development of urban infrastructure, frequently entailing the conversion of coniferous and subtropical forests.

Figure 8. Net change graph between 1992–2002, 2022–2012, and 2003–2021.

Figure 9. LULC transformation maps of the study area between 1992–2002, 2002–2012, and 2012–2022.

3.6. Multilayer Perceptron Neural Network (MLP-NN) architecture and skill evaluation

The present study utilized a MLP-NN, a well-known machine learning technology recognized for its ability to effectively capture complex patterns and behaviours (Kussul et al. 2017a; Raj and Sharma 2022). The architecture utilized in this experiment consisted of a MLPNNs with two layers: an input layer, a hidden layer, and an output layer. The specific configuration of these layers is shown in Table 5. Every layer in the system was composed of interconnected nodes, which are also known as neurons, with specific weights assigned to each connection. The input layer consisted of seven nodes representing driving variables, the hidden layer contained eight nodes, and the output layer consisted of nine nodes. A comprehensive statistical study was undertaken to evaluate the significance and influence of each driving variable on the performance of the MLP-NN model. The present investigation entailed examining the driver factors concerning the MLP-NN skill measure over 10,000 iterations (Eastman 2020; Kussul et al. 2017b). In this approach, three unique sensitivity analysis techniques were utilized, as discussed in the following sections.

Table 5. Parameter and performance.

Parameters	Performance
Input layer neurons	6
Hidden layer neurons	7
Output layer neurons	9
Requested samples per class	1500
Final learning rate	0.0001
Momentum factor	0.5
Sigmoid constant	1
Acceptable RMS	0.01
Iterations	10000
Training RMS	0.2347
Testing RMS	0.236
Accuracy rate	89.08%
Skill measure	0.8724

The main aim of training the MLP-NN was to ascertain the most favourable weights for the connections linking the input and hidden layers, as well as the hidden and output layers. The employment of optimized weights allows the neural network to proficiently categorize unfamiliar pixels by leveraging the acquired patterns and behaviours from the training dataset. The MLP-NN model had a notable performance, achieving an overall skill of 0.8724 when accounting for all variables. Additionally, the model exhibited an accuracy rate of 89.08%. The performance measures observed in this study are above the established threshold, providing evidence that the trained model is appropriate for making future projections. This is illustrated in Table 5.

3.6.1. Relative influence of independent variables on model accuracy

In our study, we performed an analysis to determine the comparative impact of different independent variables. This was achieved by maintaining the constancy of each variable while assessing the model. The results, as outlined in Table S4, demonstrated a significant trend. In particular, when the evidence likelihood was held constant, the model had a significantly reduced degree of accuracy, measuring only 33.1%. The outcomes of our study indicate that evidence likelihood is the most influential characteristic in our research environment. This aligns with previous research conducted in a similar geographical context (Leta, Demissie, and Tränckner 2021), which likewise underlined the significance of evidence probability as a significant variable.

In contrast, when we maintained a constant aspect, we noticed a rather minimal effect on the performance of the model. Based on our investigation, it can be concluded that this aspect had the least impact on the parameters examined. Nevertheless, it is imperative to recognize the distinctive biophysical characteristics of our research, which may have a substantial influence on environmental dynamics (Simwanda, Murayama, and Ranagalage 2020)

3.6.2. Isolating individual independent variables by holding all but one constant”

To examine the specific effects of individual independent variables, a research approach was utilized that involved keeping all variables constant except for one. The present study employed this methodology to gain a comprehensive understanding of the unique characteristics associated with each variable being examined. Based on the information shown in Table S5, it is apparent that there was a significant variation in accuracy and skill metrics observed in consecutive tries, except for the initial try where all variables were unlimited. The presence of interaction effects and intercorrelations among the input variables, as previously discussed by (Eastman 2020), is strongly indicated by this significant variation in outcomes.

3.6.3. Backwards stepwise variable selection

The training strategy commenced by incorporating all variables, thereafter implementing a methodical approach of isolating and maintaining each variable at a constant level to evaluate its influence on the model’s proficiency. The objective of this approach was to determine pairs of variables whose joint exclusion had minimal impact on the overall performance of the model. As a result, it was noted that there were instances where the model’s proficiency displayed a slight enhancement upon the removal of specific variables. Nevertheless, it is important to highlight that in the present investigation, all variables, as indicated in Table S6, were included in the analysis without any exclusions. The rationale behind this decision was based on the lack of significant differences in the model’s performance found throughout both the initial and subsequent stages of the investigation, as illustrated in Table S6.

3.7. Validation and performance assessment of the MLP-CA-MC model for simulating LULC patterns in 2022

To validate the MLP-CA-MC model, it was initially applied to simulate the LULC patterns in the year 2022, referred to as T’3, encompassing both soft and hard categories. The simulation utilized LULC maps from 1992 (T1) and 2012 (T2) as input data. To assess the precision of the simulations, Kappa statistics were computed, considering both the quantitative and spatial attributes and subsequently compared to a reference map representing the year 2022. The obtained statistical measures demonstrated positive outcomes, including an 83% Kappa location (K location) coefficient, an 83% Kappa location strata (K locationStrata) coefficient, and a Kappa standard (K standard) coefficient of 81%. It is important to highlight that all the Kappa index values mentioned in the study (de Sousa Oliveira et al. 2022; Eastman 2020; Hasan et al. 2020) are above the established threshold of 80%. This indicates a significant level of agreement between the simulated LULC maps and the actual ones.

Quantity disagreement refers to situations in which the number of cells belonging to a specific category in T’3 is not the same as the number of cells in T3. In a study conducted by (M. M. Gharaibeh 2021; A. Gharaibeh et al. 2020), it was observed that there was a dispute in the spatial position of the same category between T’3 and T3. The bar graph in Figure 10 illustrates the agreement and disagreement components, as well as the overall proportion of valid predictions (91.3%) for the MLP-CA-MC model.

In addition to employing Kappa statistics, the Receiver Operating Characteristic (ROC) metric was employed to evaluate the model’s predictive capacity in determining locations based on the distribution of a Boolean variable (soft prediction). The present study successfully attained an Area Under the Curve (AUC) value of 0.8507, which aligns with the benchmark set by (Mandrekar 2010).

Figure 10. Successes and correctness of the simulation.

3.8. Assessing land use and land cover dynamics: projections and trends under high urbanization scenario

This study presents a comprehensive assessment of LULC dynamics, utilizing a modelling approach, to project anticipated landscape configurations for the years 2032, 2042, and 2052 Figure 11. The projections are grounded in the High Urbanization Scenario (HUS), reflecting aggressive urbanization trends potentially driven by rapid population growth. This scenario serves as a basis for examining worst-case scenarios regarding forest loss. LULC dynamics are fundamental indicators of environmental change, with urbanization being a key driver of shifts in these patterns. Over the study period, the Coniferous Forest experienced a substantial loss of 79.82 square kilometres, signifying a 45.02% decrease (Table 6). Notably, the rate of decline in this forest class exhibited an acceleration between 2032–2042, followed by a deceleration between 2042–2052. This trend may be attributed to the implementation of conservation policies and heightened environmental awareness during the latter period. It is conceivable that governmental incentives for reforestation and stricter regulations on deforestation, such as the Billion Tree Tsunami Afforestation Program (BTFP) initiated in 2014, contributed to this observed minor increase. Subtropical Forest witnessed a significant decline of 285.65 square kilometres, equivalent to a 27.60% decrease from 2022 to 2052. The rate of change analysis underscored the consistent and severe annual loss of this forest type throughout the study period, indicating the challenges posed by persistent land use changes. Arable Land initially exhibited growth, expanding by 145.79 square kilometres (4.15%) from 2022 to 2032. However, it subsequently contracted by 46.47 square kilometres (1.26%) from 2032 to 2052. This fluctuation highlights the complexities in arable land management, influenced by evolving agricultural practices and regional factors. The Build Up category experienced a remarkable increase of 182.51 square kilometres (51.01%) from 2022 to 2032, followed by a further surge of 43.30% from 2032 to 2052. This data reveals a consistent annual expansion in urban areas throughout the study period, with a notable acceleration observed between 2022–2032. The findings emphasize the significant impact of urbanization on land cover dynamics, particularly the decline of forests, in the context of the High Urbanization Scenario (HUS). While Coniferous Forests and Subtropical Forests experienced substantial losses, the dynamics of Arable Land exhibited fluctuations, and Build-Up areas recorded consistent growth. These results underscore the importance of sustainable land use policies and conservation efforts in managing the environmental consequences of rapid urbanization.

Figure 11. Projected LULC map for 2032 (A), 2042(B), and 2052(C).

Table 6. Area (km² and %) for 2022, 2032, 2042 and 2052 and % of change of LULC for 2022–2032, 2032–2042 and 2042–2052.

	2022		2032		2042		2052		2022–2032		2032–2042		2042–2052
	Km^2	%	Km^2	%	Km^2	%	Km^2	%	Rate ∆	% of ∆	Rate ∆	% of ∆	Rate ∆	% of ∆
CF	177.3	3.4	99.7	1.9	88.7	1.7	97.5	12.9	−77.6	−43.7	−11.1	−11.1	8.8	9.9
STF	906.1	17.4	797.9	15.3	705.3	13.5	628.6	12.1	−108.2	−11.9	−92.5	−11.6	−76.7	−10.9
AL	3515.3	67.4	3661.1	70.2	3647.2	70.0	3614.7	69.3	145.8	4.1	−14.0	−0.4	−32.5	−0.9
BL	227.5	4.4	79.9	1.5	62.6	1.2	65.3	1.3	−147.5	−64.9	−17.4	−21.7	2.8	4.4
BU	357.8	6.9	540.3	10.4	676.5	13.0	774.2	14.8	182.5	51.0	136.2	25.2	97.7	14.4
WA	29.6	0.6	29.6	0.6	29.7	0.6	29.7	0.6	−0.1	−0.2	0.1	0.3	0.0	0.0
	5213.6		5208.6		5210.0		5210.0

*Rate ∆ represent rate of change between two intervals, while % of ∆ represent percent of change.

Over the period spanning from 2022 to 2052, sustained agricultural practices are projected to exhibit significant gains primarily originating from transitions in land cover categories, notably AL (969.32 km₂), STF (281.1 km²), and BL (24.11 km²), sequentially. The land cover categories associated with forests, denoted as CF and STF, are anticipated to undergo pronounced transformations primarily influenced by the expansion of AL (1.91 km² in CF and 362.36 km² in STF), respectively, alongside the encroachment of BL spanning 0.01 km² in CF and 7.71 km² IN STF, as well as BU replacing 0.11 km² in CF and 195 km² in STF, in the year 2032. Subsequently, these forest-based categories are expected to experience a decline in their extent during the ensuing years. The collective reduction of forest cover in 2032, diminishing from an initial area of 177.31 km² to 99.75 km² in the case of CF and from 906.10 km² to 797.88 km² for STF, can be primarily attributed to a confluence of factors, including urbanization, land conversion for alternative land uses, and shifts in agricultural practices. During the projected period, it is observed that Bare Land (BL) primarily replaced AL and BU land categories Table S8.

Analysis of land cover transitions spanning two decades from 2032 to 2052 reveals dynamic changes in the landscape. From 2032 to 2042, the transition matrix illustrates a shift from Coniferous Forest (CF) to Subtropical Forest (STF) (11.35 km²), signifying the expansion of subtropical forests and a subsequent transition to Agriculture Land (AL) (159.56 km²). By 2042, AL further transitions into Barren Land (BL) (1.17 km²) and Built-Up (BU) (118.91 km²) areas, reflecting evolving land use practices. Simultaneously, minimal changes are observed in Water Areas (WA), emphasizing the stability of aquatic ecosystems. The transition from 2042 to 2052 demonstrates a continued transformation of the landscape. STF expands, leading to an increase in AL, while BU areas gain prominence as shown in Figure S7. Notably, AL transitions to BL, indicating land degradation or altered land management. Water Areas remain relatively consistent, indicating the resilience of aquatic ecosystems. In the final transition from 2042 to 2052, CF undergoes a further reduction, and STF stabilizes, while AL continues to transition into BL (1.01 km²) and BU (85.82 km²) areas as shown in Table S8. BU experiences a notable gain (14.45%), reflecting urbanization and infrastructure development. These transitions reflect the evolving land use priorities and environmental dynamics, with increasing urbanization in the area. These combined transitions from 2032–2042 and 2042–2052 highlight the complex interplay of ecological and anthropogenic factors shaping land cover changes over time. The shift from forests to agriculture, urbanization, and potential land degradation underscore the evolving priorities and environmental dynamics within the study area.

4. Discussion

The results highlight the significant effects of human activities on the environment by demonstrating the changes in land use and land cover patterns that have occurred over the past few decades. The agricultural land has experienced a notable growth from 2939.07 km² in 1992 to 3515.3 km² in 2022. This expansion has been accompanied by a decrease in forested areas, specifically the subtropical forests which decreased from 1213.9 km² in 1992 to 1098.8 km² in 2002, and coniferous forests which decreased from 413.6 km² in 1992 to 336.4 km² in 2002. The decline of subtropical and coniferous forests is concerning due to habitat and biodiversity loss. Converting forests into urban areas disrupts ecosystem services like carbon sequestration and soil preservation (Cheng et al. 2023; Geng et al. 2023). Deforestation patterns have been noticed in several parts of the world, and numerous studies have emphasized the possible negative ecological consequences, such as the fragmentation and disappearance of connections across ecosystems (Adepoju and Salami 2017; Badshah et al. 2017; H. Wang and Qiu 2017).

Urbanization can lead to economic growth, but it can also cause environmental problems . As cities sprawl, pollution levels rise, natural landscapes are affected, and green spaces shrink, ultimately impacting the local climate and environment (Chen, Fang, and Liu 2023; Rahaman, Kalam, and Al-Mamun 2023; Shobairi et al. 2022; Usoltsev et al. 2020, 2022). Furthermore, the rapid urbanization witnessed in the research area corresponds to worldwide patterns, in which urban regions are growing at the cost of natural habitats, resulting in notable environmental challenges (Wilson et al. 2016; J. Xu et al. 2019). The phenomenon of deforestation to convert forests into agricultural fields is not exclusive to the specific location under investigation. Indeed, the loss of agricultural land has been recognized as a key factor contributing to deforestation in many locations (H. Wang and Qiu 2017). This trend is worsened by economic activities, as the demand for agricultural products and timber drives deforestation (Ferretti-Gallon and Busch 2014; Soe and Yeo-Chang 2019). The consequences of these alterations in land use are complex and have multiple aspects. While they satisfy the increasing need for food and resources (Aslam et al. 2022), they also contribute to environmental deterioration, reduction in biodiversity, and climate change (Andreevich et al. 2020; Mehmood, Anees, Rehman, et al. 2024). The transformation of 159.4 km² of coniferous forest and 432.5 km² of subtropical forest into other land use categories from 2002 to 2012 is not exclusive to the specific area being studied. Agricultural land losses have been highlighted as a significant cause of deforestation in several regions (Soe and Yeo-Chang 2019). The observed decrease in forested regions (Akram et al. 2022), including the subtropical and coniferous forests, has wider ramifications for climate change (Muhammad, Hamza, et al., 2023). Forests have a vital function in capturing and storing carbon dioxide, and their destruction can worsen the greenhouse effect (Birdsey et al. 2023; Moroni 2013). Moreover, deforestation and forest degradation can result in the irreversible depletion of forests, which has a significant impact on both the local and global climate (Anees, Yang, and Mehmood 2024; Hosonuma et al. 2012; Kyere-Boateng and Marek 2021; Pan et al. 2023). The findings of this study, which indicate a decrease in forest coverage and a rise in urban and agricultural regions, align with global observations where deforestation is influenced by both local needs and global demands (Kengoum et al. 2020; Lapola et al. 2023).

The study area has had a significant increase in urbanization, with built-up areas rising from 82.87 km² in 1992 to 344.10 km² in 2022, which is consistent with global patterns. Urbanization is encroaching upon natural environments, resulting in substantial environmental issues (Anees, Zhang, Khan, et al. 2022; Anees, Zhang, Shakeel, et al. 2022). The rate of growth of urban areas over a decade, which increased from 2.36% between 1992 and 2002 to 12.08% between 2012 and 2022, demonstrates the rapid speed at which cities are developing. The expansion of urban areas onto agricultural fields and woods signifies a change in land-use preferences, propelled by economic development and urbanization (Qiao, Guan, and Huang 2021). Urbanization can drive economic growth but requires responsible land use and environmental protection. Up-to-date regulations and sustainable development practices are vital for managing arable land amidst changing socio-economic trends (Abbas et al. 2020; Mahtta et al. 2022). Socioeconomic factors significantly affect land use and cover dynamics, including population growth, economic policies, and cultural practices. Integrating socioeconomic data into LULC models can aid in developing more sustainable and effective land management strategies (Filho et al. 2022; Hou et al. 2020). The conversion of arable lands into barren lands signifies the possibility of land degradation, which can have enduring consequences for both food security and the health of ecosystems (AbdelRahman, 2023). The study highlights the rising dependence on advanced machine learning techniques, specifically the Multilayer Perceptron Neural Network (MLP-NN), to analyse and predict LULC changes. The Multilayer Perceptron Neural Network (MLP-NN), consisting of input, hidden, and output layers, has proven its capacity to accurately represent complex patterns and behaviours in LULC datasets (Nguyen et al. 2020). The model’s ability to accurately simulate LULC patterns is demonstrated by achieving a skill score of 0.8724 and an accuracy rate of 89.08% (Alqadhi et al. 2021). Although the MLP-NN model demonstrated commendable performance, it is crucial to acknowledge its potential limits. The model’s dependence on physical parameters, without integrating socio-economic elements, may fail to encompass the complete range of factors influencing land use and land cover changes. It is important to mention that socioeconomic factors, such as population density, have a substantial impact on LULC changes (Girma, Fürst, and Moges 2022). Incorporating these socioeconomic variables into future algorithm iterations could improve the accuracy of predictions. The transition probability matrix and the transition area matrix, obtained through a Markov chain model, yielded valuable insights into the dynamics of transitions between various LULC categories (V. Kumar and Agrawal 2023b; Tariq and Mumtaz 2023). The matrices play a vital role in comprehending the probabilities linked to particular LULC transitions. They provide a measurable viewpoint on the potential changes in land use categories throughout time (Guarderas, Smith, and Dufrene 2022). In addition, the use of Multilayer Perceptron Neural Network (MLP-NN) and Cellular Automata Markov (CA Markov) models to predict Land Use and Land Cover (LULC) changes demonstrate the vast potential of machine learning and simulation (Joorabian Shooshtari and Aazami 2023). However, to enhance their predictive accuracy and application in policy-making, these models must evolve by incorporating more comprehensive variables such as climate change projections, economic trends, and policy shifts (Indraja, Aashi, and Krishna Vema 2024). The utilization of the CA Markov model for predicting future LULC patterns, specifically under the High Urbanization Scenario (HUS), provides insight into potential unfavourable situations. Scenario-based modelling plays a crucial role in helping policymakers and land managers anticipate and reduce potential negative consequences (Eastman 2020). Although the modelling methodologies used in this study have offered useful insights into land use and land cover changes, it is crucial to consistently improve and adjust these models. Expanding the scope of variables, including those related to physical and socio-economic factors, and conducting experiments with different scenarios can improve the predictive accuracy of the models and their usefulness in informing sustainable land management methods.

5. Conclusion

The thorough examination of the changes in LULC throughout 1992 to 2022 provides a significant understanding of the complex interaction between ecological and anthropogenic factors that influence the landscape of the region under investigation. Over the three-decade span, there has been a significant decrease in the extent of forested areas, notably in Coniferous and Subtropical Forests. This drop may largely be linked to the expansion of agricultural activities and the rapid growth of urban areas. The expansion of the Build-Up category, showing urban growth, shows the enormous impact of human activities on land cover dynamics, especially in the context of the High Urbanization Scenario (HUS). The MLP-NN and CA Markov models have been highly effective in capturing these complex patterns, resulting in noteworthy levels of accuracy. Nevertheless, the exclusion of socio-economic factors in the modelling approach implies that there are additional factors contributing to land use and land cover changes beyond just physical parameters. The future predictions, based on the Human Urban System (HUS), depict a scenario of ongoing urban growth that comes at the cost of forested regions. This highlights the importance of implementing sustainable land use policies and conservation initiatives. The study includes a comprehensive analysis of historical and current land use and land cover changes, as well as future estimates. Additionally, it emphasizes the significance of taking a holistic approach. By integrating socio-economic factors, improving modelling techniques, and addressing a wider range of future possibilities, we can improve our comprehension and control of land use and land cover dynamics. As urbanization and agricultural activities persist in transforming our landscapes, these studies become crucial in directing sustainable development and conservation efforts, guaranteeing a harmonious coexistence of human and natural ecosystems.

Author’s contribution

Khadim Hussain: conceptualization, methodology, software, formal analysis, visualization, data curation, writing – original draft, investigation, validation, writing – review and editing. Kaleem Mehmood: visualization, writing – review and editing, Sun Yujun: writing – review and editing, Supervision. Tariq Badshah: writing – review and editing. Shoaib Ahmad Anees: formal analysis, investigation, writing – review and editing. Fahad Shahzad: writing – review and editing. Nooruddin: writing – review and editing.Jamshid Ali: writing – review and editing. Bilal Muhammad: writing – review and editing. All authors have read and agreed to the published version of the manuscript.

Consent

Informed consent was obtained from all participants involved in this study. Participants were provided with detailed information about the research objectives, procedures, potential risks, and benefits before agreeing to participate. They were assured that their participation was voluntary, and they had the right to withdraw from the study at any time without facing any consequences.

All participants were informed about the confidentiality measures in place to protect their identity and personal information. Data collected during the study will be used solely for research purposes and will be securely stored.

This study was conducted in accordance with ethical standards and guidelines, and participants were encouraged to ask questions and seek clarification at any stage of the research process. If you have any further questions or concerns regarding the consent process, please contact [email protected].

Ethical approval

This research did not involve human or animal subjects; therefore, formal ethical approval was not required. The study strictly adheres to general ethical principles, and the authors are committed to upholding the highest standards of ethical research conduct. Any potential conflicts of interest that could have influenced the ethical conduct of this research have been declared.

Acknowledgements

We are grateful to the State Forestry and Grassland Administration Key Laboratory of Forest Resources and Environmental Management, Beijing Forestry University, Beijing (100083), P. R. China, for providing assistance and platforms for this research.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Correction Statement

This article has been corrected with minor changes. These changes do not impact the academic content of the article.

Additional information

Funding

This material is based on the original work conducted as part of my Ph.D. dissertation, generously supported by the State Forestry and Grassland Administration Key Laboratory of Forest Resources and Environmental Management and under the guidance of my supervisor, Yujun Sun. While no external funding sources were involved, the dedicated resources and mentorship provided by the lab and supervisor were instrumental in the completion of this research.