Method for identifying floating population settlement area based on remote sensing image

Abstract

The application relates to the technical field of remote sensing image processing, and particularly discloses a method for identifying a floating population settlement area based on remote sensing images, which comprises the following steps: fusing a time-series Sentinel-2 spectral angle mapping and a night light convolutional autoencoder, separating the spectral shift features of day-time temporary building materials and the spatial distribution of night-time low-power consumption lights in a target area, and outputting a temporary settlement preliminary screening area; introducing a GRU trajectory prediction and a generative adversarial network in the temporary settlement preliminary screening area, simulating the decoupling of space-time dependence of travel chains and jobs and residences, and judging the crowd flow stability and the stay mode through mode classification; the application separates the spectral shift features of day-time temporary building materials by using the multi-temporal spectral angle time-series variance and simultaneously separates the night-time low-power consumption cluster light distribution through residual reconstruction of the autoencoder by fusing the time-series Sentinel-2 spectral angle mapping and the night light convolutional autoencoder, realizes dynamic separation of the two types of ground objects, and significantly improves the accuracy and robustness of the temporary settlement preliminary screening.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image processing technology, and in particular to a method for identifying migrant population clusters based on remote sensing images. Background Technology

[0002] With economic development and accelerated urbanization, a large number of migrant workers have formed specific residential clusters in cities. These areas are of great significance in terms of social management, public service allocation, and urban planning. Migrant worker clusters are usually characterized by dense spatial distribution, diverse population structure, relatively concentrated and rapidly changing living conditions, making it difficult for traditional population statistics and management methods to reflect their distribution in real time and accurately.

[0003] For example, a method for identifying Wa ethnic minority residential areas, published in Chinese Patent Publication No. CN109190585A, provides a scientific and reasonable solution for distinguishing Wa ethnic minority settlements. This not only saves human resources and greatly improves work efficiency, but also helps to protect the cultural characteristics of ethnic minority areas.

[0004] In existing technologies, the spectral characteristics of temporary residences and permanent buildings of the floating population are highly overlapping in day and night remote sensing images, making it difficult to separate the daytime reflection signals of temporary building materials from the distribution of low-power lighting at night from static features. Moreover, within the initially separated temporary residential area, there is a lack of evidence to distinguish the behavioral patterns of high-density short-term renters and low-density nighttime residents. A single image cannot reflect the spatiotemporal dependence of their work and residence and the differences in their mobility stability. Furthermore, since the obtained building type features and population behavior features are independent of each other, there is a lack of a mechanism for collaboratively determining the boundaries of the real settlement area. This leads to conflicts between the spatial distribution of buildings and the dynamic patterns of behavior during the boundary generation process, making it impossible to form a unified identification decision. Therefore, this paper proposes a method for identifying floating population settlements based on remote sensing images to solve the above-mentioned problems. Summary of the Invention

[0005] To overcome the shortcomings of the prior art, the present invention provides a method for identifying migrant population settlements based on remote sensing images, which can effectively solve the problems mentioned in the prior art.

[0006] The objective of this invention can be achieved through the following technical solution: This invention provides a method for identifying migrant population clusters based on remote sensing imagery, comprising the following steps:

[0007] Step 1: Integrate temporal Sentinel-2 spectral angle mapping with nighttime light autoencoder to separate the spectral drift characteristics of daytime temporary building materials and the spatial distribution of nighttime low-power lights in the target area, outputting the initial screening area of ​​temporary housing, effectively separating temporary and permanent buildings, and improving the accuracy of the initial screening.

[0008] Step 2: Introduce GRU trajectory prediction and generative adversarial network in the initial screening area of ​​temporary residence area to simulate the travel chain to decouple the spatiotemporal dependence of work and residence. Determine the stability of population flow and residence pattern through pattern classification, decouple work and residence dependence, and finely classify population behavior patterns.

[0009] Step 3: Construct a dual-agent deep Q-network to enable the building feature agent and the behavior feature agent to make decisions alternately, and integrate the local constraints of the conditional random field to generate the initial settlement boundary, coordinate building and behavior features, and eliminate misjudgment of single feature boundaries.

[0010] Step 4: The two agents iteratively adjust their weights based on the initial settlement boundary feedback, eliminate boundary conflicts between building space and behavioral dynamics through a joint reward mechanism, and iteratively eliminate boundary conflicts to achieve collaborative judgment based on dual features.

[0011] Step 5: Utilize the difference between the pseudo-trajectory and real trajectory distribution output by the generative adversarial network to quantify the stability of the mobile population, correct the division of high-density and low-density areas within the boundary and form a stability-level patch boundary, quantify the flow stability, and achieve refined classification of settlement areas.

[0012] Step 6: Globally optimize boundary smoothness based on conditional random field, integrate the final decision of the two agents, output the spatial range and stability classification results of the settlement area, globally optimize boundary smoothness, and output the connectivity classification map.

[0013] Preferably, step 1 specifically includes:

[0014] We acquired multi-temporal Sentinel-2L2A surface reflectance image sequences of the target area, constructed a spectral reference library for temporary building materials, calculated the temporal variance of spectral angles pixel by pixel, extracted the spectral drift feature map caused by daytime temporary building materials, and effectively distinguished the spectral differences between temporary building materials and permanent buildings.

[0015] Simultaneously acquire low-power nighttime lighting images, use a convolutional autoencoder to separate the response of permanent buildings with strong light sources, and extract the distribution map of low-intensity clustered lights through decoder residual reconstruction to accurately separate low-power temporary lights from permanent buildings with strong light sources.

[0016] The spectral drift feature map and the low-intensity clustered light distribution map are subjected to pixel-level logical AND operation. After morphological closing operation and area threshold filtering, the initial screening area of ​​temporary residential area is output. The day and night features are integrated to remove static ground object interference and improve the purity of the initial screening area.

[0017] Preferably, step 2 specifically includes:

[0018] Individual travel time series are extracted from the initial screening area of ​​temporary residential areas. Using trajectories from multiple consecutive days as input, the temporal dependency relationship is encoded by gated cyclic unit encoding to predict the sequence of residence, workplace and stopover point in the travel chain of the next day. This effectively captures the temporal dependency of individual travel and improves the prediction accuracy of residence and workplace.

[0019] By analyzing the error distribution between the predicted trajectory and the actual trajectory, the spatiotemporal dependence of work and residence is decoupled, the spatiotemporal regularity of the round trip from residence to workplace and the frequency of daytime stop migration are calculated, the degree of work-residence binding is quantified, and regular commuting behavior is distinguished from random mobility behavior.

[0020] Generative adversarial networks are used to generate pseudo-trajectories to expand the samples. The prediction error of GRU and the job-housing dependence index are combined into a feature vector for subsequent classification of population flow stability and residence patterns, which enriches the diversity of training samples and enhances the generalization ability of flow pattern classification.

[0021] Preferably, step 2 further includes:

[0022] Based on the intensity of work-residence spatiotemporal dependence, individuals are classified into work-residence bound type and work-residence free type. The work-residence bound ratio and daytime work relocation rate of each grid in the initial screening area of ​​temporary housing are calculated to quantify the degree of dependence of grid population on external employment centers.

[0023] Based on the kurtosis and skewness of the distribution of continuous nighttime stay duration, high-density short-term rental patterns and low-density nighttime stay patterns are divided, and the dominant behavior categories of each grid are marked to distinguish between frequent rental rotation and temporary nighttime stay.

[0024] The stability of the flow is quantified as the temporal variation coefficient of the trajectory daytime similarity, which together with the dwelling pattern category constitutes a grid-level behavioral feature label, realizing the joint representation of the intensity of population flow and dwelling pattern.

[0025] Preferably, step 3 specifically includes:

[0026] Construct an intelligent agent for building features. The state space consists of the spectral angle drift intensity, temporary building material probability, and local texture features of each pixel. The action space consists of binary classification of whether the pixel belongs to the boundary of a residential area or not, thus achieving autonomous classification of pixel-level building features.

[0027] A behavioral feature intelligent agent is constructed, with the state space consisting of grid-level flow stability indicators, dwelling pattern probabilities, and job-housing dependence strength. It adopts the same action space and forms a parallel decision-making structure with the building intelligent agent to achieve parallel discrimination of grid-level behavioral features.

[0028] The design employs an alternating decision-making strategy, in which the building agent prioritizes outputting the building-side classification result in each round, and the behavioral agent corrects the behavioral anomaly grid based on this. The two alternately iterate until local convergence, and then output a classification probability map to eliminate boundary misjudgments driven by a single feature.

[0029] Preferably, step 3 further includes:

[0030] The classification probability maps of the alternating decision-making outputs of the two agents are used as the unary potential energy of the conditional random field. The building potential energy comes from the output of the building agent, and the behavior potential energy comes from the correction result of the behavior agent, effectively integrating the two heterogeneous features of building and behavior.

[0031] We construct a binary potential energy for a conditional random field to penalize label inconsistencies between adjacent pixels or grids. At the same time, we introduce spectral similarity and behavioral consistency constraints within the local neighborhood to suppress local noise and enhance the continuity of the boundary space.

[0032] The mean-field inference algorithm is used to solve the conditional random field. While maintaining the local consistency between buildings and behaviors, it outputs a smooth initial settlement boundary that meets local constraints, generating a smooth boundary that takes into account both spectral and behavioral consistency.

[0033] Preferably, step 4 specifically includes:

[0034] The joint reward function is defined as the harmonic average of the entropy of building features and the entropy of behavior features within the initial settlement boundary. A high reward is given when both are highly consistent, and a low reward is given when there is conflict between building and behavior, which effectively quantifies the degree of synergy between building and behavior.

[0035] The two agents update their respective deep Q-network weights based on joint rewards, record the difference in their contributions to boundary generation in each iteration, reduce the state fusion weight of the agent with excessive contribution, and dynamically balance the decision influence of the two agents.

[0036] Repeatedly alternating decision-making and conditional random field boundary generation until the joint reward no longer increases after multiple rounds, outputting the settlement boundary after dynamic collaborative determination of building space and behavior, ensuring that the boundary converges to the optimal collaborative state.

[0037] Preferably, step 5 specifically includes:

[0038] Within the boundaries of the settlement area after collaborative determination, the distribution of real travel trajectories of each grid is extracted, and the trained generative adversarial network is run to generate a corresponding number of pseudo trajectory distributions, so as to achieve benchmark alignment between real behavior and generated patterns.

[0039] Calculate the Wasserstein distance between the real trajectory distribution and the pseudo trajectory distribution within each grid, as a quantitative indicator of trajectory pattern complexity. The larger the distance, the more difficult it is to predict the actual behavior, thus quantifying the complexity and randomness of the behavior.

[0040] The Wasserstein distance is normalized and mapped to a stability index of the mobile population. The higher the index, the worse the stability, and the lower the index, the more regular and stable the clustering behavior, thus constructing a unified metric for population stability.

[0041] Preferably, step 5 further includes:

[0042] The stability index of the mobile population in each grid within the original high-density short-term rental area is compared with a preset high stability threshold. Grids below the threshold are re-labeled as stable settlements, while those above the threshold remain high-density mobile areas. Grids with regular behavior patterns within the high-density short-term rental area are identified as stable settlements.

[0043] The stability index of the mobile population in each grid within the original low-density nighttime stay area is compared with a preset low stability threshold. Grids with values ​​higher than the threshold are relabeled as extremely mobile temporary areas, while those with values ​​lower than the threshold remain low-density stable areas. Grids with random behavior in the low-density nighttime stay area are classified as extremely mobile temporary areas.

[0044] The corrected four-category region division results are output to form the boundary of stability classification patches, covering high-density stable areas, high-density flow areas, low-density stable areas and low-density flow areas. The four types of stability classification patches are output to refine the behavioral characteristics of settlement areas.

[0045] Preferably, step 6 specifically includes:

[0046] The boundary of the settlement after the collaborative judgment of the final decision output of the two agents is used as the unary potential energy of the conditional random field, and the boundary of the output four types of stability graded patches is used as the regional connectivity constraint of the binary potential energy to ensure that the boundary integrates the dual prior information of architecture and behavior.

[0047] An isolated region with an area smaller than the preset minimum connected patch is penalized by a higher-order potential energy, and is forcibly merged into the adjacent dominant category. A graph cut algorithm is used to globally optimize the boundary smoothness and spatial connectivity, eliminate noisy patches, and improve the continuity and integrity of the boundary space.

[0048] Outputs the globally optimized spatial extent vector boundary of the settlement, along with a stability classification distribution map of each settlement, including the area proportion and average flow stability index of each classification, providing refined attribute data that can be directly used for spatial analysis and mapping.

[0049] Compared with the prior art, the beneficial effects of the present invention are:

[0050] 1. This method for identifying migrant population settlements based on remote sensing imagery integrates temporal Sentinel-2 spectral angle mapping with a nighttime light convolutional autoencoder. It utilizes the temporal variance of multi-phase spectral angles to extract the spectral drift features of daytime temporary building materials. Simultaneously, it separates the distribution of low-power clustered lights at night through autoencoder residual reconstruction, achieving dynamic separation of the two types of land features and significantly improving the accuracy and robustness of the initial screening of temporary settlements.

[0051] 2. This method for identifying migrant population clusters based on remote sensing imagery introduces GRU trajectory prediction and generative adversarial network into the initial screening area. It decouples the spatiotemporal dependence of work and residence by simulating individual travel chains. At the same time, it divides the behavior patterns into high-density short-term rental and low-density nighttime stay based on the kurtosis and skewness of the nighttime stay distribution. It also combines the coefficient of variation of trajectory daytime similarity to quantify the stability of migration. This method overcomes the limitation that a single image cannot reflect the differences in population behavior patterns and provides multi-dimensional behavioral characteristics for determining the boundaries of the clusters.

[0052] 3. This method for identifying migrant population clusters based on remote sensing imagery constructs a dual-agent deep Q-network, enabling the building feature agent and the behavior feature agent to process pixel-level building spectral features and grid-level population behavior features respectively through an alternating decision-making mechanism. Furthermore, it introduces a conditional random field to fuse local spectral similarity and behavioral consistency constraints, effectively addressing the problem of independent building type features and population behavior features lacking a collaborative decision-making mechanism, and eliminating the conflict between building space and behavioral dynamics during boundary generation. Attached Figure Description

[0053] Figure 1 This is a schematic diagram illustrating the workflow of the method for identifying migrant population settlements based on remote sensing imagery according to the present invention.

[0054] Figure 2 This is a schematic diagram of the method flow for identifying migrant population clusters based on remote sensing imagery according to the present invention. Detailed Implementation

[0055] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are some embodiments of the present invention, but not all embodiments.

[0056] Example 1, please refer to Figure 1 , Figure 2 This invention provides a technical solution: a method for identifying migrant population clusters based on remote sensing imagery, comprising the following steps:

[0057] Step 1: Integrate temporal Sentinel-2 spectral angle mapping with nighttime light autoencoder to separate the spectral drift characteristics of daytime temporary building materials and the spatial distribution of nighttime low-power lights in the target area, outputting the initial screening area for temporary housing, effectively separating temporary and permanent buildings, improving the accuracy of the initial screening, acquiring multi-temporal Sentinel-2L2A surface reflectance image sequences of the target area, constructing a spectral reference library for temporary building materials, calculating the temporal variance of spectral angles pixel by pixel, extracting the spectral drift feature map caused by daytime temporary building materials, effectively distinguishing the spectral differences between temporary building materials and permanent buildings, simultaneously acquiring nighttime low-power light images, using a convolutional autoencoder to separate the response of permanent buildings with strong light sources, extracting the distribution map of low-intensity clustered lights through decoder residual reconstruction, accurately separating low-power temporary lights from permanent buildings with strong light sources, performing pixel-level logical AND operation on the spectral drift feature map and the distribution map of low-intensity clustered lights, and outputting the initial screening area for temporary housing after morphological closing operation and area threshold filtering, integrating day and night features to remove static ground object interference, improving the purity of the initial screening area;

[0058] It should be noted that the Sentinel-2 L2A surface reflectance image sequence for the target area with no cloud or low cloud cover for each month of the year was obtained, totaling 12 periods. A temporary building material spectral reference library was constructed, including three typical materials: corrugated steel sheets, asbestos sheets, and plastic tarpaulins. At least 50 pure pixel spectra were collected for each type, and their mean values ​​were taken as the standard spectral vector. The spectral angle of each pixel was calculated with respect to the standard spectra of each type in the reference library, and the variance of the spectral angle over the 12 periods was statistically analyzed. A variance threshold of 0.15 radians² was set. Pixels exceeding this threshold and with a mean spectral angle between 0.3 and 0.8 radians were identified as daytime temporary building materials. The spectral drift characteristics induced are used to generate a binarized spectral drift feature map. A pixel value of 1 indicates the presence of temporary building material spectral anomalies, while 0 indicates no anomalies. Simultaneously, NPP-VIIRS monthly synthetic nighttime light images corresponding to Sentinel-2 are acquired, and the spatial resolution is resampled to 10 meters. A convolutional autoencoder is constructed, containing four 3×3 convolutional kernel layers with a stride of 2, outputting a response map of permanent buildings with strong light sources. The decoder reconstructs the original light through upsampling and residual concatenation, combining the input and reconstructed data. The difference is used as a residual map. Pixels with gray values ​​between 2 and 15 are extracted from the residual map, and the area of ​​continuous connected regions is not less than 300 square meters. This serves as a low-intensity clustered light distribution map. The spectral drift feature map and the light distribution map are subjected to pixel-level logical AND operation to obtain candidate pixels. A circular structuring element with a size of 5×5 pixels is used for closing operation, and connected regions with an area less than 600 square meters are removed to output the initial screening area of ​​the temporary housing area. After the closing operation and area threshold filtering are completed, a vector planar layer of the initial screening area of ​​the temporary housing area is generated. Each planar pixel records its minimum bounding rectangle and centroid coordinates. The initial screening area of ​​the temporary housing area is divided into regular grids with a grid size of 100 meters × 100 meters, forming an integer multiple relationship with the actual ground resolution of Sentinel-2 (10 meters). The proportion of temporary building material pixels and the density of low-power light pixels are statistically analyzed within each grid. At the same time, the raster mask of the initial screening area is output. The spatial reference system is consistent with the original remote sensing image to ensure spatial alignment with the nighttime light and trajectory data in subsequent steps.

[0059] Step 2: In the initial screening area of ​​the temporary residential area, GRU trajectory prediction and generative adversarial network are introduced to simulate the decoupling of the spatiotemporal dependence of work and residence in the travel chain. The stability of population flow and the stay pattern are judged by pattern classification. The work-residence dependence is decoupled and the population behavior pattern is finely divided. Individual travel time series are extracted in the initial screening area of ​​the temporary residential area. The trajectory of multiple consecutive days is used as input. The temporal dependence relationship is encoded by gated cyclic unit to predict the residence, workplace and stay point sequence in the travel chain of the next day. The temporal dependence of individual travel is effectively captured, and the prediction accuracy of residence and workplace is improved. By the error distribution between the predicted trajectory and the real trajectory, the spatiotemporal dependence strength of work and residence is decoupled. The spatiotemporal regularity of the round trip from residence to workplace and the frequency of daytime stay point migration are calculated to quantify the degree of work-residence binding and distinguish between regular commuting and random flow behavior. Pseudo-trajectories are generated by generative adversarial network to expand the sample. The GRU prediction error and work-residence dependence index are combined into a feature vector for subsequent classification of population flow stability and stay pattern, enriching the diversity of training samples and enhancing the generalization ability of flow pattern classification.

[0060] It should be noted that within the initial screening area of ​​the temporary accommodation zone, individual travel time series data covering 30 consecutive natural days were first extracted. Data sources included anonymous mobile phone signaling or public transportation card swipe records, with timestamp accuracy at the minute level. Each trajectory was sliced ​​on a daily basis to form a 24-hour discretized spatiotemporal sequence. The spatial grid adopted the same 100m × 100m subdivision as in step one. The trajectories of 7 consecutive days were used as the input sequence, and a two-layer gated recurrent unit network was adopted with the hidden layer dimension set to 128. The prediction model for the residence, workplace, and stopover sequence of the next day was trained through time backpropagation. The prediction output is the category probability distribution of the grid in each time period, where residence is defined as from 22:00 to 06:00 the next day. The grid with the longest dwell time during the period was defined as the workplace, which was defined as the grid with the highest frequency of occurrence between 09:00 and 17:00 during the daytime. Dwell points were defined as grids where the user stayed continuously for more than one hour during other time periods. During training, a joint loss function of mean squared error and classification cross-entropy was used, with a learning rate of 0.001, a batch size of 64, and a training epoch of no more than 200 epochs, with early termination. After trajectory prediction, the error distribution between the predicted trajectory and the actual trajectory was calculated for each individual. Error metrics included the residential grid offset distance (expressed as Euclidean distance), the workplace identification consistency rate, and the edit distance of the daytime dwell point sequence. Based on this, the spatiotemporal dependence of work and residence was decoupled, specifically the spatiotemporal regularity of the round trip between residence and workplace, using the morning and evening commuting periods (0... The reciprocal of the coefficient of variation (COP) for trajectory repetition (7:00-09:00 and 17:00-19:00) is used. A COP below 0.2 is considered strong regularity. The frequency of daytime stopover point migration is statistically analyzed daily, using the mean and standard deviation over 30 days. Simultaneously, whether the workplace in each travel chain is outside the initial screening area is recorded, and the migration rate is calculated as a supplementary indicator of the degree of job-residence separation. The error distribution and the spatiotemporal dependence of job and residence are standardized using Z-scores to form a 12-dimensional feature vector for each individual. A generative adversarial network (GAN) is used to generate pseudo-trajectories to expand the training sample set. The generator uses a GRU with 64 hidden units as a temporal encoder. The output trajectory sequence is evaluated for authenticity by a discriminator consisting of a three-layer convolutional network using Wafer Array (WA). The Sserstein loss function and gradient penalty term are used, with the gradient penalty coefficient set to 10. The generator and discriminator are trained alternately, with the generator updated twice and the discriminator updated once per round. The learning rate for both is 0.0001. The number of generated pseudo trajectories is controlled to be twice the number of real trajectories, and they must pass the spatial validity test, that is, the generated residential and workplace grids must be located within the initial screening area or its surrounding 500-meter buffer zone. The GRU prediction errors corresponding to the real and generated trajectories are combined with the job-housing dependence index to form an enhanced feature vector set, which serves as the input for the classification of population flow stability and dwelling pattern. The classifier uses a random forest with 100 trees and a maximum depth limit of 10. The output is the stability label and dwelling pattern label for each individual.

[0061] Furthermore, step 2 also includes: classifying individuals into work-residence bound type and work-residence free type based on the spatiotemporal dependence intensity of work and residence; calculating the work-residence bound ratio and daytime work relocation rate of each grid in the initial screening area of ​​temporary residence; quantifying the degree of dependence of grid population on external employment centers; dividing high-density short-term rental pattern and low-density nighttime stay pattern based on the kurtosis and skewness of the nighttime continuous stay duration distribution; labeling the dominant behavior category of each grid; distinguishing between frequent rental rotation and temporary overnight stay; quantifying mobility stability into the temporal variation coefficient of trajectory daytime similarity; and forming grid-level behavioral feature labels together with the stay pattern category to achieve joint representation of population mobility intensity and stay pattern.

[0062] It should be noted that the "work-residence binding" type is defined as individuals whose trajectory repetition variation coefficient during morning and evening commutes is less than 0.2 and whose workplace relocation rate is less than 0.3. The rest are classified as "work-residence free." Within a 100-meter grid, the work-residence binding ratio is calculated as the ratio of the number of bound individuals to the total number of individuals in the grid. The daytime workplace relocation rate is the average percentage of days when all individuals' workplaces are located outside the initial screening area of ​​temporary residences. These indicators directly reflect the degree of dependence between the population within the grid and external employment centers. Grids with a work-residence binding ratio higher than 0.6 tend to form stable clusters, while a relocation rate exceeding 0.5 indicates high mobility. Based on the kurtosis and skewness of the nighttime continuous stay duration distribution, stay patterns are classified. For the stay duration sequence of all individuals in each grid from 22:00 to 06:00 the next day, kurtosis and skewness are calculated. The high-density short-term rental pattern is characterized by a kurtosis greater than 3 and a skewness greater than 1, indicating that the stay duration is concentrated in a short period and has a relatively long right tail. The rental behavior should be frequently rotated; the low-density nighttime stay pattern is characterized by a kurtosis of less than 2 and an absolute skewness of less than 0.5, with an approximately uniform distribution, corresponding to the temporary overnight population. The dominant behavior category of each grid is determined by the distribution pattern of stay duration. If the sample sizes of the two patterns are comparable, they are labeled as mixed type. The mobility stability is quantified as the time series variation coefficient of trajectory daytime similarity. For all individuals in each grid, the daily cosine similarity of the trajectory grid sequence from 09:00 to 17:00 for 30 consecutive days is calculated. The ratio of the standard deviation of the mean to the mean is taken as the variation coefficient. A variation coefficient less than 0.25 is defined as high stability, between 0.25 and 0.60 is medium mobility, and greater than 0.60 is high mobility. This coefficient, together with the stay pattern category, constitutes the grid-level behavioral feature label. Among them, the high-density short-term rental pattern superimposed with a high mobility coefficient represents an extremely unstable settlement area, and the low-density nighttime stay pattern superimposed with a low mobility coefficient represents a relatively stable temporary overnight area.

[0063] Step 3: Construct a dual-agent deep Q-network, enabling the building feature agent and the behavior feature agent to make decisions alternately. Integrate local constraints from conditional random fields to generate initial settlement boundaries. Coordinate building and behavior features to eliminate misjudgments of single-feature boundaries. Construct a building feature agent with a state space consisting of spectral angle drift intensity, temporary building material probability, and local texture features for each pixel, and an action space consisting of binary classification of whether the pixel belongs to the settlement boundary. This enables autonomous classification of pixel-level building features. Construct a behavior feature agent with a state space consisting of grid-level flow stability indicators, dwelling pattern probability, and job-housing dependence intensity. Using the same action space, it forms a parallel decision-making structure with the building agent to achieve parallel discrimination of grid-level behavior features. Design an alternating decision-making strategy, in which the building agent prioritizes outputting the building-side classification results in each round, and the behavior agent corrects abnormal behavior grids based on these results. The two alternate iteratively until local convergence, and then output a classification probability map to eliminate boundary misjudgments driven by single features.

[0064] It should be noted that the building feature agent adopts a deep Q-network structure. Its state space consists of three types of features for each pixel: spectral angle drift intensity, temporary building material probability, and local texture features. The spectral angle drift intensity is derived from the 12-period spectral angle temporal variance calculated in step 1. The temporary building material probability is mapped from the spectral angle mean to the 0-1 interval using a Sigmoid function. The local texture features are extracted using a 3×3 window gray-level co-occurrence matrix to obtain contrast and homogeneity parameters. The action space is a binary classification, with a value of 1 indicating that the pixel belongs to the boundary of a residential area and a value of 0 indicating that it belongs to the outside. The network contains two fully connected hidden layers with 256 and 128 neurons respectively. The output layer uses a Softm... The ax function generates classification probabilities. During training, an ε-greedy strategy is employed, with ε initially set at 0.9 and decreasing by 0.01 every 1000 steps until it reaches 0.05. The reward function is set based on the overlap between the pixel and the initial screening area of ​​the temporary accommodation zone in step 1; overlap results in a +1 reward, while non-overlap results in a -0.5 penalty. The behavioral feature agent also uses a deep Q-network, with its state space based on a 100m × 100m grid-level feature space. Each grid contains three dimensions: mobility stability index, dwelling pattern probability, and job-housing dependence strength. The mobility stability index is the cosine similarity coefficient of variation of the daytime trajectory calculated in step 2; the dwelling pattern probability is the proportion of individuals in the high-density short-term rental pattern within the grid. For example, the job-housing dependency strength is represented by a combination of the job-housing binding ratio and the rate of relocation from the workplace. A value of 1 is assigned when the former is higher than 0.6 and the latter is lower than 0.3; otherwise, it is 0. The action space is consistent with the building agent, using binary classification. The network hidden layer structure consists of 128 and 64 neurons, and the output layer uses Softmax. The reward function of the behavioral agent is set based on the consistency of behavior within the grid; that is, a +1 reward is given when more than 70% of individuals within the grid belong to the same dwelling pattern, otherwise a -0.3 reward is given. In the alternating decision-making strategy, each round first involves the building feature agent classifying all pixels and outputting a binary classification map of the building side. Subsequently, the behavioral feature agent, using a 100-meter grid as a unit, classifies the data based on this classification map. Based on this, corrections are made. The correction rule is as follows: if the building-side classification classifies a grid as inside the boundary, but the dominant dwelling pattern determined by the behavioral feature agent inside the grid differs from the average pattern of the three adjacent grids by more than one category, then the behavioral agent flips the grid label to outside the boundary. After completing one round of correction, the corrected boundary map is re-input into the building agent. The building agent adjusts its Q-value estimation according to the new boundary distribution. The alternating iteration termination condition is set as follows: the number of boundary change pixels in two consecutive rounds is less than 0.5% of the total number of pixels, or the number of iterations reaches 20 rounds to form a classification probability map, ensuring that the building spatial features and behavioral dynamic features reach consistency in the local area and eliminating boundary misjudgment driven by a single feature.

[0065] Furthermore, step 3 also includes: using the classification probability maps output by the alternating decisions of the two agents as the unary potential energy of the conditional random field, where the building potential energy comes from the output of the building agent and the behavior potential energy comes from the correction result of the behavior agent, effectively integrating the two heterogeneous features of building and behavior, constructing the binary potential energy of the conditional random field, penalizing the label inconsistency between adjacent pixels or grids, and introducing spectral similarity and behavior consistency constraints in the local neighborhood to suppress local noise and enhance the continuity of the boundary space. The conditional random field is solved using the mean field inference algorithm, and under the premise of maintaining the local consistency of building and behavior, a smooth initial settlement boundary that meets the local constraints is output, generating a smooth boundary that takes into account both spectral and behavior consistency.

[0066] It should be noted that architectural potential directly uses the boundary probability values ​​output by the architectural feature agent for each pixel, while behavioral potential originates from the probability distribution of the behavioral feature agent, which is then interpolated to the pixel after correction using a 100-meter grid. Univariate potential is calculated using the pixel as the basic unit. Architectural and behavioral potentials are stored independently and maintain a consistent spatial resolution of 10 meters. During the initialization phase of the conditional random field, the probability map files output by the architectural and behavioral agents in their final rounds are loaded, with each pixel corresponding to two probability values, serving as prior information for subsequent binary potential optimization. The construction of binary potential focuses on penalizing label inconsistencies between adjacent pixels, while introducing spectral similarity and behavioral consistency constraints within the local neighborhood. The neighborhood system is defined as a four-neighbor or eight-neighbor system. The feature vectors of adjacent pixel pairs include spectral angle difference values ​​and dwell mode category difference values. The spectral angle difference is obtained by calculating the spectral angular distance between two pixels on the Sentinel-2 image. Penalty weights are reduced when the difference is less than 0.1 radians and reduced when the difference is greater than 0.3 radians. The penalty weight is increased, and the behavioral consistency constraint compares whether the dominant dwelling patterns of the grids where two pixels are located are the same. If the patterns are different, the binary potential energy is increased by a fixed penalty value of 0.5. The spatial distance weight of adjacent pixels is set according to the reciprocal of the Euclidean distance. The closer the distance, the higher the penalty weight. All binary potential energy parameters are determined on the validation set through grid search and do not depend on subjective settings. The mean field inference approximates the label distribution of each pixel as the product of independent distributions. The label probability of each pixel is updated iteratively until convergence. In each iteration, the unary potential energy of each pixel remains unchanged, and the binary potential energy is calculated based on the expected label distribution of adjacent pixels in the current round. The convergence condition is set as the sum of the absolute values ​​of the changes in the labels of all pixels between two consecutive iterations is less than one percent of the total number of pixels, or the maximum number of iterations of 50 is reached. The label probability of each pixel output by the mean field inference is binarized after a threshold of 0.5 and then connected component analysis is performed to remove isolated areas with an area of ​​less than 300 square meters. Finally, a smooth initial settlement boundary that meets the local consistency constraints of architecture and behavior is output.

[0067] Step 4: The two agents iteratively adjust their weights based on the feedback from the initial settlement boundary. They eliminate boundary conflicts between building space and behavioral dynamics through a joint reward mechanism. By iteratively eliminating boundary conflicts, they achieve collaborative judgment of the two features. The joint reward function is defined as the harmonic average of the entropy of building features and the entropy of behavioral features within the initial settlement boundary. A high reward is given when both are highly consistent, and a low reward is given when there is a conflict between building and behavior. This effectively quantifies the degree of collaboration between building and behavior. The two agents update their respective deep Q-network weights based on the joint reward. They record the difference in contribution of the two agents to boundary generation in each iteration. The state fusion weight of the agent with excessive contribution is reduced to dynamically balance the decision influence of the two agents. The alternating decision-making and conditional random field boundary generation are repeated until the joint reward no longer increases for several consecutive rounds. The settlement boundary after collaborative judgment of building space and behavioral dynamics is output, ensuring that the boundary converges to the optimal collaborative state.

[0068] It should be noted that the joint reward function is calculated using all pixels within the initial settlement boundary as the statistical object. First, the architectural feature entropy and behavioral feature entropy are calculated separately: the architectural feature entropy is calculated based on the spectral angle drift intensity histogram of pixels within the boundary, divided into 32 equally spaced intervals, and then the Shannon entropy is calculated; the behavioral feature entropy is calculated based on the dominant dwelling pattern category distribution of each 100-meter grid within the boundary. Then, the harmonic average of the two is taken as the joint reward value. When the architectural feature entropy is below 1.5 and the behavioral feature entropy is below 0.8, it indicates high consistency between the two, and a reward of +1.0 is given; when the architectural feature entropy is above 2.5 or the behavioral feature entropy is above 1.2, it is judged as architectural-behavioral conflict, and a reward of -0.8 is given; otherwise, a reward of 0.2 is given.

[0069] The expression for calculating the entropy of building features is as follows:

[0070] ;

[0071] In the formula: The building feature entropy reflects the degree of disorder in the distribution of pixels in terms of spectral angle drift intensity within the boundary of a settlement. The higher the entropy value, the more diverse and inconsistent the building spectral features are. This represents the total number of intervals into which the histogram is divided, and its value is 32. This is the interval index of the histogram, with values ​​ranging from 1 to... ; For all pixels within the boundary, the spectral angular shift intensity value falls within the first... The proportion of pixels in each interval;

[0072] The expression for calculating behavioral feature entropy is as follows:

[0073] ;

[0074] In the formula: The behavioral feature entropy reflects the degree of disorder in the distribution of dominant dwelling mode categories in each 100-meter grid within the boundary of the settlement. The higher the entropy value, the more severe the mixing of grids with different behavioral modes. The total number of categories representing the dominant dwell time pattern, with a value of 3; This is the pattern category index, with values ​​ranging from 1 to... ; For all 100-meter grids within the boundary, the dominant dwell mode belongs to the [missing information]. The proportion of the grid space occupied by the class;

[0075] The harmonic average of architectural feature entropy and behavioral feature entropy is calculated as follows:

[0076] ;

[0077] The entropy range is mapped to the final joint reward value. Its expression is as follows:

[0078] ;

[0079] In the formula: This is the harmonic mean of architectural feature entropy and behavioral feature entropy, used to measure the degree of consistency between the two. The harmonic mean is sensitive to extreme values; it only becomes effective when both entropy values ​​are small. It will be noticeably smaller; The final joint reward value after threshold mapping is directly used to guide the weight updates and iterative convergence determination of the two agents. This indicates a high degree of synergy between architectural spatial characteristics and human behavioral characteristics. This indicates a significant conflict between the two. represents moderate consistency; represents the high consistency and conflict determination threshold of architectural feature entropy, determined by experimental statistics; represents the high consistency and conflict determination threshold of behavioral feature entropy, determined by experimental statistics.

[0080] An experience replay mechanism is used to store the state, action, reward, and next state of each round of alternating decision-making. After each round of alternating decision-making and conditional random field boundary generation is completed, the contribution of the building agent and the behavior agent to the final boundary change is recorded. The contribution is represented by the average absolute difference between the output probability map of the agent in this round and the output probability map of the previous round. For the side with too high contribution (above 0.15), its feature weight in the next round of state fusion is reduced by 20%, while for the side with too low contribution (below 0.05), it is increased by 15%. After the weight adjustment, the two agents re-execute the alternating decision-making and conditional random field inference in step 3 to generate a new round of initial boundary.

[0081] The expression for calculating the contribution of a building agent is as follows:

[0082] ;

[0083] In the formula: For the first The contribution of the building agent to the boundary changes in each iteration; the larger the value, the more significant the change in the agent's output in this round compared to the previous round. This is the iteration round index, numbered starting from 1; The total number of pixels within the study area; For cell indexes, the value range is from 1 to ; For the first In the round of iteration, the output probability graph of the building agent is the th The probability value inside the boundary of each cell; For the first In the round of iteration, the output probability graph of the building agent is the th The probability value inside the boundary of each cell;

[0084] The expression for calculating the contribution of the behavioral agent is as follows:

[0085] ;

[0086] In the formula: For the first The contribution of the agent to the boundary changes in each iteration; the larger the value, the more significant the change in the agent's output in this round compared to the previous round. This represents the total number of 100-meter grid cells within the study area. For grid indices, values ​​range from 1 to... ; For the first In the round of iteration, the output probability graph of the behavioral agent is the th The probability value inside the boundary of each grid; For the first In the round of iteration, the output probability graph of the behavioral agent is the th The probability value inside the boundary of each grid;

[0087] Repeat the above process of alternating decision-making, conditional random field boundary generation, and joint reward calculation until the absolute value of the change in the joint reward value is less than 0.05 for 5 consecutive rounds. It is determined that it will no longer increase. At this time, the boundary map optimized by the conditional random field in the last round is output as the boundary of the settlement after collaborative judgment. This boundary has eliminated the conflict between the spatial distribution of buildings and the dynamic pattern of behavior through the iterative adjustment of the joint reward mechanism, so that the boundary simultaneously satisfies the two conditions of abnormal clustering of temporary building materials spectrum and consistency of population behavior pattern. The output results include a pixel-level binary mask (10-meter resolution) and a vector area layer. Each settlement area cell is accompanied by the final value of the joint reward and the final contribution record of the two agents.

[0088] Step 5: Utilize the difference between the pseudo-trajectory and real trajectory distribution output by the generative adversarial network to quantify the stability of the mobile population, correct the division of high-density and low-density areas within the boundary and form a stability-level patch boundary, quantify the flow stability, and achieve refined classification of settlement areas.

[0089] Step 6: Globally optimize boundary smoothness based on conditional random field, integrate the final decision of the two agents, output the spatial range and stability classification results of the settlement area, globally optimize boundary smoothness, and output the connectivity classification map.

[0090] Example 2, as Figure 1 , Figure 2 As shown, based on Embodiment 1, the present invention provides a technical solution: Step 5 specifically includes: within the boundary of the settlement area after collaborative determination, extracting the real travel trajectory distribution of each grid, and simultaneously running the trained generative adversarial network to generate a corresponding number of pseudo trajectory distributions, realizing the benchmark alignment between real behavior and generated patterns, calculating the Wasserstein distance between the real trajectory distribution and the pseudo trajectory distribution in each grid as a quantitative indicator of trajectory pattern complexity. The larger the distance, the more difficult it is to predict the actual behavior, quantifying the complexity and randomness of behavior. After normalizing the Wasserstein distance, it is mapped to the stability index of the mobile population. The higher the index, the worse the stability, and the lower the index, the more regular and stable the settlement behavior, thus constructing a unified measurement scale for population stability.

[0091] It should be noted that, using the 100-meter grid output from step 2 as the basic unit, the actual travel trajectory sequence of all individuals within each grid over 30 consecutive natural days is extracted. The trajectory is sliced ​​on a daily basis, and the spatial location is mapped to the corresponding grid. Simultaneously, the generative adversarial network trained in step 2 is run. The generator uses a GRU structure with 64 hidden units to generate a pseudo-trajectory sequence for each grid with an equal number of actual trajectories. The generation process is validated for spatial validity, meaning that the generated residential and workplace grids must be located within the initial screening area of ​​temporary residences or their surrounding 500-meter buffer zone. Both actual and pseudo-trajectories are stored as grid sequences. The starting point of each trajectory is the residential grid at night, and the daily trajectory length is a 24-hour discretized time period. For each grid, the distribution of actual and pseudo-trajectories is statistically analyzed. The distribution is defined as the probability matrix of the occurrence of all trajectories within that grid through each time period and each grid, with the dimension being 24 time periods multiplied by the total number of grids. The Wasserstein distance between the two probability distributions is calculated using a one-dimensional Earth Mover's distance. The Distance approximation algorithm uses the Euclidean distance between grids as a measure of ground distance. The larger the Wasserstein distance value, the further the actual trajectory deviates from the regular trajectory pattern learned by the generative adversarial network. In other words, the actual travel behavior of people is more complex and more difficult to predict with a fixed pattern. Conversely, it indicates that the actual behavior is highly consistent with the pattern captured by the generative model.

[0092] The expression for calculating the Wasserstein distance is as follows:

[0093] ;

[0094] Its discrete approximation calculation formula is:

[0095] ;

[0096] In the formula: For grid The Wasserstein distance between the true trajectory distribution and the pseudo trajectory distribution; For grid indices, values ​​range from 1 to... ; This is the infimum operator, which represents finding the minimum transportation cost among all possible joint distributions; It is a joint probability distribution, belonging to the set ; This indicates that the unit mass is from position Transport to location The probability density; For all and It is the set of joint probability distributions of marginal distributions; The metric space is formed by all possible trajectory points; and For metric space Two trajectory points in the middle; The ground distance function is defined as the distance per unit mass from the trajectory point. Transport to track point The cost is achieved by using the Euclidean distance between the centroids of the grids containing the two trajectory points; For joint distribution At point The differential mass element at point represents the mass element from... Transport to The mass of the infinitesimal element; In the product space On transportation costs According to joint distribution Integrate the components to obtain the total cost of the transportation plan; The length of the one-dimensional vector after flattening the spatiotemporal distribution matrix; The index of the flattened one-dimensional vector, with values ​​ranging from 1 to... ; For grid The flattened vector of the true trajectory distribution, and the cumulative distribution function after sorting by ground distance in ascending order, are on the th... The value of each position; For grid The flattened vector of the pseudo-trajectory distribution, and the cumulative distribution function after being arranged according to the same sorting rule, are on the th... The value of each position; For the two cumulative distribution functions at the th The absolute difference at each position, for all Summing and dividing by This yields an approximation of the one-dimensional EarthMover's Distance, which is in... It has high computational efficiency when using Euclidean distance;

[0097] The Wasserstein distance calculated for each grid is min-max normalized over the entire region. The normalized value ranges from 0 to 1 and is defined as the stability index of the mobile population. The closer the index is to 1, the worse the stability, corresponding to populations with complex trajectory patterns and high diurnal variability. The closer the index is to 0, the better the stability, corresponding to stable populations with regular clustering behavior and strong spatiotemporal dependence on work and residence. This index is cross-validated with the diurnal similarity variation coefficient of the trajectory calculated in step 2. A correlation coefficient higher than 0.7 is considered valid. Finally, the stability index is divided into three levels: high stability, medium mobility, and high mobility, with thresholds of 0.33 and 0.67. The stability index value and level label of each grid are output.

[0098] The formula for calculating the stability index of the mobile population is as follows:

[0099] ;

[0100] In the formula: For grid The stability index of the mobile population; The minimum Wasserstein distance among all grids in the entire region is used as the lower bound for normalization. This minimum value corresponds to the grid in the entire region where the distribution of the true trajectory and the pseudo trajectory are closest, i.e., the grid with the most regular behavior. The maximum value of Wasserstein distance in all grids of the entire region is used as the upper bound of normalization, corresponding to the grid with the greatest difference between the distribution of real and pseudo trajectories in the entire region, that is, the grid with the most chaotic behavior and the most unpredictable behavior.

[0101] Furthermore, step 5 also includes: comparing the stability index of the mobile population in each grid of the original high-density short-term rental area with a preset high stability threshold; grids below the threshold are relabeled as stable settlements, while those above the threshold remain high-density mobile areas; grids with regular behavior in the high-density short-term rental area are identified as stable settlements; comparing the stability index of the mobile population in each grid of the original low-density nighttime stay area with a preset low stability threshold; grids above the threshold are relabeled as extremely mobile temporary areas, while those below the threshold remain low-density stable areas; grids with random behavior in the low-density nighttime stay area are classified as extremely mobile temporary areas; and outputting the corrected four-category area division results to form the stability classification map boundary, covering high-density stable areas, high-density mobile areas, low-density stable areas, and low-density mobile areas, outputting four types of stability classification maps to refine the behavioral characteristics of the settlements;

[0102] It should be noted that after calculating the stability index of the mobile population, the stability of the grids marked as high-density short-term rental patterns in step 2 is assessed. The stability index value of each such grid is extracted and compared with a preset high stability threshold of 0.33. Grids with a mobile population stability index below 0.33 indicate that the actual travel behavior of the people within them is highly regular, with a daytime trajectory similarity coefficient of variation of less than 0.25 and a work-residence binding ratio of more than 0.6. This does not meet the frequent rotation characteristics that high-density short-term rentals should have, and these grids are re-marked as high-density stable areas. Grids with a mobile population stability index above 0.33 have higher trajectory pattern complexity and are more prone to over-the-counter travel. The normalized in-distance value reflects the degree of behavioral variability consistent with the expected high-density short-term rental pattern, maintaining the original labeling of its high-density flow area unchanged. This threshold is determined based on the tertiary position of the overall stability index distribution, ensuring statistically significant differences in behavioral characteristics between the two types of areas. For the grids labeled as low-density nighttime stay patterns in step 2, stability is assessed, and the flow population stability index of each such grid is extracted and compared with the preset low stability threshold of 0.67. Grids with a flow population stability index higher than 0.67 have the Wasserstein distance between their true trajectory distribution and the pseudo trajectory distribution generated by the generative adversarial network at the highest tertiary position in the entire area, indicating that the flow population... The complex travel behavior and highly random daytime activity patterns lack the relatively fixed overnight characteristics implied by low-density nighttime stays. Therefore, these grids are relabeled as highly mobile temporary zones. Grids with a mobile population stability index below 0.67 exhibit relatively regular trajectory patterns, with nighttime stay duration distribution kurtosis less than 2 and absolute skewness less than 0.5, consistent with the uniform distribution characteristics of low-density nighttime stays. Their original low-density stable zone labeling is maintained. This threshold is also set based on the ternary loci of the overall mobile population stability index, forming a symmetrical division with the high stability threshold. The corrected four-zone division results are spatially integrated to output the stability grading patch boundaries. The four zones include: high-density stable zones (… The original high-density short-term rental model included grids with a stability index of less than 0.33 for the mobile population, high-density mobile areas (grids with a stability index of more than 0.33 for the mobile population in the original high-density short-term rental model), low-density stable areas (grids with a stability index of less than 0.67 for the mobile population in the original low-density nighttime stay model), and low-density mobile areas, i.e., extremely mobile temporary areas (grids with a stability index of more than 0.67 for the mobile population in the original low-density nighttime stay model). The boundaries of each patch were aggregated with a 100-meter grid as the smallest unit. Isolated patches with an area of ​​less than 900 square meters (i.e., 9 consecutive grids) were removed. The output was a vector area layer with the average stability index of each grid and the partition category code.

[0103] Step 6 specifically includes: using the boundary of the settlement after collaborative judgment of the final decision output of the two agents as the unary potential energy of the conditional random field, using the boundary of the four types of stability classification patches as the regional connectivity constraint of the binary potential energy, ensuring that the boundary integrates the dual prior information of architecture and behavior, introducing higher-order potential energy to penalize isolated regions with an area smaller than the preset minimum connected patch, forcibly merging them into the adjacent dominant category, using the graph cut algorithm to globally optimize the boundary smoothness and spatial connectivity, eliminating noisy patches, improving the continuity and integrity of the boundary space, outputting the vector boundary of the settlement space range after global optimization, and attaching a stability classification distribution map inside each settlement, including the area proportion of each classification and the average flow stability index, providing refined attribute data that can be directly used for spatial analysis and mapping;

[0104] It should be noted that the settlement boundary after collaborative determination output in step 4 is used as the unary potential energy of the conditional random field. This boundary is a pixel-level binary mask, with each 10-meter pixel corresponding to the probability value inside the boundary after collaborative determination by the building and behavior dual agents. Simultaneously, the boundaries of the four types of stability-graded patches output in step 5 are used as the regional connectivity constraints of the binary potential energy. The four types of regions include high-density stable regions, high-density flow regions, low-density stable regions, and low-density flow regions. Each type of patch uses a 100-meter grid as the smallest unit. Both the unary and binary potential energies are weighted... Sampling was performed to a spatial resolution of 10 meters to ensure consistency with the original remote sensing image spatial reference frame. To address the potential for small, isolated patches during conditional random field optimization, a higher-order potential energy was introduced for penalty processing. A minimum connected patch area threshold of 900 square meters was set, corresponding to nine consecutive 100-meter grids or 8100 10-meter pixels. Isolated patches with areas smaller than this threshold were forcibly merged into the adjacent region containing the patch with the largest dominant category area. The merging rule was based on the category priority of adjacent patches: high-density stable regions had the highest priority, followed by... The algorithm identifies high-density flow areas and low-density stable areas, with low-density flow areas having the lowest priority. A graph cut algorithm is used to globally optimize the conditional random field. The iteration termination condition is set to the number of pixels with changed labels being less than 0.5% of the total number of pixels in two consecutive iterations, or reaching the maximum number of iterations (50). The output is a settlement boundary with smooth boundaries and optimized spatial connectivity. Finally, the globally optimized spatial extent vector boundary of the settlement is output in Shapefile format. Each settlement cell includes internal stability hierarchy distribution attributes, with attribute fields including: the area proportion of high-density stable areas, the area proportion of high-density flow areas, the area proportion of low-density stable areas, the area proportion of low-density flow areas, and the arithmetic mean of the stability index of the mobile population within each zone. The overall average flow stability index of each settlement is also output, calculated by weighting the stability indices of all grids within the zone by area. The smallest mapping unit of the vector boundary is a continuous patch after aggregating 100-meter grids. The spatial reference system remains consistent with the input remote sensing image, allowing direct use for subsequent spatial statistical analysis and mapping.

[0105] The workflow of this method for identifying migrant population clusters based on remote sensing imagery will be explained in detail below.

[0106] First, monthly Sentinel-2 L2A image sequences and corresponding NPP-VIIRS nighttime light images of the target area were acquired to construct a temporary building material spectral reference library. The temporal variance of the spectral angles for 12 periods was calculated pixel-by-pixel to extract daytime spectral drift features. Simultaneously, a convolutional autoencoder was used to reconstruct residuals and extract low-intensity clustered light distributions. After logical AND operations, morphological closing operations, and area threshold filtering, the initial screening area of ​​the temporary residential area was output and divided into a 100-meter grid. Second, the travel trajectories of individuals over 30 consecutive days were extracted from the initial screening area. A two-layer GRU network was used to predict the next day's residence, workplace, and stopover sequence. The spatiotemporal dependence strength was decoupled by the prediction error and the work-residence dependence index. The samples were expanded using pseudo-trajectories generated by GAN, and individual stability and stay pattern labels were output using random forest, aggregated into grid-level behavioral features. Finally, a dual-intelligence architecture combining building features and behavioral features was constructed. A volumetric deep Q-network is used to generate an initial boundary by fusing local constraints from a conditional random field after alternating decision-making. The weights of the two agents are iteratively adjusted through a joint reward mechanism until convergence, outputting the settlement boundary after collaborative judgment. Next, the Wasserstein distance between the real trajectory and the pseudo trajectory distribution of GAN for each grid is calculated and normalized into a stability index of the mobile population. Based on this, the high-density short-term rental pattern is corrected into a high-density stable area or a high-density mobile area, and the low-density nighttime stay pattern is corrected into a low-density stable area or a low-density mobile area, forming four types of stability-level patches. Finally, the collaborative judgment boundary is used as a univariate potential energy and the four types of leveled patches are used as a binary potential energy. A higher-order potential energy is introduced to penalize small-area isolated patches, and a graph cut algorithm is used to globally optimize the boundary smoothness and spatial connectivity, outputting a settlement vector boundary with stability-level distribution attributes.

[0107] The above are merely specific embodiments of the present invention, but the scope of protection of the present invention is not limited thereto. The scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A method for identifying migrant population clusters based on remote sensing imagery, characterized in that, Includes the following steps: Step 1: Integrate temporal Sentinel-2 spectral angle mapping with nighttime light autoencoder to separate the spectral drift characteristics of daytime temporary building materials and the spatial distribution of nighttime low-power lights in the target area, and output the initial screening area of ​​temporary residential area. Step 2: Introduce GRU trajectory prediction and generative adversarial network in the initial screening area of ​​temporary accommodation to simulate the decoupling of work-residence spatiotemporal dependence of travel chain and determine the stability of population flow and stay pattern through pattern classification. Step 3: Construct a dual-agent deep Q-network, enabling the building feature agent and the behavioral feature agent to make decisions alternately, and integrate local constraints of conditional random fields to generate the initial settlement boundary; Step 4: The two agents iteratively adjust their weights based on the feedback from the initial settlement boundary, and eliminate the boundary conflict between building space and behavioral dynamics through a joint reward mechanism. Step 5: Utilize the difference between the pseudo-trajectory and real trajectory distribution output by the generative adversarial network to quantify the stability of the mobile population, correct the division of high-density and low-density regions within the boundary, and form a stability-level patch boundary. Step 6: Globally optimize boundary smoothness based on conditional random field, integrate the final decision of the two agents, and output the spatial range and stability classification results of the settlement area.

2. The method for identifying migrant population clusters based on remote sensing imagery according to claim 1, characterized in that: Step 1 specifically includes: Acquire multi-temporal Sentinel-2L2A surface reflectance image sequences of the target area, construct a spectral reference library for temporary building materials, calculate the temporal variance of spectral angles pixel by pixel, and extract spectral drift feature maps induced by daytime temporary building materials. Simultaneously acquire low-power nighttime light images, use a convolutional autoencoder to separate the response of strong light sources in permanent buildings, and extract low-intensity clustered light distribution maps through decoder residual reconstruction. The spectral drift feature map and the low-intensity clustered light distribution map are subjected to pixel-level logical AND operation, followed by morphological closing operation and area threshold filtering to output the initial screening area of ​​temporary residential area.

3. The method for identifying migrant population clusters based on remote sensing imagery according to claim 1, characterized in that: Step 2 specifically includes: Individual travel time series are extracted from the initial screening area of ​​temporary residential areas. Using trajectories from multiple consecutive days as input, gated cyclic unit encoding of temporal dependencies is used to predict the sequence of residence, workplace and stopover points in the travel chain for the next day. By analyzing the error distribution between the predicted trajectory and the actual trajectory, the spatiotemporal dependence of work and residence is decoupled, and the spatiotemporal regularity of the round trip from residence to workplace and the frequency of daytime stop point migration are calculated. Generative adversarial networks are used to generate pseudo-trajectories to expand the samples, and the GRU prediction error is combined with the job-housing dependence index into a feature vector.

4. The method for identifying migrant population clusters based on remote sensing imagery according to claim 3, characterized in that: Step 2 also includes: Individuals are classified into work-residence bound type and work-residence free type based on the intensity of work-residence spatiotemporal dependence. The work-residence binding ratio and daytime workplace relocation rate of each grid in the initial screening area of ​​temporary housing are calculated. Based on the kurtosis and skewness of the distribution of continuous nighttime stay duration, high-density short-term rental patterns and low-density nighttime stay patterns are divided, and the dominant behavioral categories of each grid are labeled. The flow stability is quantified as the temporal variation coefficient of trajectory diurnal similarity, which together with the dwell pattern category constitutes a grid-level behavioral feature label.

5. The method for identifying migrant population clusters based on remote sensing imagery according to claim 1, characterized in that: Step 3 specifically includes: Construct a building feature intelligent agent, with the state space consisting of the spectral angle drift intensity, temporary building material probability, and local texture features of each pixel, and the action space consisting of binary classification of whether it belongs to the boundary of a settlement. Construct a behavioral characteristic intelligent agent with a state space consisting of grid-level flow stability indicators, dwelling pattern probabilities, and job-housing dependence intensity. Adopt the same action space and form a parallel decision-making structure with the building intelligent agent. The design employs an alternating decision-making strategy, in which the building agent prioritizes outputting the building-side classification result in each round, and the behavioral agent corrects the behavioral anomaly grid based on this. The two alternately iterate until local convergence, and then output a classification probability map.

6. The method for identifying migrant population clusters based on remote sensing imagery according to claim 5, characterized in that: Step 3 also includes: The classification probability maps of the alternating decision-making outputs of the two agents are used as the unary potential energy of the conditional random field, where the building potential energy comes from the output of the building agent and the behavior potential energy comes from the correction results of the behavior agent. We construct a binary potential energy for a conditional random field to penalize label inconsistencies between adjacent pixels or grids, while introducing spectral similarity and behavioral consistency constraints within the local neighborhood. The mean-field inference algorithm is used to solve the conditional random field. While maintaining the local consistency between buildings and behaviors, it outputs a smooth initial settlement boundary that meets local constraints.

7. The method for identifying migrant population clusters based on remote sensing imagery according to claim 1, characterized in that: Step 4 specifically includes: The joint reward function is defined as the harmonic average of the entropy of building features and the entropy of behavior features within the initial settlement boundary. A high reward is given when both are highly consistent, and a low reward is given when there is a conflict between building and behavior. The two agents update their respective deep Q-network weights based on joint rewards, record the difference in their contributions to boundary generation in each iteration, and reduce the state fusion weight of the agent with excessive contribution. Repeated alternating decision-making and conditional random field boundary generation until the joint reward no longer increases after multiple consecutive rounds, outputting the settlement boundary after dynamic collaborative determination of building space and behavior.

8. The method for identifying migrant population clusters based on remote sensing imagery according to claim 7, characterized in that: Step 5 specifically includes: Within the boundaries of the settlement area after collaborative determination, the actual travel trajectory distribution of each grid is extracted, and the trained generative adversarial network is run to generate a corresponding number of pseudo trajectory distributions. Calculate the Wasserstein distance between the true trajectory distribution and the pseudo trajectory distribution within each grid, and use it as a quantitative indicator of trajectory pattern complexity; By normalizing the Wasserstein distance and mapping it to the stability index of the mobile population, a unified metric for population stability is constructed.

9. The method for identifying migrant population clusters based on remote sensing imagery according to claim 8, characterized in that: Step 5 further includes: The stability index of the mobile population in each grid within the original high-density short-term rental area is compared with a preset high stability threshold. Grids below the threshold are re-marked as stable residential areas, while those above the threshold are maintained as high-density mobile areas. The stability index of the mobile population in each grid within the original low-density nighttime residence area is compared with a preset low stability threshold. Grids with values ​​higher than the threshold are re-marked as extremely mobile temporary zones, while those with values ​​lower than the threshold remain low-density stable zones. The corrected four-category region division results are output to form the stability classification map boundary, covering high-density stable region, high-density flow region, low-density stable region and low-density flow region.

10. The method for identifying migrant population clusters based on remote sensing imagery according to claim 1, characterized in that: Step 6 specifically includes: The boundary of the settlement after the collaborative judgment of the final decision output of the two agents is used as the unary potential energy of the conditional random field, and the boundary of the output four types of stability hierarchical patches is used as the regional connectivity constraint of the binary potential energy. A higher-order potential energy penalty is introduced for isolated regions whose area is smaller than the preset minimum connected patch, and these regions are forcibly merged into the adjacent dominant category. A graph cut algorithm is used to globally optimize the boundary smoothness and spatial connectivity. Output the vector boundary of the settlement space after global optimization, and attach a stability classification distribution map of each settlement, including the area proportion of each classification and the average flow stability index.

Images