A multi-objective coupled geographical stratification method based on spatial continuity constraints

Abstract

The application discloses a multi-target coupling geographical stratification method and system based on spatial continuity constraint, which firstly carries out quality screening on geographical index and explanatory factor data; subsequently, an initial population is generated by adopting an NSGA-II algorithm, and each individual represents a stratification scheme defined by a hyperplane parameter; for each scheme, a geographical probe q statistic (measuring homogeneity in the layer), an adjacent homogeneity index and a blockiness index are synchronously calculated, and are collectively used as a target function of multi-target optimization; population evolution is carried out through non-dominated sorting, congestion calculation, simulated binary crossover and polynomial mutation, and a simulated annealing mechanism is introduced to enhance the global search capability; finally, spatial smoothing is carried out on the optimization result through an adaptive window model, and geographical stratification results meeting statistical significance and maintaining spatial continuity are output, so that the spatial fragmentation and fragmentation problems in the traditional method are effectively solved.

Description

Technical Field

[0001] This invention belongs to the field of geographic information science and technology, and in particular relates to a multi-objective coupled geographic layering method based on spatial continuity constraints. Background Technology

[0002] Geographic stratification, as a crucial tool for analyzing the spatial differentiation of geographical phenomena and uncovering their spatial characteristics and driving mechanisms, has wide-ranging applications in regional planning, environmental management, resource assessment, and disaster risk analysis. With the continuous development of Geographic Information Systems (GIS) and spatial data analysis methods, geographic stratification has gradually shifted from traditional qualitative description to a quantitative and model-based research paradigm. The industry is generally committed to revealing the distribution patterns and formation mechanisms of geographical elements through spatial stratification to support comprehensive geographical decision-making across multiple scales and factors. Currently, geographic stratification methods are evolving towards multi-indicator, multi-process collaborative modeling to adapt to the needs of complex geographical system simulation and multi-dimensional spatial optimization under sustainable development goals.

[0003] In existing technologies, geographic stratification methods mainly fall into two categories: single-objective and multi-objective. Single-objective models focus on a single geographic indicator, modeling its spatial stratification heterogeneity through one or more explanatory variables. Multi-objective geographic stratification methods, on the other hand, consider multiple geographic indicators simultaneously, aiming to reveal the spatial heterogeneity structure exhibited by these indicators through common or related explanatory factors. Early multi-objective methods typically employed weighted synthesis to transform multiple objectives into a single comprehensive objective, or optimized the order of objectives by manually setting their priorities, thereby achieving collaborative modeling of multi-dimensional indicators.

[0004] Although multi-objective evolutionary algorithms have achieved multi-objective collaborative hierarchical structure to some extent and effectively overcome the subjectivity problem of weight setting, existing methods still have a most obvious drawback: they do not adequately model the structural and relational characteristics inherent in geospatial space itself, and fail to explicitly incorporate key attributes such as spatial autocorrelation and spatial dependence into the optimization objective function. This often leads to problems such as fragmented spatial units and inconsistencies between boundary delineation and the actual geographical pattern in the hierarchical results, limiting their ability to deeply reveal the spatial coupling mechanisms and driving factors of geographical phenomena. Summary of the Invention

[0005] Purpose of the invention: The purpose of this invention is to provide a multi-objective coupled geographic stratification method based on spatial continuity constraints that incorporates spatial autocorrelation considerations, can reasonably stratify complex geographic phenomena, and improves the practicality and geographic interpretability of the results.

[0006] Technical solution: The multi-target coupled geographic layering method of the present invention includes the following steps:

[0007] (1) For the target study area and the geographical phenomena to be analyzed, select at least two indicator variables and one or more explanatory factors to describe the characteristics of geographical response. By constructing a dynamic mask, the indicator variables and explanatory factors are removed to obtain a geographical dataset without null values.

[0008] (2) Randomly generate an initial population containing N individuals. Each individual is encoded with a set of hyperplane parameters. Multiple hyperplanes are constructed based on each set of hyperplane parameters. According to the spatial segmentation logic corresponding to the hyperplane, each geographic unit in the target study area is divided into the corresponding stratum. Each individual corresponds to an independent candidate stratification scheme.

[0009] (3) For each candidate stratification scheme, calculate the geospatial detector q value of each indicator variable under the stratification scheme, combine the adjacency homogeneity index and block morphology index to construct a multi-objective optimization function, perform non-dominated sorting and crowding calculation on individuals in the current population, and select parent individuals based on binary tournament to form the parent population.

[0010] (4) For the selected parent individuals, perform dynamic simulated binary crossover. Based on the difference in objective function between the parent individuals, adaptively adjust the crossover range of the simulated binary crossover to generate preliminary offspring individuals. Perform polynomial perturbation on the hyperplane parameter vector corresponding to the preliminary offspring individuals to generate offspring individuals and form the offspring population.

[0011] (5) Merge the parent population and the offspring population to form a new composite population; select a new generation of parent individuals from the composite population based on a binary tournament; introduce simulated annealing local search, select the top K elite individuals in the non-dominated sorting to perform corresponding neighborhood perturbations, and adopt a global optimal solution retention mechanism.

[0012] (6) Repeat steps (3) to (5). When the number of iterations reaches the preset maximum number of iterations, terminate the evolution and output the optimal region stratification result from the non-dominated sorting result of the final generation.

[0013] (7) Adaptive window model is used to perform spatial smoothing on the regional stratification results of the iterative output to eliminate discrete geographic units.

[0014] This invention ensures the integrity and reliability of the geographic dataset by constructing a dynamic mask to remove null values, laying a solid foundation for subsequent analysis. It provides a broad search space for the optimization algorithm by randomly generating diverse candidate stratification schemes through hyperplane parameter encoding. The core of this invention lies in constructing a multi-objective function that integrates the q-value of the geographic detector, the adjacency homogeneity index, and the block morphology index. This directly incorporates spatial autocorrelation considerations and simultaneously optimizes the statistical heterogeneity of stratification and the rationality of spatial morphology, thereby achieving the scientific quantification and stratification of complex geographic phenomena. Furthermore, it employs adaptive crossover and polynomial perturbation for evolutionary operations, effectively exploring new solutions while inheriting the excellent characteristics of parent solutions, avoiding premature convergence. Simulated annealing local search and elite retention mechanisms further explore the potential of high-quality solutions and ensure the generational inheritance of optimal solutions, significantly improving the quality of solutions and the convergence efficiency of the algorithm. Through iterative evolution, a Pareto optimal solution set is selected, ensuring the comprehensive optimality of the final stratification scheme. Finally, adaptive window smoothing is used to remove small discrete units, enhancing the spatial continuity and visual aesthetics of the partitioned block structure, greatly improving the practicality and geographic interpretability of the stratification results.

[0015] Preferably, the dynamic mask described in step (1) is defined as:

[0016]

[0017] In the formula, i and j represent the row and column numbers, respectively, and the number of rows and columns is defined as... Construct a mask array .

[0018] By constructing dynamic masks, the automatic quality screening of source data for geographical indicators and explanatory factors was effectively realized, eliminating the influence of outliers, providing a high-quality data foundation for subsequent multi-objective hierarchical optimization, and enhancing the robustness of the model.

[0019] Preferably, the hyperplane parameters of each individual in step (2) are defined as follows:

[0020]

[0021] In the formula, j is the j-th hyperplane. Let b be the hyperplane parameter corresponding to the i-th explanatory factor, and b be a constant term.

[0022] For the i-th pixel, its discriminant function on the j-th hyperplane is:

[0023]

[0024] In the formula, This represents the value of the nth explanatory factor in the i-th geographic unit;

[0025] The classification rules are as follows:

[0026]

[0027] In the formula, To determine the category of the i-th pixel in the j-th hyperplane, calculate the category of the i-th pixel in all hyperplanes, and finally obtain... The category label, h represents the number of hyperplanes.

[0028] By defining hyperplane parameters and discriminant functions, the complex problem of geographic stratification is transformed into a linear partitioning problem in high-dimensional space. This allows each geographic unit to be accurately classified into different strata based on the comprehensive characteristics of multiple explanatory factors, thereby achieving effective capture and structured expression of multi-dimensional and nonlinear relationships of geographic phenomena and enhancing the mathematical rigor and geographic interpretability of the stratification model.

[0029] Preferably, the formula for calculating the q-statistic value of the geographic detector in step (3) is as follows:

[0030]

[0031] In the formula, Let q represent the geodetector q value obtained by classifying the data based on the hyperplane parameters of the i-th individual, with a value range of [0,1]. The closer a value is to 1, the stronger the explanatory power of the classification for geographical phenomena within the region.

[0032] Calculate the spatial continuity index under each stratification scheme:

[0033] Adjacent cells in a raster are compared. If two adjacent cells belong to the same category, they are considered homogeneous; otherwise, they are considered heterogeneous. The Adjacency Homogeneity Index (AH) is defined as follows:

[0034]

[0035] The range of AH is [0,1]. The closer AH is to 1, the more adjacent pixels belong to the same category, that is, the stronger the spatial continuity. The closer AH is to 0, the greater the difference between adjacent pixels, that is, the more fragmented the space.

[0036] Based on the boundary complexity, area size, and shape compactness of the blocks, the ability of similar pixels to form larger, regular, continuous blocks is determined using the Block Morphology Index (BI), calculated as follows:

[0037]

[0038] In the formula, P represents the perimeter of the block, and A represents the area of ​​the block.

[0039] By combining the q-statistic of the geographic detector to measure the explanatory power of the stratification scheme for geographic indicators, and introducing the adjacency homogeneity index and block morphology index to comprehensively assess spatial continuity, this step effectively balances statistical heterogeneity and the integrity of spatial patterns. This results in a final stratification that not only fully reveals the inherent laws of geographic phenomena, but also has good spatial aggregation and boundary rationality, significantly improving the applicability and practical geographic significance of the results.

[0040] Preferably, the objective function described in step (3) is designed as follows:

[0041]

[0042] In the formula, This represents the clustering target value of the nth y-factor. Represents the adjacency homogeneity index. Indicates the mass morphology index. This represents the value of a multi-objective optimization function, used to comprehensively reflect the overall performance of candidate stratification schemes in terms of both geographical explanatory power and spatial continuity.

[0043] This objective function integrates the geodetector q-value, adjacency homogeneity index, and block morphology index to couple the traditional stratification objective, which focuses on statistical differences, with the tightness of spatial adjacency relationships and the regularity of block morphology. This mathematically unifies and synergistically optimizes the geographical explanatory power and spatial continuity of the stratification scheme, ultimately ensuring that the resulting stratification results are not only statistically significant but also have high continuity and regularity in spatial pattern.

[0044] Preferably, the non-dominated sorting in step (3) includes:

[0045] Compare any two individuals in the initial population based on n geographical indicators:

[0046]

[0047] In the formula, Let represent the target vector of individual a on the i-th geographical indicator. Let this be the target vector of individual b on the i-th geographic indicator; if for any geographic indicator, All greater than or equal to a is said to dominate b; if for some target i, there exists Meanwhile, other objectives also exist. If the dominance relationship is not established, it indicates that the two individuals have their own advantages in different objectives and cannot establish a dominance relationship. Therefore, the relationship between individual a and individual b is determined to be non-dominant. By performing the above dominance determination process on any two individuals in the population, the dominance matrix of the entire population is obtained. The system stratifies all individuals according to the dominance matrix to form different non-dominant frontier hierarchical structures.

[0048] The congestion calculation includes:

[0049] In the process of selecting a preset number of dominant individuals by level, if the total number of individuals in a certain level exceeds the planned quota of dominant individuals for that level, then the crowding degree of all individuals within that level is quantitatively calculated:

[0050]

[0051] In the formula, Let m be the crowding degree of the j-th individual, and m be the initial total number of individuals in the population. Let q represent the geodetector value of the x-th individual on the i-th geographic indicator.

[0052] By using non-dominated sorting to finely stratify the population based on the dominance relationship of multi-objective functions, we can ensure that elite individuals with comprehensive advantages in various geographical indicators are selected first. In addition, by combining crowding calculation, we can prioritize solutions with sparser spatial distribution of objectives among individuals with similar performance. This effectively maintains the diversity and broad representativeness of the solution set, avoids the optimization process from getting trapped in local optima, and ultimately leads to a set of high-quality stratified solutions that achieve Pareto optimality and are evenly distributed in terms of explanatory power, spatial continuity and morphological regularity.

[0053] Preferably, after the non-dominated ranking and crowding calculation described in step (3), the process further includes screening offspring individuals using a binary tournament:

[0054]

[0055] In the formula, a.rank represents the non-dominant ranking rank of individual a, and b.rank represents the non-dominant ranking rank of individual b. Let 'a' represent the crowding level of individual 'a'. Rank individuals b by their crowding level; The comparison function is as follows: if individual a has a better non-dominated ranking level, then individual a is directly retained; if individual a and individual b belong to the same non-dominated ranking level, then crowding is used as a secondary criterion, and the individual with the smaller crowding value is selected to enter the subsequent population. Through pairwise competition logic, offspring individuals with both better ranking and lower crowding are gradually selected.

[0056] By introducing a binary tournament selection mechanism and comprehensively comparing and selecting individuals based on non-dominated sorting levels and crowding distance, this step ensures that the evolutionary process always prioritizes retaining elite solutions that are in better non-dominated levels and can increase population diversity. This effectively balances the convergence speed of the algorithm with the distribution breadth of the solution set, avoids premature convergence, and ultimately leads to obtaining a batch of high-quality Pareto optimal solutions that are excellent in performance and evenly distributed across multiple optimization objectives.

[0057] Preferably, the crossover mutation calculation steps in step (4) are as follows:

[0058] By aggregating multiple objective functions, a parent individual evaluation index is constructed, and the geographical explanatory power index Q is:

[0059]

[0060] In the formula, This represents the clustering target value of the first y-factor. This represents the clustering target value of the second y-factor. This represents the clustering target value for the third y-factor;

[0061] The spatial continuity index S is aggregated as follows:

[0062]

[0063] In the formula, Represents the adjacency homogeneity index. The weights for the adjacency homogeneity index; Indicates the mass morphology index. The weights for the adjacency homogeneity index;

[0064] The comprehensive evaluation index W is:

[0065]

[0066] In the formula, Q represents the geographical explanatory power index. represents the weight of the geographical explanatory power index; S represents the spatial continuity index. The weights for the spatial continuity index;

[0067] The multi-objective differences among parent individuals are:

[0068]

[0069] In the formula, This represents the comprehensive evaluation index of the stratification scheme corresponding to parent individual i. This represents the comprehensive evaluation index of the stratification scheme corresponding to the parent individual j. The D value represents the difference in overall fitness between two individuals; the larger the D value, the greater the difference in their overall performance across multiple objectives; the smaller the D value, the more similar their performance, which can be used for local refinement.

[0070] Based on the degree of difference Adaptive adjustment of cross-range control parameters The greater the difference, The larger the value, the wider the range of intersection; the smaller the degree of difference, The smaller the chromosome, the narrower the crossover range; dynamic simulation of binary crossover of two parent chromosomes generates the following offspring:

[0071]

[0072] In the formula, and These are the parameter vectors of the two parent individuals. and These are the two offspring individuals after the crossover;

[0073] After simulating binary crossover, a polynomial mutation is performed, with the mutation amount defined as follows:

[0074]

[0075] In the formula, The distribution index is the variation index when When the value is larger, the variation is smaller; when When the value is small, the variation is large; u is a random number in (0, 1);

[0076] Update individuals based on changes:

[0077]

[0078] In the formula, These are the offspring individuals resulting from dynamic simulated binary crossover. This represents the variation of a single parameter, where U and L are the upper and lower bounds for the individual, respectively, determined by the input data. This represents the individual after the change has been updated;

[0079] Perform boundary processing on the updated individuals:

[0080]

[0081] In the formula, This represents the individual after the change has been updated, where U is the upper bound of the individual and L is the lower bound of the individual. This represents the offspring individual after boundary processing, with a value range of [L, U].

[0082] By dynamically simulating binary crossover and polynomial mutation operations, this step flexibly controls the exploration range and refinement of offspring individuals in the search space based on the differences in the comprehensive fitness of the parents. When the differences between parents are small, the algorithm tends to conduct local exploration within their neighborhoods to inherit superior features and achieve refined optimization; when the differences between parents are large, it expands the range of crossover and mutation to achieve global exploration, thereby avoiding getting trapped in local optima. Combined with strict boundary handling to ensure the validity of the solution, the algorithm effectively balances its global search capability and local refinement capability during the evolutionary process, significantly enhancing population diversity and accelerating convergence to a high-performance Pareto optimal front.

[0083] Preferably, the simulated annealing local search in step (5) includes:

[0084] A subset of elite individuals are selected, and their parameters are Gaussian perturbed to constrain each parameter within a preset reasonable range. A non-dominant relationship judgment mechanism is used to compare the performance differences between the perturbated and original individuals. A binary bidding process is then used for initial screening. Individuals eliminated during the screening are accepted with probability, calculated as follows:

[0085]

[0086] In the formula, To disrupt the crowding level of individuals, Let T be the crowding level of undisturbed individuals, and T be the annealing temperature. A higher T indicates a greater probability of accepting a suboptimal solution; a lower T indicates a lower probability of accepting a suboptimal solution. The temperature cooling calculation formula is as follows:

[0087]

[0088] In the formula, This represents the annealing temperature at the k-th iteration. This represents the annealing temperature at the (k+1)th iteration. This represents the hyperparameter that controls the cooling rate.

[0089] By introducing a simulated annealing local search mechanism and conducting guided neighborhood exploration of elite individuals based on Gaussian perturbation, non-dominant relation comparison, and crowding assessment, this step can dynamically balance global broad search and local fine optimization during the evolutionary process: enhancing the ability to escape local optima and explore new regions in the high-temperature stage, focusing on the fine development of advantageous solutions in the low-temperature stage, and accepting potential inferior solutions in a probabilistic manner to maintain population diversity, thereby significantly improving the search efficiency and robustness of the algorithm in complex solution spaces, and ultimately obtaining a higher quality and better-distributed Pareto hierarchical scheme.

[0090] Preferably, the adaptive window model described in step (7) includes:

[0091] Window size decision: Calculate the local heterogeneity index HI of the geographic unit to be processed, and the calculation formula is as follows:

[0092]

[0093] In the formula, The standard deviation of all valid values ​​within the window. This is the average of all valid values ​​within the window. It is an extremely small number;

[0094] Define the adaptive window size decision function:

[0095]

[0096] In the formula, As an empirical threshold, As a local heterogeneity index, when When the region is of low heterogeneity, a large window is used; when In regions with high heterogeneity, small windows are used.

[0097] Adaptive filtering: Based on the determined window size, a corresponding window is opened with the current geographic unit as the center, the attribute values ​​of all valid geographic units in the window are counted, and the mode of the attribute value is assigned to the current geographic unit to complete the filtering and smoothing.

[0098] By employing an adaptive mechanism that dynamically adjusts the size of the filtering window based on the local heterogeneity index, this step can intelligently identify and distinguish between spatially homogeneous and heterogeneous regions: in low-heterogeneity regions, a large window is used to effectively smooth noise and enhance spatial continuity; in high-heterogeneity regions, a small window is used to finely preserve detailed features and boundary integrity; and a mode filtering strategy is used to ensure the statistical significance of the smoothed categories. Thus, while eliminating discrete geographic units, the spatial details and authenticity of the original hierarchical pattern are maintained to the maximum extent, significantly improving the spatial coherence and practical application value of the final hierarchical results.

[0099] Beneficial effects: Compared with the prior art, the present invention has the following significant advantages: 1. It takes spatial continuity as the core optimization objective and optimizes it in conjunction with the q-value of the geographic detector, so that the stratification results simultaneously meet the practical requirements of statistical significance and spatial continuity; 2. It uses a multi-objective evolutionary algorithm to automatically search for the optimal stratification scheme without manually setting weights, and efficiently provides a Pareto optimal solution set that is balanced among multiple objectives; 3. It introduces a simulated annealing mechanism to enhance the algorithm's ability to escape local optima, ensuring that the obtained solution is a globally better and more stable stratification scheme; 4. It uses adaptive filtering for post-processing, which can intelligently eliminate small patches and smooth boundaries, and finally output a zoning map with clear geographical significance and spatial continuity. Attached Figure Description

[0100] Figure 1 This is a flowchart of the method of the present invention. Detailed Implementation

[0101] The technical solution of the present invention will be further described below with reference to the accompanying drawings.

[0102] Geographical phenomena are the external forms and surface characteristics of geographical things during their occurrence, development, and change. Spatial continuity is their inherent essential attribute; the attributes of adjacent geographical units are not random and discrete, but rather naturally related. In this invention, the geographical phenomena studied are complex types formed by the combined effects of multiple factors. The multiple geographical indicators characterizing them not only exhibit complex relationships such as synergy and conflict, but also, due to the constraint of spatial continuity, show significant spatial correlation characteristics. The indicator values ​​of adjacent units often have transitional and similar characteristics, forming a continuous spatial pattern as a whole. Based on this, the geographical phenomena involved in this invention specifically include the following three types:

[0103] (1) Complex geographical phenomena, such as the synergistic occurrence of soil erosion and vegetation degradation, and the fragmentation process of wetland ecosystems;

[0104] (2) Spatial interpolation error and verification of geographical elements, such as spatial prediction bias of soil organic matter content;

[0105] (3) Analysis and evaluation of human geography complex systems, such as urban planning, ecological zoning, and environmental governance.

[0106] In this invention, explanatory factors are geoscientific elements capable of explaining the spatial differentiation of the studied geographical phenomenon; geographical indicators refer to geoscientific elements capable of characterizing the studied geographical phenomenon. For example, if the incidence of respiratory diseases in a certain area is selected as the phenomenon to be studied, potential driving factors such as temperature, precipitation, vegetation index, DEM, GDP, surface temperature, population density, and traffic density are selected as explanatory factors, and air pollution and health problems are selected as indicators to characterize the incidence of respiratory diseases, that is, the spatial differentiation of the incidence of respiratory diseases characterized by air pollution and health problems is explained by explanatory factors such as temperature, precipitation, vegetation index, DEM, and GDP.

[0107] Example 1

[0108] like Figure 1 As shown, this embodiment provides a multi-objective coupled geographic layering method based on spatial continuity constraints, including the following steps:

[0109] (1) For specific geographical phenomena in the target geographical area, select appropriate explanatory factors and multiple geographical indicators that characterize the phenomenon, construct a dynamic mask for the explanatory factors and target factors and prune them to obtain a geographical dataset without null values.

[0110] Dynamic masks are defined as:

[0111]

[0112] In the formula, i and j represent the row and column numbers, respectively, and the number of rows and columns is defined as... Construct a mask array ;

[0113] (2) Randomly generate an initial population containing several individuals, each defined by a set of hyperplane parameters; construct hyperplane equations using hyperplane parameters and explanatory factors, each hyperplane dividing geographic units within a geographic region into two types of strata, and each individual corresponding to a regional stratification scheme;

[0114] The parameters of the hyperplane constituting each individual are defined as follows:

[0115]

[0116] In the formula, j is the j-th hyperplane. Let be the hyperplane parameters corresponding to the i-th explanatory factor. For constant terms;

[0117] For the i-th pixel, its discriminant function on the j-th hyperplane is:

[0118]

[0119] The classification rules are as follows:

[0120]

[0121] In the formula, To determine the category of the i-th pixel in the j-th hyperplane, calculate the category of the i-th pixel in all hyperplanes, and finally obtain... The category label, h represents the number of hyperplanes.

[0122] (3) Calculate the geodetector q value and spatial continuity index of each geographic indicator using the stratification results of the target geographic region, and use the two together as the objective function value of the individual; screen the initial population by combining non-dominated sorting with crowding sorting to obtain the parent population.

[0123] (31) Calculate the geodetector q-value for each geographic indicator:

[0124]

[0125] In the formula, Let q represent the geodetector q value obtained by classifying the data based on the hyperplane parameters of the i-th individual, with a value range of [0,1]. The closer a value is to 1, the stronger its explanatory power for geographical phenomena within the region.

[0126] (32) Calculate the spatial continuity index under each stratification scheme:

[0127] The adjacency homogeneity index is defined as:

[0128]

[0129] In the formula, homogeneity means that two adjacent pixels belong to the same category.

[0130] Calculate the block morphology index for each individual:

[0131]

[0132] In the formula, P represents the perimeter of the block, and A represents the area of ​​the block.

[0133] (33) The objective function is designed as follows:

[0134]

[0135] In the formula, This represents the clustering objective of the nth y-factor. Represents the adjacency homogeneity index. This indicates the blocky morphology index.

[0136] (33) Non-dominated sorting includes:

[0137] For any two individuals a and b in the initial population, compare them based on n geographical indicators:

[0138]

[0139] In the formula, Let represent the target vector of individual a on the i-th geographical indicator. Let this be the target vector of individual b on the i-th geographic indicator; if for any geographic indicator, All greater than or equal to A is said to dominate B.

[0140] (34) Crowding ranking includes:

[0141] When selecting a set number of dominant individuals, if the number of individuals in a certain layer exceeds the number of dominant individuals to be selected in that layer, the crowding degree of each individual in that layer is calculated:

[0142]

[0143] In the formula, Let m be the crowding degree of the j-th individual, and m be the initial total number of individuals in the population. This represents the geodetector q-value of the x-th individual on the ith geographic index;

[0144] (35) The binary tournament selects offspring individuals:

[0145]

[0146] In the formula, a.rank is the non-dominated ranking of individual a, and b.rank is the non-dominated ranking of individual b. Let 'a' represent the crowding level of individual 'a'. Sort individuals b by crowding. If individual a has a better non-dominated ranking, retain individual a directly; if individuals a and b belong to the same non-dominated ranking level, use crowding as a secondary criterion, prioritizing individuals with lower crowding values ​​to enter the subsequent population, and gradually select offspring individuals with both better ranking and lower crowding through pairwise competition logic.

[0147] Based on non-dominant sorting and crowding sorting, a binary tournament is used to rank the individuals in the population as a whole to obtain the population ranking relationship. Dominant individuals are selected in turn until the number of dominant individuals reaches a set number.

[0148] (4) Direction-aware simulated binary crossover and polynomial mutation operations are performed on the hyperplane parameters of the parent individuals in the parent population to generate offspring individuals. First, the distribution parameters of the simulated binary crossover are dynamically adjusted according to the degree of difference between the parent individuals and the multi-objective function. When the difference between the parent individuals is large, the crossover coefficient is increased to expand the search range and enhance the global exploration capability. When the difference between the parent individuals is small, the crossover coefficient is decreased to strengthen the local development capability, thereby avoiding getting trapped in local optima while maintaining the inheritance of excellent genes. Subsequently, in the polynomial mutation stage, a random probability following the interval [0,1] is generated for each set of hyperplane parameters to determine whether to perform the mutation operation. This allows for local development near the parent to inherit excellent characteristics, and also enables global exploration through larger mutations to escape local optima. Through the above mechanism of combining difference-driven crossover and polynomial mutation, adaptive exploration and fine convergence of the solution space are achieved, ensuring that the offspring individuals have multi-objective performance characteristics.

[0149] (5) At regular intervals of a certain number of iterations, a local reconstruction operation based on the annealing mechanism is triggered. Individuals with high fitness are selected from the population as candidate solutions, and their hyperplane parameters are subjected to Gaussian perturbation. The perturbation results are restricted to a preset reasonable range to avoid exceeding the limit. The relative merits of the perturbed individuals and the original individuals are evaluated using non-dominated sorting rules. Initial screening is performed using binary bidding rules. Individuals eliminated in the screening are accepted with probability, calculated as follows:

[0150]

[0151] In the formula, To disrupt the crowding level of individuals, The density of undisturbed individuals is represented by T, which is the annealing temperature. A larger T indicates a higher probability of accepting a poorly conceived solution, leading to broader exploration. Conversely, a smaller T indicates a lower probability of accepting a poorly conceived solution, favoring local refinement. The temperature cooling calculation formula is as follows:

[0152]

[0153] In the formula, Both T and k are hyperparameters, and k is the number of iterations.

[0154] (6) Merge the parent and child individuals into a population and perform screening to obtain a new parent population, thereby ensuring that the optimal solution is not lost and improving the convergence speed of the algorithm.

[0155] (7) Determine whether the maximum number of iterations has been reached. If it has, stop the iteration and output the layered results of the target geographic region.

[0156] (8) Use an adaptive window model to further optimize the output hierarchical results to ensure the spatial continuity of the final result data.

[0157] Adaptive window models include:

[0158] Window size decision: Calculate the local heterogeneity index HI of the geographic unit to be processed, which is the mean absolute difference of the 8 neighboring pixel values. The calculation formula is as follows:

[0159]

[0160] The geographic unit is the geographic element to be classified within the target geographic area, and the data type is raster data;

[0161] Define the adaptive window size decision function:

[0162]

[0163] In the formula, As an empirical threshold, As a local heterogeneity index, when In regions with low heterogeneity, a large window is used to enhance isotropic smoothing; when In regions with high heterogeneity, small windows are used to better preserve edge features;

[0164] Adaptive filtering: Based on the determined window size, a corresponding window is opened with the current geographic unit as the center, the attribute values ​​of all valid geographic units in the window are counted, and the mode of the attribute value is assigned to the current geographic unit to complete the filtering and smoothing.

[0165] The following is a detailed explanation of a specific implementation scenario, taking Beijing as the target geographical area and the urban heat island effect as the geographical phenomenon to be studied. We select elevation, urban building height, and nighttime light index as explanatory factors, and use the temperature data of three pressure layers (780 Pa, 850 Pa, and 920 Pa) from the MOD07 atmospheric profile product as three geographical indicators.

[0166] Using geographical elevation (DEM), nighttime light index, and urban building height as three explanatory factors, this study explains the spatial differentiation characteristics of temperature across different pressure layers in Beijing, thereby revealing the driving role of each factor in the spatial pattern of the urban heat island effect. Both constrained and unconstrained models were employed. Based on the above data, geographic stratification was performed, and the corresponding q-values ​​of the geographic detectors for different numbers of hyperplanes are shown in the table below.

[0167] Table 1 Comparison of Results for the Three Hyperplanes

[0168] 780pa 850pa 920pa Unconstrained model q-value 0.40351 0.48254 0.49246 Constraint model q-value 0.41199 0.49814 0.50781 Target value growth percentage 2.10% 3.23% 3.12%

[0169] Table 2 Comparison of Results for the Four Hyperplanes

[0170] 780pa 850pa 920pa Unconstrained model q-value 0.43699 0.52442 0.52768 Constraint model q-value 0.45664 0.53342 0.53538 Target value growth percentage 4.50% 1.72% 1.46%

[0171] In both constrained and unconstrained models, the geospatial detector q-values ​​for 920 hPa are the largest, while those for 780 hPa are the smallest. This result is highly consistent with actual geographical patterns. From a physical perspective, the temperature data from the 920 hPa barylayer, being closer to the surface, exhibits a more significant correlation with surface environmental factors. It shows higher correlations with topographic relief reflected by the DEM, the intensity of human activity represented by the nighttime light index, and the land cover characteristics reflected by building height, thus providing a more prominent explanatory power for geographic stratification. In contrast, the 780 hPa barylayer, being farther from the ground, is less affected by surface factors, and its explanatory power is naturally lower.

[0172] As shown in Tables 1 and 2, a further comparison of the two models reveals that, regardless of whether there are 3 or 4 hyperplanes, the constrained model consistently outperforms the unconstrained model in terms of q values ​​at the three pressure levels of 780 hPa, 850 hPa, and 920 hPa. This demonstrates that by introducing spatial constraints, the constrained model can effectively strengthen the correlation between explanatory factors and geographical indicators, thereby enhancing the rationality and effectiveness of geographical stratification.

[0173] The objective q-values ​​of the three hyperplane-constrained models were improved by approximately 2.82% compared to the unconstrained models. The results show that the constrained models outperformed the unconstrained models in all pressure scenarios, and the optimization effect showed a positive correlation with pressure. The lower the pressure, the smaller the improvement in the q-value of the constrained models; the higher the pressure, the more significant the improvement. This pattern is essentially determined by the strength of the correlation between the geographical characteristics of each pressure layer and surface factors. Under high-pressure environments, the correlation between geographical phenomena and surface explanatory factors is more complex, and the constraints can more specifically strengthen this correlation, thus achieving a more significant improvement in explanatory power. Under low-pressure environments, the correlation is weaker, and the optimization space is relatively limited. This data fully validates the actual effectiveness of the constraints, provides a basis for model selection and optimization of the explanatory power of geographical phenomena under different pressure scenarios, and highlights the unique value of constrained models in uncovering complex geographical correlations.

[0174] In the four hyperplane scenarios, the target q-value of the constrained model was improved by an average of approximately 2.56% compared to the unconstrained model. Furthermore, the improvement gradually decreased from the 780 hPa to the 920 hPa pressure layer, a phenomenon highly consistent with the classification optimization pattern resulting from the increased number of hyperplanes. The 780 hPa pressure layer, being far from the ground, is less affected by surface explanatory factors, leading to a lower base q-value for the unconstrained model at this layer, indicating greater optimization potential. The constrained model, by introducing spatial constraints, can more accurately uncover the potential correlation between explanatory factors and the temperature of this pressure layer, thus achieving a more significant q-value improvement. The 920 hPa pressure layer, being closer to the surface, has a stronger correlation with surface explanatory factors. The base q-value of the unconstrained model at this layer is already at a high level, approaching the theoretical optimal state of geographic stratification explanatory power, leaving relatively limited optimization space. Therefore, the q-value improvement brought by the constrained model is less than that at the 780 hPa layer. Although the increase in the number of hyperplanes leads to an increase in the number of classification categories, the difference in the optimization effect of the constraint model on different pressure layers is essentially determined by the correlation strength between the geographical characteristics of each pressure layer and surface factors. Among them, the improvement advantage of the 780hPa pressure layer, which is "weakly correlated and has low basic explanatory power", is more prominent, which fully reflects the unique value of the constraint model in exploring the spatial differentiation patterns of weakly correlated geographical phenomena.

[0175] Example 2

[0176] This embodiment provides a multi-objective coupled geographic layering system based on spatial continuity constraints, corresponding to the multi-objective coupled geographic layering method based on spatial continuity constraints in Embodiment 1 above, including:

[0177] Data preprocessing module: For the target study area and the geographical phenomena to be analyzed, at least two indicator variables and one or more explanatory factors are selected to describe the characteristics of the geographical response. A dynamic mask is constructed to remove null values ​​from the multi-source factor data, resulting in a geographical dataset without null values.

[0178] Spatial stratification construction module: Randomly generates an initial population containing N individuals, each individual being encoded as a set of hyperplane parameters. Each set of parameters can construct multiple hyperplanes. Based on the spatial segmentation logic formed by the hyperplanes, each geographic unit within the target geographic region is assigned to the corresponding stratum, and each individual corresponds to an independent set of candidate stratification schemes;

[0179] Population optimization and screening module: For each stratification scheme, the q-value of each geographic indicator under the stratum is calculated. Combined with the adjacency homogeneity index and the block morphology index, the explanatory power and spatial continuity of the stratification scheme are comprehensively reflected. The objective function of multi-objective optimization is constructed to perform non-dominated ranking and crowding calculation on individuals in the current population, and the parent individuals are screened based on the binary tournament.

[0180] Parameter evolution module: Constructs evaluation indicators, dynamically adjusts the crossover range of simulated binary crossover based on the differences between parent individuals, and generates preliminary offspring individuals; then performs polynomial perturbation on the hyperplane parameter vector of the offspring individuals piecewise to generate offspring individuals;

[0181] Iterative optimization module: Merges the parent and offspring populations to form a new composite population, performs a selection process to select a new generation of parent populations from the composite population; introduces simulated annealing local search during the evolution process, performs corresponding neighborhood perturbation on some elite individuals, and adopts a global optimal solution retention mechanism.

[0182] Output and Judgment Module: Used to repeatedly call the iterative optimization module and determine whether the maximum number of iterations has been reached. If yes, terminate the iteration and output the region stratification result that satisfies the spatial continuity constraint; otherwise, return to the parameter evolution module to continue iterating.

[0183] Results Optimization Module: Utilizes a dynamically adjustable window analysis model to optimize the final results. Larger windows are used to enhance spatial cohesion in areas where similar land features are distributed in contiguous patches, while smaller windows are used at type boundaries to retain necessary detail differences, ensuring that the spatial distribution pattern of the final results closely matches the natural differentiation patterns of actual geographical phenomena.

Claims

1. A multi-objective coupled geographic layering method based on spatial continuity constraints, characterized in that, Includes the following steps: (1) For the target study area and the geographical phenomena to be analyzed, select at least two indicator variables and one or more explanatory factors to describe the characteristics of geographical response. By constructing a dynamic mask, the indicator variables and explanatory factors are removed to obtain a geographical dataset without null values. (2) Randomly generate an initial population containing N individuals. Each individual is encoded with a set of hyperplane parameters. Multiple hyperplanes are constructed based on each set of hyperplane parameters. According to the spatial segmentation logic corresponding to the hyperplane, each geographic unit in the target study area is divided into the corresponding stratum. Each individual corresponds to an independent candidate stratification scheme. (3) For each candidate stratification scheme, calculate the geospatial detector q value of each indicator variable under the stratification scheme, combine the adjacency homogeneity index and block morphology index to construct a multi-objective optimization function, perform non-dominated sorting and crowding calculation on individuals in the current population, and select parent individuals based on binary tournament to form the parent population. (4) For the selected parent individuals, perform dynamic simulated binary crossover, and adaptively adjust the crossover range of simulated binary crossover according to the difference in objective function between parent individuals to generate preliminary offspring individuals; A polynomial perturbation is performed on the hyperplane parameter vector corresponding to the initial offspring individuals to generate offspring individuals and form an offspring population; (5) Merge the parent population and the offspring population to form a new composite population; select a new generation of parent individuals from the composite population based on a binary tournament; introduce simulated annealing local search, select the top K elite individuals in the non-dominated sorting to perform corresponding neighborhood perturbations, and adopt a global optimal solution retention mechanism. (6) Repeat steps (3) to (5). When the number of iterations reaches the preset maximum number of iterations, terminate the evolution and output the optimal region stratification result from the non-dominated sorting result of the final generation. (7) An adaptive window model is used to perform spatial smoothing on the regional stratification results of the iterative output to eliminate discrete geographic units; The formula for calculating the q-statistic value of the geographic detector mentioned in step (3) is as follows: In the formula, Let q represent the geodetector q value obtained by classifying the data based on the hyperplane parameters of the i-th individual, with a value range of [0,1]. The closer a value is to 1, the stronger the explanatory power of the classification for geographical phenomena within the region. Calculate the spatial continuity index under each stratification scheme: Adjacent cells in a raster are compared. If two adjacent cells belong to the same category, they are considered homogeneous; otherwise, they are considered heterogeneous. The Adjacency Homogeneity Index (AH) is defined as follows: The range of AH is [0,1]. The closer AH is to 1, the more adjacent pixels belong to the same category, that is, the stronger the spatial continuity. The closer AH is to 0, the greater the difference between adjacent pixels, that is, the more fragmented the space. Based on the boundary complexity, area size, and shape compactness of the blocks, the ability of similar pixels to form larger, regular, continuous blocks is determined using the Block Morphology Index (BI), calculated as follows: In the formula, P represents the perimeter of the block, and A represents the area of ​​the block.

2. The multi-objective coupled geographic layering method according to claim 1, characterized in that, The dynamic mask mentioned in step (1) is defined as follows: In the formula, i and j represent the row and column numbers, respectively, and the number of rows and columns is defined as... Construct a mask array .

3. The multi-objective coupled geographic layering method according to claim 1, characterized in that, The hyperplane parameters for each individual mentioned in step (2) are defined as follows: In the formula, j is the j-th hyperplane. Let b be the hyperplane parameter corresponding to the i-th explanatory factor, and b be a constant term. For the i-th pixel, its discriminant function on the j-th hyperplane is: In the formula, This represents the value of the nth explanatory factor in the i-th geographic unit; The classification rules are as follows: In the formula, To determine the category of the i-th pixel in the j-th hyperplane, calculate the category of the i-th pixel in all hyperplanes, and finally obtain... The category label, h represents the number of hyperplanes.

4. The multi-objective coupled geographic layering method according to claim 1, characterized in that, The objective function described in step (3) is designed as follows: In the formula, This represents the clustering target value of the nth y-factor. Represents the adjacency homogeneity index. Indicates the mass morphology index. This represents the value of a multi-objective optimization function, used to comprehensively reflect the overall performance of candidate stratification schemes in terms of both geographical explanatory power and spatial continuity.

5. The multi-objective coupled geographic layering method according to claim 1, characterized in that, The non-dominated sorting described in step (3) includes: Compare any two individuals in the initial population based on n geographical indicators: In the formula, Let represent the target vector of individual a on the i-th geographical indicator. Let this be the target vector of individual b on the i-th geographic indicator; if for any geographic indicator, All greater than or equal to a is said to dominate b; if for some target i, there exists Meanwhile, other objectives also exist. If the dominance relationship is not established, it indicates that the two individuals have their own advantages in different objectives and cannot establish a dominance relationship. Therefore, the relationship between individual a and individual b is determined to be non-dominant. By performing the above dominance determination process on any two individuals in the population, the dominance matrix of the entire population is obtained. The system stratifies all individuals according to the dominance matrix to form different non-dominant frontier hierarchical structures. The congestion calculation includes: In the process of selecting a preset number of dominant individuals by level, if the total number of individuals in a certain level exceeds the planned quota of dominant individuals for that level, then the crowding degree of all individuals within that level is quantitatively calculated: In the formula, Let m be the crowding degree of the j-th individual, and m be the initial total number of individuals in the population. Let q represent the geodetector value of the x-th individual on the i-th geographic indicator.

6. The multi-objective coupled geographic layering method according to claim 1, characterized in that, After the non-dominated ranking and crowding calculation described in step (3), the process also includes screening offspring individuals using a binary tournament: In the formula, a.rank represents the non-dominant ranking rank of individual a, and b.rank represents the non-dominant ranking rank of individual b. Let 'a' represent the crowding level of individual 'a'. Rank individuals b by their crowding level; The comparison function is as follows: if individual a has a better non-dominated ranking level, then individual a is directly retained; if individual a and individual b belong to the same non-dominated ranking level, then crowding is used as a secondary criterion, and the individual with the smaller crowding value is selected to enter the subsequent population. Through pairwise competition logic, offspring individuals with both better ranking and lower crowding are gradually selected.

7. The multi-objective coupled geographic layering method according to claim 1, characterized in that, The calculation steps for crossover variation described in step (4) are as follows: By aggregating multiple objective functions, a parent individual evaluation index is constructed, and the geographical explanatory power index Q is: In the formula, This represents the clustering target value of the first y-factor. This represents the clustering target value of the second y-factor. This represents the clustering target value for the third y-factor; The spatial continuity index S is aggregated as follows: In the formula, Represents the adjacency homogeneity index. The weights for the adjacency homogeneity index; Indicates the mass morphology index. The weights for the adjacency homogeneity index; The comprehensive evaluation index W is: In the formula, Q represents the geographical explanatory power index. The weights for geographical explanatory power indicators; S is an index of spatial continuity. The weights for the spatial continuity index; The multi-objective differences among parent individuals are: In the formula, This represents the comprehensive evaluation index of the stratification scheme corresponding to parent individual i. This represents the comprehensive evaluation index of the stratification scheme corresponding to the parent individual j. The D value represents the difference in overall fitness between two individuals; the larger the D value, the greater the difference in their overall performance across multiple objectives; the smaller the D value, the more similar their performance, which can be used for local refinement. Based on the degree of difference Adaptive adjustment of cross-range control parameters The greater the difference, The larger the value, the wider the range of intersection; the smaller the degree of difference, The smaller the chromosome, the narrower the crossover range; dynamic simulation of binary crossover of two parent chromosomes generates the following offspring: In the formula, and These are the parameter vectors of the two parent individuals. and These are the two offspring individuals after the crossover; After simulating binary crossover, a polynomial mutation is performed, with the mutation amount defined as follows: In the formula, The distribution index is the variation index when When the value is larger, the variation is smaller; when When the value is small, the variation is large; u is a random number in (0, 1); Update individuals based on changes: In the formula, These are the offspring individuals resulting from dynamic simulated binary crossover. This represents the variation of a single parameter, where U and L are the upper and lower bounds for the individual, respectively, determined by the input data. This represents the individual after the change has been updated; Perform boundary processing on the updated individuals: In the formula, This represents the individual after the change has been updated, where U is the upper bound of the individual and L is the lower bound of the individual. This represents the offspring individual after boundary processing, with a value range of [L, U].

8. The multi-objective coupled geographic layering method according to claim 1, characterized in that, The simulated annealing local search described in step (5) includes: A subset of elite individuals are selected, and their parameters are Gaussian perturbed to constrain each parameter within a preset reasonable range. A non-dominant relationship judgment mechanism is used to compare the performance differences between the perturbated and original individuals. A binary bidding process is then used for initial screening. Individuals eliminated during the screening are accepted with probability, calculated as follows: In the formula, To disrupt the crowding level of individuals, Let T be the crowding level of undisturbed individuals, and T be the annealing temperature. A higher T indicates a greater probability of accepting a suboptimal solution; a lower T indicates a lower probability of accepting a suboptimal solution. The temperature cooling calculation formula is as follows: In the formula, This represents the annealing temperature at the k-th iteration. This represents the annealing temperature at the (k+1)th iteration. This represents the hyperparameter that controls the cooling rate.

9. The multi-objective coupled geographic layering method according to claim 1, characterized in that, The adaptive window model described in step (7) includes: Window size decision: Calculate the local heterogeneity index HI of the geographic unit to be processed, and the calculation formula is as follows: In the formula, The standard deviation of all valid values ​​within the window. This is the average of all valid values ​​within the window. It is an extremely small number; Define the adaptive window size decision function: In the formula, As an empirical threshold, As a local heterogeneity index, when When the region is of low heterogeneity, a large window is used; when In regions with high heterogeneity, small windows are used. Adaptive filtering: Based on the determined window size, a corresponding window is opened with the current geographic unit as the center, the attribute values ​​of all valid geographic units in the window are counted, and the mode of the attribute value is assigned to the current geographic unit to complete the filtering and smoothing.

Images