A Calculation and Optimization Method for the Theoretical Reclamation Timing of Cropland

: In mining areas with high groundwater tables, mining subsidence can lead to the inundation of cropland by water, causing damage to cropland and posing a threat to national food security. The implementation of concurrent mining and reclamation techniques can effectively enhance the reclamation rate of cropland. This technique requires engineers to initiate reclamation measures before cropland waterlogging occurs. Therefore, when mining a panel underground, an accurate calculation of the time when cropland becomes waterlogged, known as the theoretical reclamation timing, is crucial. To address this issue, this study proposes a computational method for the theoretical reclamation timing of cropland under the conditions of single-panel mining based on intelligent optimization algorithms. In addition, this paper also proposes an optimization method for the theoretical reclamation timing of cropland within a district based on an intelligent optimization algorithm. Utilizing this method makes optimizing the layout of multiple panels possible, thereby delaying the theoretical reclamation timing for cropland within a district. This approach aims to shorten the duration of reclamation projects and minimize their interference with agricultural activities. Through experimental validation, this paper demonstrates the reliability of these two methods. This study is beneficial for the rational planning of reclamation projects.


Introduction
China exhibits a considerable demand for coal [1,2].In 2022, coal consumption constituted 56.2% of China's total energy consumption [3].Given its massive population, China also presents a substantial need for grain.The significance of cropland in ensuring food security cannot be overstated.Significant spatial overlap is observed between the coal mining and grain production bases in eastern China.This overlapping area is commonly called the coal-grain overlapping area [4,5].However, underground mining can cause damage to the mining area environment, including the destruction of buildings and structures [6], reduction of agricultural output [7], soil contamination [8], surface ponding and pollution of surface waters [9,10], vegetation damage [11,12], and alteration of landscapes [13].In coal-grain overlapping areas, mining subsidence destroys cultivated soil's physical and chemical environments, ultimately reducing the grain yield [7,[14][15][16][17]. In addition, mining subsidence-induced cropland ponding has flooded significant tracts of farmland, leading to irreversible damage [18][19][20][21][22].
In China, reclamation measures are typically implemented for damaged cropland in mining areas.The technique of deep digging and shallow filling is a commonly used method for cropland reclamation [23].This technique involves digging the future deep subsidence area in advance to salvage more topsoil.The salvaged topsoil then fills the shallow subsidence area to restore more cropland.The digging areas will retain waterlogged areas, which can be transformed into vast fish ponds, thereby increasing the land's value for utilization [18,24].In addition, the creation of waterlogged areas can also positively affect the biodiversity of the landscape [25].The traditional reclamation approach applies reclamation techniques to cropland where subsidence has already ceased.At this point, the cropland has already been severely damaged, making it difficult to protect the cropland, resulting in a low reclamation ratio [1,2,26].To address the aforementioned concerns, Hu et al. proposed a concurrent mining and reclamation technology (CMR) [26].This technology advocates using mining subsidence control technology to excavate coal resources underground and implement reclamation measures on the surface before or during mining [4,26].Relevant studies have shown that this technology can improve the land reclamation rate [4,18,[27][28][29].This technique requires engineers to implement reclamation measures before cropland waterlogging occurs to salvage more soil.If the reclamation timing is set too early, the premature disturbance of cropland may ensue, adversely impacting farmers' cultivation practices despite the potential for topsoil preservation.On the contrary, if the reclamation timing is delayed, farmland may be submerged in water, reducing the reclamation rate.The reclamation timing of cropland is subject to influence by a myriad of factors, including but not limited to geological and mining conditions as well as human-driven factors.The ideal reclamation timing for cropland is when the cropland reaches the critical ponding state.At this time, more topsoil can be protected to restore more cropland.Furthermore, conducting reclamation during this period also helps to minimize the disruption to agricultural activities.This ideal reclamation timing is defined as the theoretical reclamation timing (TRT) [4].Thus, the precise computation of TRT for cropland is pivotal to the efficacious application of CMR technology in mining areas.
The hot research topics surrounding the "study of TRT for cropland" can be divided into two categories: (1) how to calculate the TRT for cropland when only a single longwall panel is mined [4,[30][31][32][33] and (2) how to optimize the layout of panels within the entire district to optimize the reclamation timing for cropland.For mining areas with high groundwater tables, cropland is already submerged in water after mining a single, wide, longwall panel.Therefore, reclamation of the cropland is required during the mining of a single panel.There are multiple panels in a district.Consequently, the reclamation period for the entire district is longer, which affects farmers' agricultural activities.Planning the layout of multiple panels within a district can effectively reduce the reclamation period for cropland.Hu et al. have proposed a "skip-mining" technique [23].This technique commonly involves reducing the width of multiple panels and altering their extraction sequence to control cropland subsidence.Controlling cropland subsidence can delay the timing of cropland waterlogging, thereby postponing the overall reclamation timing for cropland within the entire district [23,27].The detailed principles of the skip-mining technique will be described in Section 2.2.The objective of using the skip-mining technique is to delay the reclamation timing for cropland and shorten the overall reclamation period within a district.By reducing the duration of the reclamation period, the impact of reclamation projects on agricultural activities can be minimized.Optimizing the layout of panels within a district implies optimizing the TRT for cropland within this district as well.The optimization of the TRT for cropland within a district has emerged as a recent research focus [23].
However, calculating and optimizing the TRT is a deeply nonlinear problem with constraints.Currently, there is a lack of analytical expressions for solving the TRT.To address the two aforementioned issues, this study introduced, for the first time, the utilization of intelligent optimization algorithms to calculate and optimize the TRT for cropland.The method proposed in this study also applies to land reclamation projects in other coalproducing countries worldwide, such as Australia, the Czech Republic, and Poland.This study aims to establish a scientific foundation that could facilitate the effective implementation of CMR technology.

Calculation and Optimization Method for TRT
As shown in Figure 1, crop growth begins to be affected when the subsidence basin is at a certain height above the groundwater table, and the cropland is in a critical ponding state at this stage.

Calculation and Optimization Method for TRT
As shown in Figure 1, crop growth begins to be affected when the subsidence basin is at a certain height above the groundwater table, and the cropland is in a critical ponding state at this stage.The ideal reclamation timing for cropland is when the cropland reaches the critical ponding state.At this time, more topsoil can be protected to restore more cropland.This ideal reclamation timing is defined as the theoretical reclamation timing (TRT).In other words, the moment the cropland subsidence Wmax meets the following condition (Equation (1)) is deemed as the TRT for the cropland: where Hwater indicates the buried depth of the groundwater before cropland subsidence occurs; hb indicates the minimum soil thickness required for crop root growth.This paper aims to address two issues: (1) how to determine the TRT for cropland when only a single panel is mined underground and (2) how to optimize the layout of panels within the entire district to optimize the reclamation timing for cropland.Optimizing the reclamation timing for cropland is beneficial for reducing the reclamation period within the entire district and minimizing the impact of reclamation projects on agricultural activities.

Calculation Methodology for TRT under Single-Panel Mining Conditions
As shown in Figure 1, within the context of single-panel mining, the cropland reaches a critical ponding state when the panel advances to a certain distance.This distance (μ) is called the starting distance of critical cropland ponding (SD).The moment when a panel advances to the SD is the TRT of the cropland under the condition of single-panel mining.Indeed, determining the SD is equivalent to determining the TRT.
The maximum subsidence of cropland, denoted as Wmax(μ), can be calculated by using a specific mining subsidence prediction equation with the advancement distance of the panel μ as an independent variable.When Wmax(μ) = Hwater − hb, μ is the SD.Wmax(μ) is a dynamic mining subsidence prediction equation.Taking the probability integration method based on the Knothe time function as an example, this equation is a transcendental equation that involves the Gaussian error function [34].Due to the extremely nonlinear and complex nature of the Wmax(μ) equation, no analytical expression is available for accurately calculating the SD.The problem of calculating the SD can be considered a nonlinear engineering optimization problem with constraints.Intelligent optimization algorithms possess robust capabilities for inverting and resolving nonlinear model parameters.This section will propose a calculation method for the SD under the condition of single-panel mining based on intelligent optimization algorithms.The basic principle and calculation process of this method are as follows: The ideal reclamation timing for cropland is when the cropland reaches the critical ponding state.At this time, more topsoil can be protected to restore more cropland.This ideal reclamation timing is defined as the theoretical reclamation timing (TRT).In other words, the moment the cropland subsidence W max meets the following condition (Equation ( 1)) is deemed as the TRT for the cropland: where H water indicates the buried depth of the groundwater before cropland subsidence occurs; h b indicates the minimum soil thickness required for crop root growth.This paper aims to address two issues: (1) how to determine the TRT for cropland when only a single panel is mined underground and (2) how to optimize the layout of panels within the entire district to optimize the reclamation timing for cropland.Optimizing the reclamation timing for cropland is beneficial for reducing the reclamation period within the entire district and minimizing the impact of reclamation projects on agricultural activities.

Calculation Methodology for TRT under Single-Panel Mining Conditions
As shown in Figure 1, within the context of single-panel mining, the cropland reaches a critical ponding state when the panel advances to a certain distance.This distance (µ) is called the starting distance of critical cropland ponding (SD).The moment when a panel advances to the SD is the TRT of the cropland under the condition of single-panel mining.Indeed, determining the SD is equivalent to determining the TRT.
The maximum subsidence of cropland, denoted as W max (µ), can be calculated by using a specific mining subsidence prediction equation with the advancement distance of the panel µ as an independent variable.When W max (µ) = H water − h b , µ is the SD.W max (µ) is a dynamic mining subsidence prediction equation.Taking the probability integration method based on the Knothe time function as an example, this equation is a transcendental equation that involves the Gaussian error function [34].Due to the extremely nonlinear and complex nature of the W max (µ) equation, no analytical expression is available for accurately calculating the SD.The problem of calculating the SD can be considered a nonlinear engineering optimization problem with constraints.Intelligent optimization algorithms possess robust capabilities for inverting and resolving nonlinear model parameters.This section will propose a calculation method for the SD under the condition of single-panel mining based on intelligent optimization algorithms.The basic principle and calculation process of this method are as follows:

Fundamental Mathematical Model
The basic principle of using intelligent optimization algorithms to solve for the SD of a panel was as follows: Initially, the computer program generated several random initial solutions, denoted as µ 1 , µ 2 , µ 3 . ..µ i .The program checked if these solutions satisfy the constraint condition W max (µ i ) ≤ H water − h b .The solutions that did not meet the constraint condition were eliminated during the algorithm's iteration process.The solutions that satisfied the constraint condition underwent specific processing.The objective of the search was to find the maximum value of µ i within the range [L lower , L upper ] that satisfied the constraint condition.This maximum value represents the SD.This method typically requires multiple iterations to achieve the desired result.The process of solving the SD is essentially solving equations.
To sum up, the fundamental mathematical model of the SD calculation method, which is based on an intelligent optimization algorithm, can be expressed as follows: where µ i represents the advancing distance of the panel, [L lower , L upper ] denotes the search range of the SD, and W max (µ i ) ≤ H water − h b serves as a constraint in the calculation.From Equation (2), the fundamental principle of the SD calculation method based on intelligent optimization algorithms is to use specific intelligent optimization algorithms to search for the maximum advancing distance of the panel that satisfies the constraint conditions.There are various types of intelligent optimization algorithms [35], which can be classified into swarm intelligence optimization algorithms, evolutionary algorithms, and so on, based on their principles.Readers can select the relevant intelligent optimization algorithm to calculate the SD as per their requirements.This section calculated the SD using the classic particle swarm optimization (PSO) algorithm.The PSO algorithm is an evolutionary algorithm that simulates the foraging behavior of birds.The foraging behavior of birds shares similarities with the process of finding the optimal solution to a specific problem.Consequently, this algorithm is frequently employed to address optimization problems.The PSO algorithm has few control parameters, is easy to implement, and runs quickly.For a detailed explanation of the PSO algorithm, refer to Reference [36].

Computational Process of the SD Based on the PSO Algorithm
The steps for solving the SD based on the PSO algorithm were as follows: Step 1. Initialization of the algorithm's parameters: This step involves the initialization of the number of particles N, the maximum number of iterations G, inertia factor ω, particles' velocities v i , and positions µ i , learning factors c 1 and c 2 , and constraint factor r p .The aforementioned parameters are fixed parameters in the PSO algorithm.Before starting the iterative computation of the algorithm, it is necessary to initialize these parameters.
First, the search range for the SD was determined, denoted as [L lower , L upper ].Within this search range, N particles, denoted as µ 1 , µ 2 , µ 3 , µ i . ..µ n , were randomly generated using a computer program.According to the definition of the PSO algorithm, µ i represents the position of a particle and serves as the initial solution for SD.Secondly, each particle needed to be assigned an initial velocity.The velocity of each particle was randomly generated by the computer.Based on experience, the velocity range for each particle was [−(L upper − L lower ) × 0.1, (L upper − L lower ) × 0.1].The maximum absolute value of the particle's velocity, V max , was equal to (L upper − L lower ) × 0.1.
Step 2. Calculation of the fitness value: According to the definition of the PSO algorithm, the fitness value f (µ i ) of a particle is the criterion used to evaluate whether it represents the optimal solution.First, assessing whether each particle satisfies the constraint conditions is necessary.The algorithm's con-straint condition is W max (µ i ) ≤ H water − h b .W max (µ i ) represents the maximum subsidence value of cropland when the panel advances to a certain length µ i .The W max (µ i ) can be obtained by calculating it using a dynamic prediction equation for mining subsidence.There are many such equations, and readers can choose a dynamic prediction equation for mining subsidence that suits the research area.This section used the probability integration method based on the Knothe time function to calculate the W max (µ i ).This calculation method has been widely used in China.The basic principles of this method have been detailed in the literature [34].Due to space limitations, this paper did not provide a detailed explanation of the calculation method for W max (µ i ).
For any arbitrary particle µ i , its fitness value f was given by: From Equation (3), the fitness value of particles not satisfying the constraint condition was set to zero for easy elimination.This is because these particles cannot become the optimal solution for the SD.If a particle satisfied the constraint conditions, its fitness value was set as µ i (which is also the position of the particle).A larger value of µ i indicates that the particle is closer to the optimal solution for the SD.
Step 3. Comparison of the fitness value of each particle with its local best fitness value: The local best fitness value refers to the best fitness value achieved by this particle in history.If the fitness value of the particle was greater than its local best fitness value, its fitness value was set as this particle's local best fitness.The position of the particle corresponding to the local best fitness value can be represented as p l.
Step 4. Comparison of the fitness value of each particle with the global best fitness value, which represents the maximum fitness value attained by the entire particle swarm: If the fitness value of a particle was greater than the global best fitness value, the particle's current fitness value was set as the global best fitness value of the particle swarm.The position of the particle corresponding to the global optimal fitness value can be represented as p g .
Step 5.The velocity v i and position µ i of each particle should be updated using Equation (4): where r 1 and r 2 are two random numbers with a range of [0, 1].
Step 6.As shown in Equation ( 5), the updated particle velocity and position must be constrained within the ranges of [−V max , V max ] and [L lower , L upper ], respectively.The value of V max was usually set to 0.1 × (L upper − L lower ).
Step 7. If the maximum number, G, of iterations was reached, the algorithm was terminated.Then, the position of the particle corresponding to the global best fitness was outputted.This particle's position represents the optimal solution for the SD.Otherwise, Step 2 is repeated, and from there, the process continues.
By following the above steps, relevant computer programs can be developed to calculate the SD.The flowchart for the program design is shown in Figure 2. The moment when the advancing distance of the longwall panel reaches the SD is considered the TRT for cropland.
The aforementioned content described a computational method for the TRT of cropland under the conditions of single-panel mining based on the intelligent optimization algorithm.

Optimization of the TRT of Cropland within a Certain District
In mining areas with high groundwater levels, cropland ponding can be induced when excavating the first panel.Therefore, when a single panel is extracted, the cropland must be reclaimed.A district typically comprises multiple longwall panels.Therefore, the reclamation timing of cropland in the entire district is generally longer.However, a longer reclamation timing affects farmers' cultivation.To tackle this problem, some studies have suggested optimizing the dimensions of multiple panels in a district and formulating a continuation plan for panel excavation to postpone the reclamation timing by delaying surface ponding [4,23,27].This measure is conducive to shortening the cycle of land reclamation projects and reducing their impact on farmers' farming activities.Hu et al. adopted the skip-mining method to postpone the reclamation timing (see Figure 3) [4].This technique segregates the mining duration of a district into two phases, namely skip mining and back mining.The panels mined during the skip-mining and back-mining periods are, respectively, defined as the skip-panel and back-panel.Firstly, the skip-panel is excavated, followed by the back-panel.The sequential excavation was conducted on panel Nos. 1, 3, 5, and 7 during the skip-mining period, while panel Nos.2, 4, and 6 were excavated sequentially during the back-mining period (see Figure 3).During the skip-mining period, the surface subsidence was controlled within the critical subsidence value by optimizing the mining width.As a result, the cropland ponding was postponed to the backmining period, thereby postponing the reclamation timing, which is beneficial for mitigating the impact of land reclamation on farmers' agricultural operations.The purpose of using this mining technique is to optimize the TRT of cropland within a certain district.The moment when the advancing distance of the longwall panel reaches the SD is considered the TRT for cropland.
The aforementioned content described a computational method for the TRT of cropland under the conditions of single-panel mining based on the intelligent optimization algorithm.

Optimization of the TRT of Cropland within a Certain District
In mining areas with high groundwater levels, cropland ponding can be induced when excavating the first panel.Therefore, when a single panel is extracted, the cropland must be reclaimed.A district typically comprises multiple longwall panels.Therefore, the reclamation timing of cropland in the entire district is generally longer.However, a longer reclamation timing affects farmers' cultivation.To tackle this problem, some studies have suggested optimizing the dimensions of multiple panels in a district and formulating a continuation plan for panel excavation to postpone the reclamation timing by delaying surface ponding [4,23,27].This measure is conducive to shortening the cycle of land reclamation projects and reducing their impact on farmers' farming activities.Hu et al. adopted the skip-mining method to postpone the reclamation timing (see Figure 3) [4].This technique segregates the mining duration of a district into two phases, namely skip mining and back mining.The panels mined during the skip-mining and back-mining periods are, respectively, defined as the skip-panel and back-panel.Firstly, the skip-panel is excavated, followed by the back-panel.The sequential excavation was conducted on panel Nos. 1, 3, 5, and 7 during the skip-mining period, while panel Nos.2, 4, and 6 were excavated sequentially during the back-mining period (see Figure 3).During the skip-mining period, the surface subsidence was controlled within the critical subsidence value by optimizing the mining width.As a result, the cropland ponding was postponed to the back-mining period, thereby postponing the reclamation timing, which is beneficial for mitigating the impact of land reclamation on farmers' agricultural operations.The purpose of using this mining technique is to optimize the TRT of cropland within a certain district.The key to successfully implementing the skip-mining method is to design skip-panel and coal pillar widths between adjacent skip-panels (back-panel widths).If the width of the skip-panel is large or the coal pillar widths between adjacent skip-panels are small, significant surface subsidence values will still occur, leading to surface ponding.On the other hand, if the coal pillar width is overly large, pronounced undulating subsidence may occur in the cropland, thus adversely affecting agricultural operations (see Figure 4a).The key to successfully implementing the skip-mining method is to design skip-panel and coal pillar widths between adjacent skip-panels (back-panel widths).If the width of the skip-panel is large or the coal pillar widths between adjacent skip-panels are small, significant surface subsidence values will still occur, leading to surface ponding.On the other hand, if the coal pillar width is overly large, pronounced undulating subsidence may occur in the cropland, thus adversely affecting agricultural operations (see Figure 4a).The key to successfully implementing the skip-mining method is to design skip-panel and coal pillar widths between adjacent skip-panels (back-panel widths).If the width of the skip-panel is large or the coal pillar widths between adjacent skip-panels are small, significant surface subsidence values will still occur, leading to surface ponding.On the other hand, if the coal pillar width is overly large, pronounced undulating subsidence may occur in the cropland, thus adversely affecting agricultural operations (see Figure 4a).Considering that the design of the panel width is a nonlinear engineering optimization problem, this section proposed an intelligent optimization algorithm-based design method for skip-panel/back-panel widths.

Basic Principles
As shown in Figure 3, b represents the width of the skip-panel, while a denotes the spacing width between two adjacent skip-panels (the width of the back-panel is a).The optimal design values for a and b must possess the following characteristics: (1) Search range: The design values for a and b should not exceed the range restricted by the coal mining technology.For example, the maximum value for the panel width is typically not more than 400 m.This range can be shown as [D lower , D upper ].The design values for a and b must fall within this range.
(2) Constraint conditions: After mining according to the designed panel width, ensuring that the cropland is not in a critical ponding state is vital after completing the skip-mining phase.This condition imposes a constraint on the subsidence of the cropland.The constraint for cropland subsidence during the skip-mining phase can be expressed as W max (a,b) ≤ H water − h b (see Figure 3).
Another constraint is preventing the occurrence of undulating subsidence basins in cropland after the completion of excavation for each skip-panel, which is imperative to avoid disrupting agricultural activities.
As shown in Figure 4a, the undulating subsidence basin refers to the bottom of the surface subsidence basin exhibiting an obvious undulating waveform.The number of extreme points n j in an undulating subsidence basin is greater than one.Preventing the occurrence of undulating subsidence basins in cropland after the completion of excavation for each skip-panel is imperative to avoid disrupting agricultural activities.Thus, another constraint on this method can also be expressed as n j = 1 (as shown in Figure 4b).However, not all undulating subsidence basins will affect agricultural activities.Only when the height difference between the wave crest and trough exceeds a certain threshold can the surface subsidence basin be considered an undulating subsidence basin.In this study, the threshold was set to 10 mm.Suppose the surface subsidence basin still exhibits an undulating waveform during the skip-mining phase, but the height difference between the wave crest and trough is less than 10 mm.In this case, the number of extreme points n j in the surface subsidence basin is still considered to be one.
(3) Objective function: If the width of the skip-panel is overly small, the duration of the skip mining will be longer, which can negatively impact the efficiency of underground mining.The design aims to maximize the width of the skip-panel.The mathematical expression for this condition was max b.This condition represents the optimization objective for this design optimization.
The mathematical principle of the design method for the width of the skip-panel and back-panel based on an intelligent optimization algorithm can be expressed by Equation (6).
In Equation (6), max b indicates the objective function; g i represents the constraints.n j represents the number of extreme points in the surface subsidence basin.
From Equation ( 6), the general process for designing the width of the skip-panel and the back-panel based on intelligent optimization algorithms is as follows: Firstly, the computer randomly generates several initial solutions, represented as [a i , b i ].The first component of each solution represents the width of the back-panel, while the Then, each set of solutions is assessed to determine if it satisfies the two constraint conditions simultaneously.The fitness value of the solutions that meet the constraints is set to b i .The solutions that do not satisfy the constraints simultaneously are assigned a fitness value of 0 to eliminate them.Subsequently, a specific method is used to update the initial solutions.The updating methods vary depending on the different intelligent optimization algorithms.The solutions are then re-evaluated for compliance with the constraint conditions.At this stage, engineers are required to carry out several iterations of the calculation.The optimal design result is output when the termination condition for the iterative calculation is met.In summary, the basic principle of the design method for skip-panel/back-panel width based on intelligent optimization algorithms is to use specific intelligent optimization algorithms to design the widths of the skip-panel and back-panel to meet the given constraints.
There are many types of intelligent optimization algorithms, and readers can choose any specific algorithm to design panel widths [35].This section presented an instance of using the chicken swarm optimization (CSO) algorithm to design panel widths.The CSO algorithm was proposed by Meng in 2014 [37].This algorithm simulates the hierarchical mechanism and competitive behavior observed in a flock of chickens during foraging.The foraging behavior of the chickens is analogous to the process of finding optimal solutions to engineering problems.As a result, this algorithm is commonly employed to solve various engineering optimization problems.

Design Process of Skip-Panel and Back-Panel Widths Based on the CSO Algorithm
The basic flow of designing the panel widths for the skip-panel and back-panel based on the CSO algorithm can be outlined as follows: Step 1. Initializing the parameters of the algorithm: To begin, a group of chickens with a population size of "pop" is randomly generated.The variable "pop" represents the total number of chickens in a population.A flock consists of roosters, hens, and chicks.The number of roosters, hens, chicks, and mother hens with offspring in the chicken swarm are denoted as N R , N H , N C , and N M , respectively.The maximum number of iterations of the algorithm is set to GT.
The position of each chicken is represented as [a i , b i ].The first component of the position represents the width of the back-panel, while the second component represents the width of the skip-panel.Then, their respective fitness values are calculated based on the given criteria.The fitness of each chicken was computed as follows: First, it was necessary to determine whether each chicken satisfied the constraint condition.The first constraint condition of the algorithm is W max (a, b) ≤ H water − h b , where W max (a, b) refers to the maximum subsidence value in the reclamation area after completing the skip-mining phase.Usually, the subsidence value of each point in the reclamation area is calculated using the mining subsidence prediction equation, and the maximum subsidence value is then determined by aggregating these values.Various equations are available for mining subsidence prediction, and readers may choose the one that suits their study area.In this paper, the subsidence value of cropland was calculated using the probability integration method.The calculation formula for the probability integration method can be found in Reference [38].
The subsidence factor, q 0 , is a parameter in the probability integration method.The subsidence factor q 0 corresponding to different panel dimensions is not the same.We can assume that the width of a skip-panel is b, and its length is L. The calculation method for q 0 can be found in Equation (7) [30]: q 0 = 0.97n 2 − 0.07n + 0.39 q, 0.1 < n ≤ 0.83 q, n > 0.83 ( 7) where r represents the main influence radius.r = H/tanβ.H represents the mining depth of the panel m.The tangent of the main influence angle, tanβ, is a fixed parameter in the probability integration method.q represents the subsidence factor under the condition of critical mining.
The second constraint condition of the algorithm aims to avoid the occurrence of conspicuous undulating subsidence patterns at the bottom of the surface subsidence basin after the completion of excavation on each skip-panel (see Figure 4a).This is achieved by requiring that the number of extreme points, n j , in the subsidence basin equal 1 (see Figure 4b).It is generally sufficient to compute the number of extreme points, n j , on the subsidence curve of the major section of the subsidence basin to facilitate the programmatic determination of whether the subsidence basin exhibits undulating subsidence patterns.
The fitness value of chickens that do not satisfy the constraints should be set to 0, leading to their gradual elimination during evolution.In the previous section, we introduced the objective function of this design method as max b i .Under the premise of satisfying the constraints, a larger skip-panel width corresponds closer to the optimal design value.Therefore, if the position of each chicken satisfies the constraints, its fitness value is b i .
Step 2. Updating the relationships between the chickens in the flock every G generation: Based on the descending order of fitness values, each chicken was assigned as a rooster, hen, or chick.The chicken population was divided into several subgroups.Each subgroup consists of one rooster, several hens, and chicks.Once the hierarchy system of the chicken population is established, it remains unchanged until every Gth generation, when it is updated.
Step 3. Updating the position of the roosters using the following formula: We can suppose that the position of the i-th rooster at iteration tg is x i , t (tg); the updated formula for the rooster's position is as follows: (10) where randn (0, σ 2 ) represents a random number that follows a normal distribution, with an expectation of 0 and a variance of σ 2 .f i and f k represent the fitness values of the i-th rooster and k-th randomly selected rooster from the rooster group, respectively, where k ̸ = i.ε denotes a tiny value close to zero to avoid having a denominator of 0.
Step 4. Updating the position of the hens using the following formula: The updated formula for the hen's position is as follows: where rr 1 represents the rooster within the same sub-flock as the i-th hen.rr 2 indicates any other rooster or hen in the flock, excluding the i-th hen and the rooster rr 1 from the same sub-flock as the i-th hen.
Step 5. Updating the position of the chicks using the following formula: Step 6. Calculating the fitness values corresponding to the updated positions: If the updated fitness value was greater than the current best fitness value, then the position and fitness value of the individual needed to be replaced.Otherwise, there was no need to replace the position of the individual.
Step 7. Searching for the global maximum fitness value and its corresponding spatial position: If the current iteration number was greater than the maximum iteration number GT, the global best fitness value and its corresponding position were outputted.This position represents the optimal design value for the working width.Otherwise, Steps 2 to 7 are repeated.
Based on the aforementioned process, a relevant computer program can be developed to design the widths of the skip-panel and back-panel.The process of optimizing the design of panel width based on the CSO algorithm is shown in Figure 5.
where rr1 represents the rooster within the same sub-flock as the i-th hen.rr2 indicates any other rooster or hen in the flock, excluding the i-th hen and the rooster rr1 from the same sub-flock as the i-th hen.
Step 5. Updating the position of the chicks using the following formula: where mc indicates the hen corresponding to the i-th chick.FL is a random number in the interval [0, 2].After each update, each component of the new position must be restricted within the range [Dlower, Dupper].
Step 6. Calculating the fitness values corresponding to the updated positions: If the updated fitness value was greater than the current best fitness value, then the position and fitness value of the individual needed to be replaced.Otherwise, there was no need to replace the position of the individual.
Step 7. Searching for the global maximum fitness value and its corresponding spatial position: If the current iteration number was greater than the maximum iteration number GT, the global best fitness value and its corresponding position were outpu ed.This position represents the optimal design value for the working width.Otherwise, Steps 2 to 7 are repeated.
Based on the aforementioned process, a relevant computer program can be developed to design the widths of the skip-panel and back-panel.The process of optimizing the design of panel width based on the CSO algorithm is shown in Figure 5.  Theoretically, if the panel is excavated according to the width of the panel designed by this method, the cropland will not be ponding during the skip-mining period.As a result, the TRT for cropland can be postponed to the back-mining stage, which shortens the cycle of cropland reclamation projects and reduces their impact on farmers.
The aforementioned content described an optimization method for the TRT of cropland within a certain district.

Calculation Result of TRT for Cropland with Single Panel Simulated Mining in the Nantun Mining Area
Taking the geological and mining conditions of the Nantun mining area in Shandong Province, China, as an example, we simulated the extraction of a longwall panel (see Figure 6).The area above the panel is cropland.The geological and mining conditions of this panel are listed in Table 1.When the cropland reached a critical ponding state, the subsidence was 1300 mm (Hwater − hb = 1300 mm).Thus, the critical subsidence of the cropland was 1300 mm.Cropland subsidence reaches this critical value when this panel advances to a specific distance.At this point, reclamation is required for the cropland.Calculating this distance also implies calculating the TRT of the cropland.This section will use the method provided in Section 2.1 to calculate the TRT.A computation program was developed using Matlab (R2021a) software for this purpose.
In Step 2 of Section 2.1.2,the probability integration method based on the Knothe time function was used to calculate the dynamic maximum subsidence value Wmax(µi) of cropland.The values of the parameters in this equation, including the subsidence factor q0, tangent of main influence angle tanβ, deviation of inflection point S, influence  When the cropland reached a critical ponding state, the subsidence was 1300 mm (H water − h b = 1300 mm).Thus, the critical subsidence of the cropland was 1300 mm.Cropland subsidence reaches this critical value when this panel advances to a specific distance.At this point, reclamation is required for the cropland.Calculating this distance also implies calculating the TRT of the cropland.This section will use the method provided in Section 2.1 to calculate the TRT.A computation program was developed using Matlab (R2021a) software for this purpose.
In Step 2 of Section 2.1.2,the probability integration method based on the Knothe time function was used to calculate the dynamic maximum subsidence value W max (µ i ) of cropland.The values of the parameters in this equation, including the subsidence factor q 0 , tangent of main influence angle tanβ, deviation of inflection point S, influence transference angle θ 0 , and time influence coefficient c, are listed in Table 2.In addition, the relevant parameters for the PSO algorithm are listed in Table 2. Based on the computed results, the change curve of the fitness value is shown in Figure 7.As shown in Figure 7, the curve converged when the algorithm was iterated to the third generation, which means that when the number of iterations was three, the algorithm found the optimal calculation result.Based on the algorithm output, the SD was 422 m; at this point, the maximum cropland subsidence was 1300 mm, which had not yet exceeded its critical subsidence threshold.To validate the accuracy of the computed results, this study calculated the maximum dynamic subsidence of the cropland when the panel advanced by 423 m.The obtained value was 1303 mm, surpassing the critical subsidence threshold.This indicated that 422 m was the optimum solution for SD.This indicated that the method proposed in this paper accurately calculated the SD.When the panel advanced to a distance of 422 m, it was the theoretical reclamation timing for cropland.

Optimization of the TRT of Cropland within a District in Nantun Mining Area
The geological and mining conditions of the Nantun mining area were used as a basis for simulating the extraction of a district, as shown in Figure 8.The district's dimensions were 1250 m in length and 900 m in width.The geological conditions of this district are the same as those stated in the case study of Section 3.1 (refer to Table 3).As shown in Figure 7, the curve converged when the algorithm was iterated to the third generation, which means that when the number of iterations was three, the algorithm found the optimal calculation result.Based on the algorithm output, the SD was 422 m; at this point, the maximum cropland subsidence was 1300 mm, which had not yet exceeded its critical subsidence threshold.To validate the accuracy of the computed results, this study calculated the maximum dynamic subsidence of the cropland when the panel advanced by 423 m.The obtained value was 1303 mm, surpassing the critical subsidence threshold.This indicated that 422 m was the optimum solution for SD.This indicated that the method proposed in this paper accurately calculated the SD.When the panel advanced to a distance of 422 m, it was the theoretical reclamation timing for cropland.

Optimization of the TRT of Cropland within a District in Nantun Mining Area
The geological and mining conditions of the Nantun mining area were used as a basis for simulating the extraction of a district, as shown in Figure 8.The district's dimensions were 1250 m in length and 900 m in width.The geological conditions of this district are the same as those stated in the case study of Section 3.1 (refer to Table 3).
The critical subsidence value for cropland was calculated to be 1300 mm using Equation (1).Multiple panels could be arranged within this district.To postpone the TRT, the skip-mining method was adopted to extract the coal resources in this district.As shown in Figure 8, several skip-panels were planned to be arranged in this district (No. 1 skip-panel; No. 2 skip-panel. ..).The advancement length of each panel was set to a fixed value of 900 m.The optimization of the width b of the skip-panel and the width a of the back-panel was necessary to avoid ponding on cropland during the skip-mining period.This section designed the width of the skip-panel and back-panel using the method presented in Section 2.2.A computation program was developed using Matlab (R2021a) software for this purpose.
In Step 1 of Section 2.2.2, the maximum subsidence of cropland, W max (a, b), was calculated using the probability integration method.The subsidence factor q under critical mining conditions in this mining area was 0.9.Due to the relatively small width of the skip-panel, the subsidence factor q 0 for each skip-panel should be adjusted according to its dimensions.The subsidence factor q 0 of each skip-panel should be corrected according to Equations ( 7)- (9).The values of other parameters (tanβ, θ 0 , and S) for the probability integration method are listed in Table 4.The relevant parameters of the CSO algorithm are listed in Table 4. culated using the probability integration method.The subsidence factor q under critical mining conditions in this mining area was 0.9.Due to the relatively small width of the skip-panel, the subsidence factor q0 for each skip-panel should be adjusted according to its dimensions.The subsidence factor q0 of each skip-panel should be corrected according to Equations ( 7)-( 9).The values of other parameters (tanβ, θ0, and S) for the probability integration method are listed in Table 4.The relevant parameters of the CSO algorithm are listed in Table 4.     Based on the computed results, the change curve of fitness value is shown in Figure 9.As shown in Figure 9, when the algorithm was iterated to the 174th generation, the algorithm found the optimal design result.The algorithm results showed that the optimal width of each skip-panel was 130 m, with a distance of 69 m between adjacent panels (the optimal width of back-panel).As shown in Figure 10, based on the calculation results, a total of six skip-panels were arranged within the entire district.Figure 11 shows the sub- As shown in Figure 9, when the algorithm was iterated to the 174th generation, the algorithm found the optimal design result.The algorithm results showed that the optimal width of each skip-panel was 130 m, with a distance of 69 m between adjacent panels (the optimal width of back-panel).As shown in Figure 10, based on the calculation results, a total of six skip-panels were arranged within the entire district.Figure 11 shows the subsidence basin after the extraction of each skip-panel.The subsidence curve on the major section of the subsidence basin along the dip direction is shown in Figure 12.As shown in Figure 9, when the algorithm was iterated to the 174th generation, the algorithm found the optimal design result.The algorithm results showed that the optimal width of each skip-panel was 130 m, with a distance of 69 m between adjacent panels (the optimal width of back-panel).As shown in Figure 10, based on the calculation results, a total of six skip-panels were arranged within the entire district.Figure 11 shows the subsidence basin after the extraction of each skip-panel.The subsidence curve on the major section of the subsidence basin along the dip direction is shown in Figure 12.As shown in Figures 11 and 12, after each skip-mining panel's extraction was completed, no noticeable undulating surface subsidence was observed within the subsidence basin.After the extraction of all six skip-mining panels, the maximum subsidence of cropland was 1299 mm, which remained below the critical subsidence threshold.Therefore, mining was carried out according to the width of the skip-panel designed by this method, and the TRT of cultivated land was successfully postponed to the back-mining stage, which is beneficial for mitigating the impact of land reclamation on farmers' agricultural operations.As shown in Figures 11 and 12, after each skip-mining panel's extraction was completed, no noticeable undulating surface subsidence was observed within the subsidence basin.After the extraction of all six skip-mining panels, the maximum subsidence of cropland was 1299 mm, which remained below the critical subsidence threshold.Therefore, mining was carried out according to the width of the skip-panel designed by this method, and the TRT of cultivated land was successfully postponed to the back-mining stage, which is beneficial for mitigating the impact of land reclamation on farmers' agricultural operations.

Influence of Panel Size on the Subsidence of Cropland
According to the general rules of mining subsidence, the subsidence of the cropland depends on the mining degree along the direction of the panel's advancement and the mining degree along the direction of the panel width [39].The advancement distance of the panel is generally quite long, which easily results in reaching a state of critical mining.However, along the direction of the panel width, the cropland is generally in subcritical mining.Therefore, in the context of longwall mining, the subsidence of the cropland primarily depends on the panel width.When the panel width is small, some key strata in the overburden do not reach their limit spans and will not fracture, thus controlling surface subsidence.Engineers have utilized this principle to develop a strip-mining technology to control surface subsidence [39].The subsidence control principle of the skip-mining method is similar to that of the strip-mining technology, which involves controlling surface subsidence by reducing the panel width.
The cases in Sections 3.1 and 3.2 are examples of simulated mining, with identical geological conditions for both cases.In the case study of Section 3.1, the panel width is 210 m.When this panel has advanced 422 m, the surface subsidence has already reached the critical subsidence value for cropland (1300 mm).In the case study of Section 3.2, the width of the skip-panel is only 130 m, and even when the panel has advanced 900 m, the cropland subsidence does not exceed its critical subsidence value.This is because the width of the skip-panel is smaller, resulting in a lower mining degree along the direction of its width.The coal pillars on both sides of the skip-panel can provide stable support for the overburden, preventing the large-scale destruction of some key strata, hence keeping the subsidence of the cropland relatively low.Additionally, the phenomenon of smaller cropland subsidence during the skip-mining stage can also be explained from the perspective of subsidence prediction.When employing the probability integration method to predict cropland subsidence, the subsidence factors corresponding to panels with different widths are different under a subcritical mining state.In the case study of Section 3.1, with a panel width of 210 m, the resulting subsidence factor q 0 is approximately 0.65, determined by applying the calculation methods from Equations ( 7)- (9).In comparison, in the case study of Section 3.2, with a panel width of 130 m, the subsidence factor q 0 is about 0.53.The lower the value of the subsidence factor, the less the predicted surface subsidence.This further explains why reducing the panel width to 130 m does not cause the subsidence of cropland to exceed its critical value, even when multiple skip-panels are arranged.Therefore, optimizing the width of a skip-panel for controlling subsidence and delaying the TRT of cropland is of great importance.
The width of coal pillars (the back-panel width) also affects the subsidence of cropland.This is because the cropland subsidence caused by the extraction of two skip-panels on either side of a coal pillar will superimpose.In the case study of Section 3.2, when mining the first skip-panel, the maximum subsidence value of the cropland is 907 mm; when mining the second skip-panel, the maximum surface subsidence increases to 1234 mm (see Figure 12).The increase in the maximum surface subsidence is due to the cumulative effect of subsidence caused by the extraction of the two skip-panels.When mining the fifth and sixth skip-panels, the maximum subsidence value of the surface remains at 1299 mm, no longer significantly increasing and not exceeding the critical subsidence value for cropland.Therefore, if the width of the coal pillars is small, the cumulative cropland subsidence value can easily exceed the critical subsidence threshold of the cropland.Conversely, when the width of the coal pillars is large, it can easily lead to the formation of undulating subsidence basins in farmland, affecting the farmers' cultivation.Therefore, optimizing the width of the coal pillars between skip-panels is also very important.

Outlook and Limitations of the Study
Calculating and optimizing the TRT of cropland is a relatively complex engineering optimization problem that lacks reliable analytical expressions to solve for the TRT.To address this issue, this study introduced, for the first time, the use of intelligent optimization algorithms to solve the planning and design problems in concurrent mining and reclamation projects.The experiments in Section 3 also demonstrate the reliability of this method.Certainly, this study encourages readers to explore more efficient intelligent optimization algorithms for calculating and optimizing the TRT.The utilization of intelligent optimization algorithms is not limited to the PSO and CSO algorithms employed in this study.Additionally, this study encourages readers to select prediction methods with a higher accuracy for predicting cropland subsidence, and there is no need to be confined to the probability integration method used in this study.In summary, this study merely provides engineers with a fresh perspective for calculating and optimizing the TRT.
The TRT studied in this paper is an idealized timing for cropland reclamation.Nevertheless, the influence of variables such as precipitation and anthropogenic factors on the reclamation timing has not been accounted for, which is a limitation of our research.

Conclusions
This study introduced, for the first time, the application of intelligent optimization algorithms in the context of concurrent mining and reclamation engineering.A computational method based on the PSO algorithm was developed to determine the TRT for cropland under the condition of single panel extraction.The reliability of this method was demonstrated through simulation experiments.Furthermore, this study employed the CSO

Figure 2 .
Figure 2. Calculation method for SD based on PSO algorithm.

Figure 2 .
Figure 2. Calculation method for SD based on PSO algorithm.
the width of the skip-panel.Both components must take values within the range [D lower , D upper ].

Land 2024, 13 , 638 11 of 19 where
mc indicates the hen corresponding to the i-th chick.FL is a random number in the interval [0, 2].After each update, each component of the new position must be restricted within the range [D lower , D upper ].

Figure 5 .
Figure 5. Optimization process for widths of the skip-panel and back-panel based on the CSO algorithm.

Figure 5 .
Figure 5. Optimization process for widths of the skip-panel and back-panel based on the CSO algorithm.

Table 1 .Figure 6 .
Figure 6.The position and size of the panel.

Figure 6 .
Figure 6.The position and size of the panel.

Table 1 .
The geological and mining conditions.
Based on the computed results, the change curve of the fitness value is shown in Figure 7.

Figure 7 .
Figure 7.The change curve of the fitness value based on the PSO algorithm.

Figure 7 .
Figure 7.The change curve of the fitness value based on the PSO algorithm.

Figure 8 .
Figure 8. Schematic diagram of a layout of skip-panel and back-panel within the district.
computed results, the change curve of fitness value is shown in Figure9.

Figure 8 .
Figure 8. Schematic diagram of a layout of skip-panel and back-panel within the district.

Land 2024 , 20 Figure 9 .
Figure 9.The change curve of the fitness value based on the CSO algorithm.

Figure 9 .
Figure 9.The change curve of the fitness value based on the CSO algorithm.

Figure 9 .
Figure 9.The change curve of the fitness value based on the CSO algorithm.

Figure 10 .
Figure 10.Optimized layout of the skip-panel and back-panel.Figure 10.Optimized layout of the skip-panel and back-panel.

Figure 10 .
Figure 10.Optimized layout of the skip-panel and back-panel.Figure 10.Optimized layout of the skip-panel and back-panel.Land 2024, 13, x FOR PEER REVIEW 16 of 20

Figure 12 .
Figure 12.Subsidence curve on the major section of the subsidence basin along the dip direction.

Figure 12 .
Figure 12.Subsidence curve on the major section of the subsidence basin along the dip direction.

Table 2 .
Parameters for subsidence prediction equation and PSO algorithm.

Table 3 .
The geological conditions of the district.

Table 4 .
Parameters for probability integration method and CSO algorithm.

Table 3 .
The geological conditions of the district.

Table 4 .
for probability integration method and CSO algorithm.