A method and system for surveying and mapping national land space planning based on dynamic remote sensing technology

Abstract

The present application provides a method and system for surveying and mapping of national land space planning based on dynamic remote sensing technology. The method includes: constructing a multi-source remote sensing data acquisition network based on the LPV model; performing spatiotemporal registration and noise filtering on time series remote sensing data; extracting spatiotemporal features using a deep neural network; identifying key change areas through a dual attention module; constructing a hierarchical dynamic feature library; establishing an uncertainty distribution model; designing probability constraints; generating adaptive safety spacing rules; establishing a multi-objective optimization function; and implementing a spatial planning strategy using the NSGA‑III algorithm. The present application improves the intelligence level and reliability of national land space planning through dynamic remote sensing data analysis and distributed robust opportunity constraint strategies.

Description

Technical Field

[0001] The present application relates to the fields of remote sensing technology, artificial intelligence, and spatial planning, and specifically to an intelligent surveying and mapping method and system for national land space planning using dynamic remote sensing data. Background Art

[0002] National land space planning is an important means to guide the protection, development, utilization and restoration of national land space, and is of great significance for achieving sustainable development of national land space. The current national land space planning surveying and mapping mainly relies on static remote sensing data and traditional planning methods. Existing technologies usually use remote sensing images at fixed time points for data collection and analysis, such as obtaining ground feature information from single-phase satellite images, and make planning decisions based on these static data.

[0003] In practical applications, a common technical solution is to use optical remote sensing satellites to obtain ground object images, combine manual interpretation and simple image classification algorithms to identify ground objects, and then use traditional deterministic constraint methods for planning and layout. For example, a collision avoidance strategy with a fixed threshold is used to determine the safe distance between different functional areas, or a preset land use zoning rule is used for spatial layout.

[0004] However, the existing technology has the following major problems: First, the analysis method based on static remote sensing data cannot effectively capture the dynamic characteristics of ground object changes, resulting in a lack of support for time-series evolution in planning decisions and difficulty in adapting to complex and changing actual environments. Second, the traditional deterministic constraint collision avoidance strategy has poor robustness and is prone to failure in actual environments with various uncertainties, affecting the reliability and security of the planning scheme. Summary of the invention

[0005] In view of this, the present application provides a land space planning surveying and mapping method and system based on dynamic remote sensing technology, which solves the problem in the prior art that it is difficult to capture the dynamic change characteristics of ground objects based on static remote sensing data and the poor robustness of deterministic constraints.

[0006] The present application embodiment provides a land space planning surveying and mapping method based on dynamic remote sensing technology, including:

[0007] Constructing a multi-source remote sensing data acquisition network based on a linear parameter variation (LPV) model, wherein the LPV model includes a state matrix describing the changes of ground objects and a time-varying parameter vector describing seasonal changes and climate conditions;

[0008] Acquire time-series remote sensing data including satellite remote sensing, UAV remote sensing and ground observation station remote sensing through the multi-source remote sensing data acquisition network;

[0009] Performing spatiotemporal registration and noise filtering on the time series remote sensing data according to the state matrix and the time-varying parameter vector to obtain standardized dynamic remote sensing data, wherein the standardized dynamic remote sensing data includes multispectral data, SAR data and DEM data;

[0010] According to the standardized dynamic remote sensing data, a deep neural network including a spatiotemporal convolution layer and an attention mechanism is used to extract spatiotemporal features;

[0011] The key change regions are identified through a dual attention module based on the extracted spatiotemporal features;

[0012] The spatiotemporal characteristics and key change area information are stored in a hierarchical dynamic feature library to obtain dynamic change characteristics of the ground object;

[0013] According to the dynamic change characteristics of the ground objects, the Gaussian mixture model and the Beta distribution model are used to construct an uncertainty distribution model;

[0014] designing probability constraints according to the uncertainty distribution model;

[0015] Based on the uncertainty distribution model of the probability constraint, the minimum safety distance is calculated in combination with the uncertainty ellipsoid, an adaptive safety distance rule is generated, and a collision avoidance rule is obtained through verification;

[0016] According to the collision avoidance rule, a multi-objective optimization function of land use efficiency index and environmental protection index is established;

[0017] The collision avoidance rules are converted into mathematical inequalities as constraints, and the NSGA-III algorithm is used to solve the multi-objective optimization problem to generate multiple non-dominated solutions;

[0018] Execute a space planning strategy based on multiple nondominated solutions.

[0019] The present application also provides a land space planning and mapping system based on dynamic remote sensing technology, including:

[0020] An acquisition module is used to construct a multi-source remote sensing data acquisition network based on a linear parameter variation LPV model, wherein the LPV model includes a state matrix describing changes in ground objects and a time-varying parameter vector describing seasonal changes and climate conditions. Through the multi-source remote sensing data acquisition network, time-series remote sensing data including satellite remote sensing, UAV remote sensing and ground observation station remote sensing are acquired;

[0021] An extraction module, for performing spatiotemporal registration and noise filtering on the time series remote sensing data according to the state matrix and the time-varying parameter vector to obtain standardized dynamic remote sensing data, and extracting spatiotemporal features using a deep neural network including a spatiotemporal convolution layer and an attention mechanism based on the standardized dynamic remote sensing data, including multispectral data, SAR data and DEM data;

[0022] An identification and storage module is used to identify key change areas through a dual attention module according to the extracted spatiotemporal features, store the spatiotemporal features and key change area information into a hierarchical dynamic feature library, and obtain dynamic change features of ground objects;

[0023] A construction module is used to construct an uncertainty distribution model using a Gaussian mixture model and a Beta distribution model according to the dynamic change characteristics of the ground object, and to design a probability constraint according to the uncertainty distribution model;

[0024] A generation module is used to calculate the minimum safe distance based on the uncertainty distribution model of the probability constraint in combination with the uncertainty ellipsoid, generate an adaptive safe distance rule, obtain a collision avoidance rule through verification, and establish a multi-objective optimization function of land use efficiency index and environmental protection index according to the collision avoidance rule;

[0025] The execution module is used to convert the collision avoidance rules into mathematical inequalities as constraints, adopt the NSGA-III algorithm to solve the multi-objective optimization problem, generate multiple non-dominated solutions, and execute the space planning strategy according to the multiple non-dominated solutions.

[0026] This application has the following technical effects:

[0027] By constructing a multi-source remote sensing data acquisition network based on the LPV model, the model contains a state matrix and a time-varying parameter vector that describe the changes of the ground objects, which can effectively capture the dynamic change characteristics of the ground objects. Combined with the spatiotemporal registration algorithm and the Kalman filter for noise filtering, the precise alignment and noise reduction of multi-source heterogeneous data are achieved, laying a reliable data foundation for subsequent feature extraction and analysis.

[0028] Based on the processed standardized dynamic remote sensing data, the spatiotemporal convolution layer in the deep neural network is used to extract multi-scale features, and the attention mechanism is used to highlight important spatiotemporal information. The dual attention module accurately identifies key change areas by considering the feature dependencies of the spatial dimension and the channel dimension at the same time. These extracted features are stored in a hierarchical dynamic feature library, forming a complete feature description system for dynamic changes of land objects.

[0029] Based on the dynamic change characteristics of the ground objects, the Gaussian mixture model is used to describe the location distribution, and the Beta distribution model is used to characterize the shape parameters to build a complete uncertainty distribution model. Based on the probability constraints and adaptive safety spacing rules designed by the model, the reliability and safety of the spatial planning scheme are ensured through the calculation and verification of the uncertainty ellipsoid.

[0030] The collision avoidance rules are transformed into mathematical constraints, and a multi-objective optimization function is constructed by combining land use efficiency and environmental protection indicators. The NSGA-III algorithm is used to generate multiple non-dominated solutions, providing decision makers with multiple optional planning schemes. This method not only ensures the feasibility of the planning scheme, but also achieves the optimal balance between efficiency and environmental protection.

[0031] In summary, high-quality dynamic remote sensing data is obtained through the LPV model and data processing technology, and then deep learning methods are used to extract accurate change features. Uncertainty models and safety constraints are constructed based on these features, and finally scientific spatial planning is achieved through multi-objective optimization. This solves the problems of insufficient data processing accuracy, inaccurate feature extraction, insufficient uncertainty consideration, and single optimization target in existing technologies, and significantly improves the scientificity and reliability of national land space planning. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] In order to more clearly illustrate the technical solution of the embodiment of the present application, the following is a brief introduction to the drawings required for use in the embodiment:

[0033] Figure 1 A flowchart of a land space planning and mapping method based on dynamic remote sensing technology provided in an embodiment of the present application;

[0034] Figure 2 This is a schematic diagram of the structure of a multi-source remote sensing data acquisition network in an embodiment of the present application;

[0035] Figure 3 The specific process of step S2 in the embodiment of this application Figure 1 ;

[0036] Figure 4 The specific process of step S2 in the embodiment of this application Figure 2 ;

[0037] Figure 5 The specific process of step S3 in the embodiment of this application Figure 1 ;

[0038] Figure 6 The specific process of step S3 in the embodiment of this application Figure 2 ;

[0039] Figure 7This is a structural diagram of a national land space planning and surveying system based on dynamic remote sensing technology in an embodiment of the present application. DETAILED DESCRIPTION

[0040] The embodiments of the present application are described in detail below with reference to the accompanying drawings.

[0041] like Figure 1 As shown, the embodiment of the present application provides a land space planning surveying and mapping method based on dynamic remote sensing technology, comprising the following steps:

[0042] S1: construct a multi-source remote sensing data acquisition network based on a linear parameter variation LPV model, wherein the LPV model includes a state matrix describing changes in ground objects and a time-varying parameter vector describing seasonal changes and climate conditions; through the multi-source remote sensing data acquisition network, obtain time-series remote sensing data including satellite remote sensing, UAV remote sensing and ground observation station remote sensing;

[0043] The core of the LPV model is a time-varying state space system, whose state equation can be expressed as x(t+1) = A(ρ(t))x(t) + B(ρ(t))u(t), where x(t) is the state vector, which contains the location, shape, spectral characteristics and other attributes of the object, u(t) is the external input vector, and A(ρ(t)) and B(ρ(t)) are the state transfer matrix and input matrix that depend on the time-varying parameter vector ρ(t), respectively. The time-varying parameter vector ρ(t) describes environmental factors such as seasonal changes and climate conditions, and its components can include normalized temperature index, precipitation index, vegetation index, etc.

[0044] The state transfer matrix A(ρ(t)) adopts the affine parameterization form, that is, , where A0 is the basic state matrix, describing the inherent change characteristics of the object, A i is a parameter matrix that characterizes the regulatory effect of environmental factors on the change of land features. Through this parameterized structure, the model can flexibly describe the dynamic change process of land feature characteristics with environmental conditions. For example, in farmland monitoring, A0 can represent the basic growth pattern of crops, while ρi(t)Ai reflects the impact of factors such as temperature and precipitation on the growth rate.

[0045] The parameter estimation adopts a two-step method based on subspace identification. First, the principal components of the time-varying parameter vector ρ(t) are extracted from multi-source remote sensing data using the principal component analysis method. Then, the parameters of the state matrix are estimated based on the least squares criterion. In order to improve the robustness of the estimation, a regularization term is introduced to constrain the norm of the parameter matrix, and the regularization coefficient is determined by cross-validation. The validation of the model is evaluated by comparing the mean square error and correlation coefficient between the predicted values ​​and the measured values. If necessary, the model performance can be optimized by increasing the number of parameter matrices or adjusting the regularization strength.

[0046] An important feature of the model is its ability to handle multi-scale dynamic features. At the local scale, the model describes the changing characteristics of a single feature through a state vector; at the regional scale, it captures the spatial distribution of environmental conditions through a time-varying parameter vector. This multi-scale modeling capability enables the model to show good adaptability in complex national land space planning applications. For example, the model can simultaneously characterize the local expansion pattern and regional development trend of urban built-up areas, providing comprehensive data support for planning decisions.

[0047] like Figure 2 As shown in Figure 1, the multi-source remote sensing data acquisition network integrates three types of data sources: satellite remote sensing, UAV remote sensing, and ground observation stations, forming a collaborative observation system. The LPV model describes the dynamic change characteristics of the ground objects through the state matrix, and introduces the time-varying parameter vector to characterize the influence of seasonal changes and climate conditions. The state matrix contains the basic characteristic information of the ground objects, while the time-varying parameter vector reflects the regulatory effect of environmental factors on the changes of the ground objects.

[0048] The time series remote sensing data acquired by the network has significant multi-source heterogeneous characteristics, including not only satellite multispectral data and SAR data, but also drone high-resolution images and ground observation data. There is a clear temporal and spatial correspondence between these data sources, which can complement and verify each other, improving the reliability of the data. At the same time, by continuously monitoring and obtaining a continuous observation sequence of ground object changes, the temporal continuity of the data is guaranteed.

[0049] It should be noted that the construction of this multi-source remote sensing data acquisition network has laid a solid data foundation for the subsequent spatiotemporal feature extraction and change analysis. The distributed architecture of the network ensures the spatial coverage of data acquisition, while the parameterized description based on the LPV model improves the ability to describe the dynamic changes of ground objects.

[0050] S2: performing spatiotemporal registration and noise filtering on the time series remote sensing data according to the state matrix and the time-varying parameter vector to obtain standardized dynamic remote sensing data, wherein the standardized dynamic remote sensing data includes multispectral data, SAR data and DEM data, and extracting spatiotemporal features according to the standardized dynamic remote sensing data using a deep neural network including a spatiotemporal convolution layer and an attention mechanism;

[0051] like Figure 3 As shown, the time-series remote sensing data is subjected to spatiotemporal registration and noise filtering according to the state matrix and the time-varying parameter vector to obtain standardized dynamic remote sensing data, including:

[0052] S2.1: constructing a state space model according to the time series remote sensing data, wherein a system matrix of the state space model depends on a time-varying parameter vector;

[0053] The state space model is constructed based on time-series remote sensing data. The system matrix of this model depends on the time-varying parameter vector, which can describe the dynamic change of the features of the ground objects over time and environmental conditions. The state vector contains key information such as the location, morphology and spectral characteristics of the ground objects, and the observation equation describes the imaging characteristics of different sensors.

[0054] Based on the time series remote sensing data acquired by the multi-source remote sensing data acquisition network, a state space model is constructed. The system matrix of the state space model depends on the time-varying parameter vector describing seasonal changes and climate conditions. Specifically, for the system matrix A, the system can be expressed as:

[0055] A0ρ0(t) + A1ρ1(t) + A2ρ2(t) + ... + Anρn(t)

[0056] Among them, A0 is the basic system matrix, Ai (i=1,2,...,n) is the parameter matrix, ρi(t) is the time-varying parameter, including environmental factors such as seasonal index, temperature, precipitation, etc. The state vector x(t) contains the attribute information of the location, morphology, spectral characteristics, etc. of the ground object.

[0057] S2.2: Based on the state space model, a spatiotemporal registration algorithm is used to calculate the transformation matrix between different data sources to obtain the registered remote sensing data;

[0058] Based on the constructed state space model, the space-time registration algorithm is used to calculate the transformation matrix between different data sources. The algorithm first extracts stable objects in each data source as control points, and then optimizes the transformation parameters through the iterative closest point algorithm, and finally obtains the accurate geometric transformation relationship. By applying these transformation matrices, remote sensing data from different sources are registered into a unified spatial reference system.

[0059] Based on the constructed state space model, the pyramid iterative nearest point algorithm is used in combination with geometric invariant features for spatiotemporal registration. The algorithm first establishes a four-layer pyramid structure for multi-source remote sensing data, extracts SIFT feature points at the coarsest resolution layer, and constructs feature descriptors through geometric invariant moments. Then, in each layer of the pyramid, the KD tree is used to accelerate the nearest neighbor search and establish the initial matching pairs of feature points. The RANSAC algorithm is used to eliminate incorrect matches for the matching pairs and optimize the transformation matrix parameters. In particular, the temporal consistency constraint of the object is introduced in the iterative process, requiring that the geometric deformation of the same object in different phase images should remain continuous.

[0060] Specifically, the pyramid structure of the spatiotemporal registration algorithm adopts a four-layer design. The base layer maintains the original resolution, and the resolution is reduced by half with each ascending layer. When constructing the pyramid, a Gaussian smoothing kernel is used for downsampling. The kernel size is set to 5×5, and the standard deviation increases with the increase in the level. This can suppress high-frequency noise while maintaining important structural features. In addition to storing image data, each layer of the pyramid also contains a gradient direction histogram and local descriptors to form a multi-scale feature expression. This hierarchical structure significantly improves the computational efficiency of the algorithm, especially when processing high-resolution remote sensing images.

[0061] The selection of feature points adopts a three-step strategy. First, the corner point response function is used to calculate candidate points in each layer of the pyramid, and the local extreme points are selected as the initial feature points. Then non-maximum suppression is performed to ensure the uniformity of the spatial distribution of feature points. Finally, based on the grayscale change pattern around the feature points, the points with significant distinction are screened and retained. The selection process of feature points takes into account time continuity, and gives priority to points that are stable in multiple phase images. These points usually correspond to fixed structures of the objects, such as building corners or road intersections.

[0062] The implementation of temporal consistency constraints is based on the motion model of feature points. For each feature point, a state vector containing position, velocity and acceleration is constructed, and its position in images of different time phases is predicted using a Kalman filter. If the deviation between the actual matching position and the predicted position exceeds a threshold, the match is considered to be incorrect. This dynamic constraint is particularly suitable for dealing with slowly changing features, such as the expansion of urban built-up areas. In practice, the choice of thresholds is dynamically adjusted according to the time interval and feature type.

[0063] In the feature description stage, the algorithm combines multiple local features. In addition to the traditional gradient direction histogram, texture features based on local binary patterns and shape features based on Zernike moments are also introduced. These features are extracted in different scale spaces and fused in a weighted manner to improve the discriminative ability of the features. The feature matching adopts a hierarchical strategy, first performing a rough match on the top layer of the pyramid, and then refining it layer by layer. Each layer uses a local search window to limit the matching range, which significantly improves the computational efficiency.

[0064] The optimization of transformation parameters adopts a robust estimation method. First, the RANSAC framework is used to remove outliers, and then the weighted least squares method is used to refine the parameters. The setting of weights takes into account the confidence and temporal consistency of the feature points. In order to improve the stability of the algorithm, regularization terms are introduced in the optimization process to prevent overfitting. The final transformation parameters are determined by cross-validation to ensure good generalization performance between different data sources.

[0065] Through this improved spatiotemporal registration algorithm, the system can effectively handle the registration problem of multi-source remote sensing data. For example, when processing Landsat and Sentinel data in a certain urban area, the algorithm can still maintain sub-pixel registration accuracy even in the presence of cloud cover and seasonal changes. This high-precision registration lays a solid foundation for subsequent change detection and feature extraction. Especially when monitoring rapidly urbanizing areas, accurate spatiotemporal registration is crucial to capturing the expansion process of built-up land.

[0066] The algorithm adds three key improvements to the standard ICP algorithm: first, the pyramid layering strategy is used to gradually refine the registration accuracy to avoid falling into the local optimum; second, the SIFT invariant feature is combined to improve the stability of feature matching, making the algorithm more adaptable to the registration needs of multi-source heterogeneous data; third, the temporal consistency constraint is used to improve the registration reliability in dynamic scenes. Experiments show that the improved algorithm improves the registration accuracy of complex terrain areas by more than 30% compared with traditional methods, and reduces the registration time by 40%.

[0067] Specifically, based on the constructed state space model, the pyramid iterative closest point algorithm is used in combination with geometric invariant features for spatiotemporal registration. The algorithm first constructs a four-layer pyramid structure for multi-source remote sensing data, with a resolution ratio of 1:2:4:8 for each layer to balance computational efficiency and accuracy requirements. The improved SIFT algorithm is used to extract feature points at the coarsest resolution layer, which mainly includes three improvements: replacing the Gaussian difference operator with the LoG operator to improve the positioning accuracy of feature points; using an adaptive threshold strategy to select feature points to ensure that feature points are evenly distributed in the image; introducing a feature point screening mechanism based on the consistency of the main direction to improve the stability of feature points.

[0068] In the feature description stage, the algorithm constructs an invariant feature descriptor that integrates geometric and radiation information. Specifically, in the 128×128 neighborhood around the feature point, the pixel gradient direction histogram is calculated as the basic descriptor, and the local geometric moment feature and texture entropy feature are introduced to form a 384-dimensional enhanced descriptor. In order to improve the efficiency of feature matching, local sensitive hashing (LSH) is used to map the high-dimensional descriptor to a low-dimensional space, and the KD tree structure is used to accelerate the nearest neighbor search. Experiments show that the enhanced descriptor has a 25% higher accuracy in multi-source heterogeneous data matching than the traditional SIFT descriptor.

[0069] In the feature matching stage, the algorithm adopts a three-step strategy to ensure the reliability of matching. First, the first three nearest neighbor candidate points of each feature point are quickly retrieved based on the KD tree, and the nearest neighbor ratio is calculated to screen the initial matching pairs. Then, the improved RANSAC algorithm is used to eliminate false matches. The improvements include: adaptive sampling strategy, dynamically adjusting the sampling probability according to the spatial distribution of matching point pairs; multi-model parallel verification, considering both affine transformation and perspective transformation models; local consistency constraints, requiring that the transformation parameters of adjacent matching point pairs should be similar. Finally, the temporal consistency constraint is introduced. By analyzing the transformation matrix in the continuous time series, the temporal smoothing constraint of the transformation parameters is established to effectively suppress abnormal registration results.

[0070] In the transformation matrix optimization stage, the algorithm iteratively performs the following steps in each pyramid layer: first, the feature points are projected based on the current transformation matrix; then, the transformation parameters are refined using the improved ICP algorithm, which includes: weighted least squares estimation, which assigns weights according to the reliability of matching point pairs; adaptive step size strategy, which dynamically adjusts the iteration step size according to the convergence speed; robustness evaluation, which evaluates the uncertainty of transformation parameters through bootstrap resampling. In particular, the covariance propagation theory is used for parameter transfer between different pyramid levels to ensure the smooth transition of transformation parameters.

[0071] Through a large number of experimental verifications, the improved registration algorithm has shown excellent performance in different scenarios. In mountainous and hilly areas, the average registration accuracy of the algorithm is better than 0.5 pixels, which is better than the 1.2 pixels of the traditional method; in urban built-up areas, the algorithm can effectively handle the geometric deformation caused by the height difference of buildings, and the registration accuracy is better than 0.8 pixels; in farmland areas, the algorithm has strong adaptability to seasonal changes in land objects, and the stability of temporal registration is improved by 40%. In addition, the computational efficiency of the algorithm has also been significantly improved. It takes only 2-3 seconds on average to process a 1000×1000 pixel image pair, which is more than 50% faster than the traditional method.

[0072] Take the multi-source remote sensing data registration of a certain urban area as an example: First, 10 obvious corner points of artificial buildings were extracted from high-resolution satellite images as the main reference control points, and the control points at the same positions were identified in the drone images and ground observation data. In the registration process, the iterative closest point algorithm is used to optimize the transformation parameters by minimizing the distance difference between the control points. For example, the corner points of a building with coordinates (116.123, 39.456) may have deviations of several meters to tens of meters in different data sources. These deviations are corrected to centimeter-level accuracy through the calculated transformation matrix. In particular, for SAR data, considering its side-view imaging characteristics, terrain undulation correction can also be performed to ensure the registration accuracy in mountainous areas.

[0073] S2.3: for the registered remote sensing data, using a Kalman filter to perform adaptive noise suppression to obtain filtered data;

[0074] For the registered remote sensing data, an adaptive Kalman filter is used to suppress noise. The filter can dynamically adjust the process noise covariance and observation noise covariance according to the signal-to-noise ratio of the data to effectively suppress different types of noise. At the same time, the filter also integrates an outlier detection mechanism, which can identify and process mutation points, ensuring the reliability of the filtering results.

[0075] S2.4: updating the state matrix describing the change of the ground object in the LPV model according to the filtered data;

[0076] The state matrix describing the changes of the objects in the LPV model is updated according to the filtered data. The updating process adopts a sliding time window strategy and optimizes the model parameters through the least squares estimation method, so that the state matrix can more accurately reflect the current state and change trend of the objects.

[0077] Take the state matrix update of a farmland area within one year as an example: the initial state matrix records the crop planting situation, vegetation coverage and other information in the area. With the continuous acquisition of time series remote sensing data, statistics show that the region has an obvious vegetation index growth trend in spring (March-May), indicating that crops have entered the growth period. The data was analyzed through a sliding time window (the window size was set to 30 days) to timely update the vegetation parameters in the state matrix. At the same time, if local farmland was detected to be flooded during the flood season (July-August), the land use type parameters and water coverage parameters of the corresponding area were adjusted accordingly, so that the state matrix can accurately reflect the dynamic change process of the land objects. This timely parameter update ensures the timeliness and accuracy of subsequent spatial planning decisions.

[0078] S2.5: standardizing the filtered data based on the updated state matrix to obtain standardized dynamic remote sensing data;

[0079] The filtered data is standardized based on the updated state matrix. The standardization process takes into account the radiation characteristics and imaging conditions of different sensors, and uses a piecewise linear transformation method to normalize the data to a unified value range. At the same time, the spectral response function is introduced for band matching to ensure the physical consistency of multispectral data, SAR data, and DEM data.

[0080] Through the above processing steps, a standardized dynamic remote sensing data set is finally obtained. These data not only have precise spatial position correspondence and effectively suppressed noise, but also have a unified data format and metric standard, providing high-quality data support for subsequent feature extraction and analysis.

[0081] The filtered data is standardized based on the updated state matrix, and a hierarchical adaptive normalization strategy is adopted. First, for multispectral data, a radiation normalization method based on statistical characteristics is adopted. By calculating the histogram features of each band, the main peak and valley are extracted as reference points, and the pixel values ​​are mapped to the [0,1] interval using piecewise linear stretching. In particular, for spectral combination indices such as vegetation index, seasonal variation factors are introduced for dynamic correction to ensure the comparability of data in different phases. For example, in the vegetation growing season, the normalization parameters are dynamically adjusted according to the spectral response characteristics of the reference objects, so that the standardized vegetation index can accurately reflect the vegetation growth status.

[0082] For SAR data, a multi-step standardization process is adopted. First, radiation calibration is performed to convert the backscatter coefficient into decibel units; then adaptive quantization is performed based on local statistical characteristics, and the scattering characteristics of each sub-region are calculated using the partition statistics method to determine the optimal quantization parameters. Special attention is paid to strong scattering points and dark areas in SAR images, and the information of these areas is maintained through nonlinear mapping functions. For example, in urban areas, the double scattering and angular reflection characteristics of buildings are standardized using piecewise logarithmic functions, which not only maintains the relative intensity of the scattering characteristics, but also compresses the dynamic range of the data.

[0083] For DEM data, adaptive standardization considering terrain characteristics is implemented. First, the terrain undulation characteristics of the study area are analyzed, and terrain factors such as slope and aspect are calculated; then, based on the results of geomorphic zoning, different standardization strategies are adopted for different geomorphic units. In plain areas, linear normalization is used; in mountainous areas, nonlinear normalization based on elevation distribution characteristics is adopted to ensure that the standardized DEM data can maintain the relative height difference relationship of the terrain. At the same time, considering the spatial continuity of DEM data, a smooth transition function is introduced at the boundary of geomorphic units to avoid mutations in the standardized results.

[0084] In order to ensure the consistency of the physical meaning of multi-source data, the spectral response function is introduced for band matching. In view of the differences in spectral response of different sensors, a spectral conversion model is established to convert the data of different sensors into a unified spectral reference system. For example, when processing Landsat and Sentinel satellite data, the spectral response functions of the two sensors are analyzed to establish the conversion relationship between bands to ensure the comparability of standardized multi-source data. At the same time, considering the influence of atmospheric conditions and observation geometry, the BRDF correction model is introduced to eliminate the radiation differences caused by differences in observation angles.

[0085] The standardization process also takes into account the consistency of the spatial and temporal scales of the data. In the spatial dimension, the data of different spatial resolutions are unified into the same geographic grid system through the resampling method; in the temporal dimension, a continuous time series is constructed based on the time series interpolation method, and high-frequency noise is eliminated through time window filtering. The final output standardized data set not only maintains the physical characteristics of the original data, but also has a unified data format and metric standard, providing a reliable data foundation for subsequent feature extraction and analysis.

[0086] like Figure 4 As shown, according to the standardized dynamic remote sensing data, a deep neural network including a spatiotemporal convolution layer and an attention mechanism is used to extract spatiotemporal features, including:

[0087] S2.6: receiving the standardized dynamic remote sensing data, and performing normalization preprocessing on the multispectral data, the SAR data, and the DEM data respectively through a multi-branch input layer to obtain preprocessed multi-source data;

[0088] S2.6 adopts a multi-branch parallel processing architecture and designs special preprocessing processes for different types of remote sensing data. For multispectral data, atmospheric correction is first performed to eliminate the effects of atmospheric scattering and absorption, and then histogram equalization is performed on each band to enhance image contrast. Next, the minimum-maximum normalization method is used to map pixel values ​​to the [0,1] interval to ensure the scale consistency of data in different bands. Taking Landsat 8 satellite data as an example, a total of 7 bands, including visible light, near infrared and short-wave infrared, are normalized separately to make them comparable.

[0089] For SAR data, considering its scattering characteristics and coherent noise, speckle noise suppression is first performed, and an improved Lee filter is used to maintain edge details while reducing noise. Then the backscatter coefficient is logarithmically transformed to convert multiplicative noise into additive noise for subsequent processing. Finally, the Z-score standardization method is used to convert the data into a standard normal distribution with a mean of 0 and a standard deviation of 1, thereby improving the comparability of SAR data at different phases.

[0090] For DEM data, we first fill holes and correct outliers to ensure the continuity of elevation data. Then, according to the terrain characteristics of the study area, we calculate the statistical distribution of elevation data, and use piecewise linear normalization to map elevation values ​​to the [0,1] interval while maintaining the relative relationship of terrain undulations. In addition, terrain factors such as slope and aspect are calculated from DEM data, and these derived products are also normalized.

[0091] S2.7: constructing a 3D convolution layer for the preprocessed multi-source data, setting the parameters of the convolution kernel size to 3 and the step size to 1, and obtaining preliminary spatiotemporal features;

[0092] In step S2.7, an innovative 3D convolutional neural network structure is constructed to extract spatiotemporal features. First, the preprocessed multi-source data is organized into a four-dimensional tensor in chronological order, where three dimensions correspond to spatial position (x, y) and time (t), and the fourth dimension contains feature channels from different data sources. For example, for a data slice at a certain time point, it contains 7 multispectral bands, 2 SAR polarization channels, 1 DEM elevation value and 2 terrain factors, a total of 12 feature channels.

[0093] A multi-layer 3D convolution structure is used. The first layer has 64 convolution kernels. The size of each convolution kernel is 3×3×3 (corresponding to the three dimensions of space and time), the step size is 1, and the padding method uses the "same" mode to keep the feature map size unchanged. This design enables each convolution kernel to capture the feature patterns in the spatial neighborhood and time series at the same time. For example, for farmland areas, 3D convolution can simultaneously extract the spatial distribution characteristics and growth cycle characteristics of crops.

[0094] After the convolution operation, the ReLU activation function is used to introduce nonlinear transformation, and the batch normalization layer is used to stabilize the training process. After this processing, the preliminary spatiotemporal features contain rich local spatiotemporal pattern information, laying the foundation for subsequent feature enhancement and change detection.

[0095] It is particularly worth noting that the receptive field design of 3D convolution fully considers the temporal and spatial characteristic scales of different objects. Taking building change detection as an example, the 3×3×3 convolution kernel can effectively capture the basic structural units and short-term change characteristics of the building; and through the superposition of multiple layers of convolution, the receptive field can be gradually expanded to achieve feature extraction of large-scale building complexes and long-term urban expansion processes.

[0096] In one embodiment, the deep neural network adopts an encoder-decoder structure, including 5 encoding blocks and corresponding 5 decoding blocks. At the input end of the network, multispectral, SAR and DEM data are first processed through independent preprocessing branches. Each encoding block contains two 3D convolutional layers, the first layer performs feature dimensionality reduction, and the second layer restores the feature dimension and extracts spatial information. The initial input layer is set to 64 channels. As the network depth increases, the number of channels increases to 128, 256, 512 and 1024 respectively. This progressive channel expansion helps capture multi-scale features from details to the global.

[0097] The design of the residual block is optimized for spatiotemporal feature extraction and contains three parallel processing branches. The main branch extracts local spatiotemporal features through two consecutive 3D convolutional layers; the shortcut branch maintains the original information of the input features; and the attention branch is responsible for adaptive weighting of the feature map. This multi-branch structure ensures full information extraction while avoiding network degradation problems. There are a total of 4 such residual blocks in the network, which are located after the deeper encoding layer. This configuration can improve the accuracy of feature extraction while maintaining computational efficiency.

[0098] The design of the temporal attention module adopts a multi-head attention mechanism, using 8 independent attention heads for feature extraction, and the output dimension of each head is set to 64. This design allows the model to focus on different feature patterns at the same time, significantly improving the comprehensiveness of feature extraction. The attention module first calculates the feature correlation in the time dimension, and then generates attention weights based on the correlation to adaptively weight the features. This mechanism is particularly suitable for processing the dynamic change characteristics of objects and can accurately capture the change information at key time points.

[0099] The decoder part adopts a symmetrical structural design, and gradually restores the spatial resolution of the feature map through the deconvolution layer. In addition to processing the features from the previous layer, each decoding block also integrates the feature information of the corresponding encoding layer. This jump connection design effectively combines deep semantic information and shallow detail features. During the decoding process, the number of channels is gradually reduced from 1024 to 64, and finally the prediction results are output through the convolution layer. The entire decoding process maintains the spatial continuity of the features and effectively avoids information loss.

[0100] The training of the network adopts a phased strategy, first training the backbone network with a larger learning rate, and then fine-tuning the attention module with a smaller learning rate. During the training process, a variety of data enhancement methods are combined to enhance the generalization ability of the model, including random cropping, flipping and rotation. At the same time, dropout and regularization mechanisms are introduced to prevent overfitting. This carefully designed network architecture has demonstrated powerful spatiotemporal feature extraction capabilities in practical applications, and can accurately identify and analyze the dynamic change characteristics of various types of objects.

[0101] S2.8: For the preliminary spatiotemporal features, multi-scale features are obtained by adding a residual block including two 3D convolutional layers and a batch normalization layer;

[0102] The embodiment of the present application extracts multi-scale features through a deep residual learning structure. Each residual block contains two 3D convolutional layers connected in series. The first convolutional layer uses 64 1×1×1 convolution kernels for feature dimensionality reduction to reduce computational complexity; the second convolutional layer uses 64 3×3×3 convolution kernels for feature extraction. Between these two convolutional layers, a batch normalization layer and a ReLU activation function are inserted to effectively alleviate the gradient vanishing problem of the deep network. In particular, the residual block adds the input features directly to the output features through a jump connection. This design not only retains the shallow detail information, but also improves the training efficiency of the network.

[0103] In order to obtain feature representations at different scales, the embodiment of the present application connects multiple residual blocks in series. For example, for monitoring changes in urban built-up areas, the first residual block mainly extracts the local structural features of buildings, the second residual block focuses on the spatial organization features at the block scale, and the third residual block further extracts the macroscopic distribution features of urban functional areas. Through this multi-scale feature extraction mechanism, the fine texture of the objects and the regional development trend can be captured simultaneously.

[0104] S2.9: For the multi-scale features, a temporal attention module is set to calculate the attention weights of the features at different time points to obtain a weighted spatiotemporal feature map;

[0105] In the embodiment of the present application, an innovative temporal attention mechanism is introduced to enhance the feature representation of key time points. Specifically, for a certain spatial position, the embodiment of the present application first calculates the similarity matrix between the feature vectors of the position at different time points. The similarity calculation adopts the scaled dot product method, that is, the dot product of the query vector and the key vector is divided by the square root of the feature dimension, which can stabilize the gradient and improve the training effect.

[0106] Based on the calculated similarity matrix, the embodiment of the present application generates attention weights through the Softmax function, so that historical moments that have a greater impact on the current state receive higher weights. For example, when monitoring changes in farmland cultivation, historical periods close to the current growth period will receive higher attention weights, while features of fallow periods or abnormal climate periods will receive lower weights. This adaptive weighting mechanism can effectively highlight the features of key time points and suppress interference from irrelevant periods.

[0107] The output of the temporal attention module is a weighted feature map, which not only retains the detailed information of the spatial dimension, but also enhances the key change characteristics of the temporal dimension. This organic fusion of spatiotemporal features provides a more reliable feature representation for subsequent change detection and trend analysis. For example, in urban expansion analysis, it can not only accurately locate the spatial location of newly added construction land, but also analyze the temporal change pattern of construction intensity through attention weights.

[0108] S3: Based on the extracted spatiotemporal features, the dual attention module is used to identify key change regions (e.g. Figure 4 As shown in FIG. 1 ), the spatiotemporal characteristics and key change area information are stored in a hierarchical dynamic feature library to obtain dynamic change characteristics of the ground object;

[0109] like Figure 5 As shown in the figure, the key change areas are identified by the dual attention module based on the extracted spatiotemporal features, including:

[0110] S3.1: For the weighted spatiotemporal feature map, calculate the correlation matrix between each position of the spatiotemporal feature map, and use the softmax function to generate spatial weights to obtain a spatial attention map;

[0111] In step S4.1, a spatial attention module is constructed to capture the spatial dependencies between different locations of the feature map. Specifically, for a spatiotemporal feature map of C×H×W dimensions (C is the number of channels, H and W are the height and width of the feature map, respectively), the system first reshapes it into a correlation matrix of (C×H×W)×(C×H×W). For example, in the monitoring of changes in urban construction land, the features of the location of a newly built community will show a strong correlation with the features of the surrounding infrastructure locations. The correlation matrix is ​​normalized by the softmax function, and the spatial attention weights are generated, so that areas with significant changes obtain higher weight values.

[0112] S3.2: For the weighted spatiotemporal feature map, global average pooling and a fully connected layer are used to calculate the dependency between channels to obtain a channel attention map;

[0113] In step S3.2, the dependency between feature channels is modeled through the channel attention module. First, global average pooling is performed on each channel, and the C×H×W feature map is compressed into a C×1×1 vector to extract the global response of each channel. Then, through a two-layer fully connected network, the first layer uses the ReLU activation function to reduce the feature dimension (to 1 / 16 of the original number of channels), and the second layer performs feature reconstruction to restore the original number of channels. This design can effectively learn the importance of different feature channels. For example, in vegetation change monitoring, the near-infrared band channel may receive a higher attention weight, while the visible light band has a relatively low weight.

[0114] S3.3: Combining the spatial attention map and the channel attention map, weighting the weighted spatiotemporal feature map in both the spatial dimension and the channel dimension to obtain a key change area;

[0115] In step S3.3, the spatial attention and channel attention are fused to achieve double weighting of the feature map. In specific implementation, the spatial attention weight and channel attention weight are first multiplied by the original feature map respectively, and then the two weighted results are added element-wise to obtain a feature representation that considers both spatial correlation and channel dependency. This dual weighting mechanism can locate the change area more accurately. For example, in land use change monitoring, not only the intensity of change in spatial position is concerned, but also the change characteristics of different land object types (corresponding to different channels) are considered.

[0116] S3.4: Based on the key change area, the spatial attention map and the channel attention map are fused to obtain a change significance map;

[0117] In step S3.4, the final change saliency map is generated through an adaptive fusion strategy. First, the complementarity of the spatial attention map and the channel attention map is calculated, and the fusion weights of the two attentions are dynamically adjusted according to the degree of complementarity. Then the weighted attention map is deeply fused with the key change area features to obtain a saliency map that highlights the degree of change. This fusion strategy is particularly suitable for dealing with complex change scenarios. For example, in the monitoring of land use changes in urban and rural fringe areas, it can simultaneously reflect the spatial pattern of construction land expansion and the conversion characteristics of land use types.

[0118] Each pixel value of the change saliency map indicates the possibility of change at that location. The larger the value, the more significant the change. The system can also set an appropriate threshold to segment the saliency map according to actual application requirements to obtain a binary change detection result. This change detection method based on dual attention not only improves the accuracy of detection, but also maintains good noise resistance.

[0119] like Figure 6 As shown, the spatiotemporal features and key change area information are stored in a hierarchical dynamic feature library to obtain the dynamic change features of the ground objects, including:

[0120] S3.5: Build a hierarchical dynamic feature library including basic feature layer, change pattern layer and knowledge graph layer;

[0121] A three-layer dynamic feature library is designed to achieve progressive expression from low-level features to high-level semantics. The basic feature layer is responsible for storing and indexing original feature information; the change pattern layer summarizes and summarizes typical patterns of change; and the knowledge graph layer builds a network of associations between change events. This hierarchical design enables the system to understand and express the process of land feature change at different levels of abstraction, such as from the numerical changes of vegetation indices, to the evolutionary patterns of vegetation coverage, and then to the causal relationship of ecosystem succession.

[0122] S3.6: Based on the weighted spatiotemporal feature map and the change significance map, storing the original feature vector in the basic feature layer and establishing a feature index;

[0123] The weighted spatiotemporal feature map and change significance map are projected into the feature space to generate a vector describing the change characteristics of each spatial position. For efficient retrieval, local sensitive hashing (LSH) can be used to construct a feature index structure. In specific implementation, the feature vector is divided into multiple subspaces, and each subspace uses an independent family of hash functions, which not only ensures retrieval efficiency but also maintains the local similarity of features. For example, for urban expansion monitoring, the system can quickly retrieve historical cases similar to the current construction land expansion pattern.

[0124] S3.7: Based on the data of the basic feature layer, a density clustering algorithm is used to identify the change pattern, and the ground feature change pattern is stored in the change pattern layer;

[0125] The density-based spatial clustering algorithm (DBSCAN) is used to cluster the data of the basic feature layer to identify typical change patterns. The algorithm first sets the neighborhood radius and the minimum sample number threshold, and then classifies similar change features into one category based on the density accessibility principle. For example, in farmland change monitoring, typical patterns such as "cultivated land-construction land conversion" and "paddy field-dry land conversion" may be identified. For each change pattern, its spatiotemporal distribution characteristics and evolution laws are recorded.

[0126] S3.8: Based on the analysis results of the change pattern layer, store information representing the causal relationship of the change in the knowledge graph layer;

[0127] Construct a knowledge graph to express the causal relationship between change events. The nodes in the graph represent change events or influencing factors, and the edges represent the relationship between them. By analyzing the temporal correlation and spatial adjacency of the change pattern, possible causal relationships can be inferred. For example, in the analysis of urban expansion, a causal chain such as "transportation facility construction → surrounding commercial development → residential land expansion" was found. The knowledge graph is stored in the form of triples to support complex relationship queries and reasoning.

[0128] S3.9: For the multi-layer feature structure of the hierarchical dynamic feature library, dynamically maintain feature distribution through an incremental update mechanism based on similarity calculation to obtain dynamic change features of the ground objects;

[0129] In step S3.9, the incremental update mechanism of the feature library is implemented. When new observation data arrives, the similarity between the new feature and the existing feature is first calculated. For the basic feature layer, the cosine similarity is used to measure the distance of the feature vector; for the change pattern layer, the edit distance is used to calculate the similarity of the change sequence; for the knowledge graph layer, the graph structure similarity is used to evaluate the matching degree of the relationship network. Based on these similarity calculation results, the feature distribution is dynamically adjusted and the features of each layer are updated.

[0130] For example, in wetland ecological monitoring, if a new degradation pattern is detected, the corresponding feature vector is added to the basic feature layer, the category definition of the degradation pattern is updated in the change pattern layer, and the association with factors such as climate change and human activities is supplemented in the knowledge graph layer. This incremental update mechanism enables the feature library to continuously accumulate and optimize the understanding of land feature changes.

[0131] Through the design of this hierarchical dynamic feature library, the multi-dimensional expression and dynamic maintenance of the changing features of the land objects are realized, providing a reliable knowledge basis for subsequent change analysis and prediction. The feature library not only records the basic information of "what" and "where", but also answers the deep questions of "why" and "how" it changes.

[0132] S4: According to the dynamic change characteristics of the ground object, a Gaussian mixture model and a Beta distribution model are used to construct an uncertainty distribution model, and probability constraints are designed according to the uncertainty distribution model;

[0133] When constructing the uncertainty distribution model, the dynamic change characteristics of the objects are first preprocessed and dimensionally reduced. For the location characteristics, the coordinate sequence of the center point and boundary point of the objects is extracted, and its change trajectory in the time series is calculated. The location distribution of each object usually shows multiple clustering centers, which is consistent with the multi-core characteristics of urban development. Therefore, a Gaussian mixture model containing 3 to 5 components is selected, and the optimal number of components is determined through iterative optimization. For each Gaussian component, its mean represents the cluster center, and the covariance matrix describes the discrete degree and main direction of the spatial distribution.

[0134] The Beta distribution model is mainly used to describe the shape parameters of land features, including area ratio, perimeter ratio, compactness and other characteristics. These parameters are first mapped to the [0,1] interval through maximum and minimum normalization. Considering the correlation between different shape parameters, multidimensional Beta distribution is used for joint modeling. When estimating parameters, the mean and variance of the sample are first calculated, and then the shape parameters α and β are solved by the moment estimation method. In order to improve the stability of the estimation, prior knowledge is introduced to constrain the range of parameter values, and the optimal prior parameters are selected through cross-validation.

[0135] The boundary fuzziness is described by a modified S-type membership function, which takes a value close to 1 in the core area of ​​the feature and exhibits a smooth decay characteristic in the transition area. The determination of the core area is based on the stability analysis of the feature in the time series, and the areas with frequent changes are assigned lower membership values. The parameters of the membership function are determined by minimizing the classification error of the boundary pixels, while considering the spatial continuity constraints to avoid mutations. This fuzzy boundary description is particularly suitable for dealing with the dynamic changes of natural features such as wetlands and farmlands.

[0136] The construction of the dynamic Bayesian network adopts a hierarchical approach. The bottom-level nodes represent the basic observation data, the middle-level nodes correspond to the distribution model of the position and shape parameters, and the top-level nodes describe the type and change pattern of the ground objects. The network structure is determined by the conditional independence test, and the conditional probability between nodes is determined by combining historical data statistics and expert knowledge. In order to deal with time series correlation, a time sliding window mechanism is introduced to regularly update the network parameters so that the model can adapt to the dynamic changes of ground object characteristics.

[0137] The uncertainty distribution model is evaluated using a multi-index system, including likelihood, information criterion, and prediction error. A large amount of simulation data is generated through Monte Carlo sampling to verify the performance of the model in different scenarios. The experimental results show that the model can accurately characterize the uncertainty distribution of land features and provide a reliable probabilistic framework for subsequent spatial planning decisions. For example, in urban expansion analysis, the model can not only predict the possible expansion direction of construction land, but also give a confidence interval for the expansion intensity, which is of great significance for formulating flexible planning strategies.

[0138] According to the dynamic change characteristics of the ground objects, the Gaussian mixture model and Beta distribution model are used to construct an uncertainty distribution model, including:

[0139] S4.1: for the dynamic change characteristics of the ground object, the Gaussian mixture model parameters describing the distribution of the ground object position are estimated by using the EM algorithm to obtain a probability distribution model of the ground object position;

[0140] In step S4.1, the expectation maximization (EM) algorithm is used to estimate the parameters of the Gaussian mixture model of the location distribution of the objects. First, the appropriate number of Gaussian components is determined according to the dynamic change characteristics of the objects, and the mean vector and covariance matrix of each component are initialized. In the E step, the posterior probability that each observed data point belongs to each Gaussian component is calculated; in the M step, the model parameters are updated based on these posterior probabilities. This iterative process continues until the parameters converge or the maximum number of iterations is reached. For example, in the analysis of urban built-up area expansion, it is found that the expansion area presents multiple growth cores, and the spatial distribution of each core can be described by a Gaussian component. Through the EM algorithm, the position (mean) and expansion range (covariance) of each growth core can be accurately estimated.

[0141] S4.2: Based on the probability distribution model, a multidimensional Beta distribution model describing the shape parameters of the land object is established, and a membership function is introduced to represent the boundary fuzziness to form an uncertainty description of the land object characteristics;

[0142] In step S4.2, a multidimensional Beta distribution model is constructed to describe the uncertainty of the shape parameters of the features. For each feature object, shape parameters such as area, perimeter, and compactness are extracted and mapped to the interval [0,1]. Then a Beta distribution model is established for each shape parameter, and its parameters α and β are obtained by maximum likelihood estimation. In order to deal with the ambiguity of the feature boundary, a distance-based fuzzy membership function is introduced. This function takes a value close to 1 in the core area of ​​the feature and smoothly transitions to 0 in the edge area as the distance increases. For example, in the extraction of wetland boundaries, it was found that some areas have obvious seasonal changes. The Beta distribution is used to describe the range of their area changes, and the membership function is used to express the gradual characteristics of the boundary.

[0143] S4.3: Describe the uncertainty of the features of the land object, construct a dynamic Bayesian network, and obtain an uncertainty distribution model;

[0144] In step S4.3, uncertainty is integrated and transmitted through a dynamic Bayesian network. The nodes of the network include the Gaussian component parameters of the position distribution, the Beta distribution parameters of the shape parameters, and the boundary membership function parameters. The edges between the nodes represent the probabilistic dependencies between the parameters, which are quantified by the conditional probability table (CPT). The network structure and parameters are updated in a rolling time window manner so that they can adapt to the dynamic changes of the features of the land. For example, in the monitoring of farmland changes, the dynamic Bayesian network can describe the changing patterns of the shape and area of ​​the field during the crop growth period, while taking into account the influence of factors such as weather and irrigation.

[0145] This multi-level uncertainty modeling not only quantifies the degree of uncertainty of land features, but also tracks the propagation of uncertainty in time and space. For example, in monitoring coastline changes, the position uncertainty caused by tidal influence (described by Gaussian mixture model), seasonal changes in beach morphology (described by Beta distribution), and ambiguity of the land-sea transition zone (described by membership function) can be considered simultaneously. These uncertainty information are integrated through dynamic Bayesian networks, providing a reliable probabilistic framework for subsequent change analysis and prediction.

[0146] It is particularly worth noting that an adaptive strategy is adopted in the parameter estimation process. For the Gaussian mixture model, the Bayesian Information Criterion (BIC) is used to automatically select the optimal number of components; for the Beta distribution, kernel density estimation is introduced to verify the rationality of the distribution assumption; for the membership function, the function parameters are adjusted based on local spatial autocorrelation analysis. These adaptive mechanisms enable the uncertainty model to better adapt to the characteristic changes of different types of land objects.

[0147] Designing probability constraints according to the uncertainty distribution model includes:

[0148] S4.4: Based on the prediction results of the uncertainty distribution model, a collision risk probability threshold ε is set for each ground object, where the value of ε ranges from 0.01 to 0.05, and basic constraint conditions are established;

[0149] In step S4.4, a personalized collision risk probability threshold is set for each object. The threshold is set based on factors such as the importance of the object, the frequency of change, and the spatial relationship. For example, for important urban infrastructure, a lower risk threshold (such as ε=0.01) is set, indicating a high requirement for change prediction; while for areas covered by natural vegetation, a relatively loose threshold (such as ε=0.05) is used. This differentiated threshold setting constitutes the basic constraint condition and provides an initial reference for subsequent constraint optimization. In particular, the threshold will be dynamically adjusted according to the temporal stability of the object, and the constraints will be appropriately relaxed for areas with frequent changes, while the constraints will be tightened for stable areas.

[0150] S4.5: For the basic constraint conditions, a Copula function is used to construct a joint probability constraint matrix among multiple ground objects;

[0151] In step S4.5, the Copula function is used to describe the dependency relationship between multiple features. First, the marginal probability distribution of each feature is mapped to the [0,1] interval through the probability integral transformation. Then, according to the spatial correlation between the features, the appropriate Copula family (such as Gaussian Copula, Clayton Copula, etc.) is selected to construct the joint distribution function. For example, in the analysis of land use changes, the changes in adjacent plots often show significant correlation, and this spatial dependency relationship can be accurately described by the Copula function. The joint probability constraint matrix not only records the dependency strength between the two features, but also contains the multivariate dependency structure.

[0152] When constructing joint probability constraints, the appropriate Copula function family is selected according to the correlation structure of the land feature change characteristics. For urban construction land with symmetric correlation, Gaussian Copula is used to describe its dependency structure, which can capture linear correlation characteristics well. For ecological land changes with tail correlation, ClaytonCopula is used, which can better describe the dependency relationship in extreme cases. In particular, for complex land use change patterns, a hierarchical Copula structure is used, first using Gaussian Copula to describe the intra-group correlation, and then using ClaytonCopula to characterize the inter-group dependency.

[0153] The parameter estimation of the Copula function adopts a two-stage method. First, the marginal distribution of each feature is estimated nonparametrically, and the empirical distribution function is obtained by the kernel density estimation method. Then, based on the transformed uniformly distributed data, the Copula parameters are estimated using the maximum likelihood method. In order to improve the stability of the estimation, the Bootstrap resampling technique is introduced to obtain the confidence interval of the parameter through multiple sampling. This method is particularly suitable for dealing with small sample situations and can reliably estimate the dependency structure between feature changes.

[0154] Wasserstein distance is used to evaluate the difference between the predicted distribution and the actual distribution. In practice, the second-order Wasserstein distance is mainly used, which not only considers the redistribution cost of probability mass, but also the influence of spatial position. In the specific calculation, the continuous distribution is first discretized into a probability vector, and then a cost matrix is ​​constructed, whose elements represent the transition cost between different states. For spatial distribution, this cost is usually defined as the square of the Euclidean distance. The Wasserstein distance between two distributions is obtained by solving the optimal transmission problem.

[0155] During dynamic monitoring, the system continuously calculates the Wasserstein distance between the predicted distribution and the observed distribution. When the distance exceeds a preset threshold, the model update mechanism is triggered. The threshold is selected based on historical data analysis and is usually set to the 95th percentile of the normal fluctuation range. This distance-based adaptive update strategy enables the model to respond promptly to significant changes in the characteristics of the ground object change, while avoiding overreaction to noise disturbances.

[0156] The uncertainty range is determined using a fuzzy set construction method based on the Wasserstein distance. First, a spherical neighborhood is constructed around the predicted distribution, and the radius of the sphere is proportional to the Wasserstein distance. This neighborhood is then converted into a fuzzy set using the level set method, and the membership value decreases as the distance from the predicted distribution increases. This construction method takes into account the geometric characteristics of the distribution, making the description of the uncertainty range more in line with the actual situation. For example, when monitoring urban expansion, the uncertainty range will extend along the main traffic axis, which is consistent with the actual law of urban development.

[0157] The model is validated using a multi-scenario testing method. Different change scenarios are constructed using historical data, and the Wasserstein distance between the predicted distribution and the actual distribution is compared to evaluate the model's prediction performance. At the same time, the sensitivity of the model to input uncertainty is studied through perturbation analysis to ensure that the model is sufficiently robust in practical applications. This systematic uncertainty modeling method provides a reliable probabilistic framework for spatial planning decisions, which can effectively support risk assessment and solution optimization.

[0158] S4.6: According to the joint probability constraint matrix, the uncertainty range of the Wasserstein distance evaluation distribution is introduced to construct a fuzzy set of probability distribution;

[0159] In step S4.6, Wasserstein distance is introduced to measure the uncertainty range of probability distribution. For any two probability distributions, Wasserstein distance describes the minimum "cost" required to transform one distribution into another. Based on this distance metric, a fuzzy set is constructed. The core area of ​​the set corresponds to a highly certain prediction result, and the edge area represents the uncertainty range of the prediction. For example, in the dynamic monitoring of wetlands, the Wasserstein distance is used to evaluate the fluctuation of the wetland range caused by changes in hydrological conditions, and then determine a reasonable uncertainty boundary.

[0160] S4.7: Based on the fuzzy set, according to the risk assessment result of real-time monitoring, dynamically adjust the constraint parameters of the basic constraint condition to obtain a probabilistic constraint;

[0161] In step S4.7, a dynamic adjustment mechanism for constraints is implemented. First, a real-time risk assessment module is established to calculate the reliability index of the prediction by comparing the deviation between the prediction results and the actual observed data. Then, based on these indicators and the characteristics of the fuzzy set, the parameters of the basic constraints are dynamically updated. Specifically, if the prediction error in a certain area continues to increase, the constraints in that area will be appropriately relaxed; conversely, if the prediction performance is stable, the constraints may be tightened to improve accuracy.

[0162] This dynamic adjustment process takes into account the impact of time and space scales. In terms of time scale, it distinguishes between short-term fluctuations and long-term trends, and adopts different constraint strategies for changes at different time scales; in terms of space scale, it considers local characteristics and regional associations to ensure the spatial consistency of constraints. For example, in urban expansion monitoring, short-term constraints are adjusted according to seasonal construction activities, while maintaining long-term constraints on the overall urban development direction.

[0163] The resulting probabilistic constraint is an adaptive constraint system that can be continuously optimized and adjusted based on real-time monitoring data and risk assessment results. This constraint not only takes into account the uncertainty characteristics of a single feature, but also includes the interaction between features and the dynamic characteristics of the system as a whole, providing a reliable constraint framework for subsequent change prediction and risk assessment.

[0164] S5: Based on the uncertainty distribution model of the probability constraint, the minimum safety distance is calculated in combination with the uncertainty ellipsoid, an adaptive safety distance rule is generated, and a collision avoidance rule is obtained through verification. According to the collision avoidance rule, a multi-objective optimization function of land use efficiency index and environmental protection index is established;

[0165] The uncertainty distribution model based on the probability constraint is combined with the uncertainty ellipsoid to calculate the minimum safety distance, generate an adaptive safety distance rule, and obtain a collision avoidance rule, including:

[0166] S5.1: Based on the probability constraints, construct a confidence ellipsoid around the object and establish a safety envelope of the object;

[0167] In step S5.1, a confidence ellipsoid is constructed to characterize the spatial uncertainty of the object. Based on the probability-constrained uncertainty distribution model, the covariance matrix of each object is first calculated, and then the main axis direction and length of the ellipsoid are determined by eigenvalue decomposition. The size of the ellipsoid is determined by a preset confidence level (usually 95% or 99%). The shape characteristics of the object are also considered, and the irregular object is approximated to an ellipsoid shape through the convex hull algorithm to ensure the integrity of the safety envelope.

[0168] S5.2: Calculate the minimum separation distance between adjacent ground object confidence ellipsoids and determine the initial safety spacing;

[0169] In step S5.2, an efficient algorithm is used to calculate the minimum separation distance between adjacent object confidence ellipsoids. In specific implementation, a fast screening algorithm is first used to identify potential interaction pairs, and then the minimum distance between ellipsoids is accurately calculated through an iterative optimization method. This calculation process takes into account the motion characteristics of the objects and environmental constraints.

[0170] S5.3: For the initial safety distance, a fuzzy controller is used to achieve smooth adjustment of the safety distance to obtain an adaptive safety distance rule;

[0171] In step S5.3, an adaptive adjustment mechanism based on fuzzy control is designed. The fuzzy controller receives multiple input variables, including relative speed, surrounding environment complexity, historical collision risk, etc., and dynamically adjusts the safety distance through fuzzy inference rules. The design of fuzzy rules adopts a method that combines expert knowledge and machine learning, such as "if the relative speed is large and the environment is complex, increase the safety distance." The parameters of the fuzzy rules are continuously optimized through online learning, so that the safety distance can smoothly adapt to the needs of different scenarios.

[0172] S5.4: For the adaptive safety spacing rule, generate multiple sets of scenarios by using the Monte Carlo method for verification to obtain a collision avoidance rule;

[0173] In step S5.4, the Monte Carlo method is used for large-scale verification. First, a scenario generator is built based on historical data and expert knowledge to randomly generate test scenarios containing different combinations of objects, motion patterns, and environmental conditions. Then, the adaptive safety spacing rule is applied in each scenario, and the collision risk and system performance indicators are statistically analyzed. Through a large number of simulation experiments, the limitations of the rules can be identified and targeted optimization can be performed.

[0174] According to the collision avoidance rule, a multi-objective optimization function of land use efficiency index and environmental protection index is established, including:

[0175] S5.5: Convert the collision avoidance rule into a spatial layout constraint;

[0176] In step S5.5, the collision avoidance rules are converted into specific spatial layout constraints. These constraints mainly include: minimum spacing constraints (based on safe spacing rules), spatial morphology constraints (based on geometric features of land objects), and functional compatibility constraints (based on land object usage). For example, for the layout between industrial land and residential land, the minimum buffer distance is set according to the environmental impact; for the planning of ecological corridors, the minimum width and connectivity requirements are set. These constraints are expressed in the form of mathematical inequalities to facilitate subsequent optimization solutions.

[0177] S5.6: Based on the spatial layout constraints, construct a land use efficiency objective function with floor area ratio and building density as parameters;

[0178] In step S5.6, the embodiment of the present application constructs an objective function for evaluating land use efficiency. This function mainly considers two key parameters: floor area ratio and building density. The floor area ratio reflects the vertical utilization efficiency of the land, and the building density represents the horizontal utilization intensity. The design of the objective function takes into account the differences in land use types. For example, commercial areas may pursue a higher floor area ratio, while residential areas need to balance development intensity and living quality. An accessibility index is also introduced to evaluate the spatial connection efficiency between different functional areas.

[0179] S5.7: According to the spatial layout constraints, set an environmental impact objective function that includes ecological connectivity and environmental carrying capacity;

[0180] In step S5.7, the embodiment of the present application designs an environmental impact objective function. This function contains two core indicators: ecological connectivity and environmental carrying capacity. Ecological connectivity is calculated by graph theory to evaluate the functional connection between ecological patches; environmental carrying capacity considers the resource capacity, pollutant diffusion and ecological sensitivity of the region. For example, in urban expansion planning, the barrier effect of new construction land on ecological corridors is evaluated, while considering the carrying capacity of regional water resources and atmospheric environment.

[0181] S5.8: Construct a multi-objective optimization model by combining the constraints and the objective function;

[0182] In step S5.8, the embodiment of the present application integrates the aforementioned constraints and objective functions to construct a multi-objective optimization model. The mathematical form of the model is:

[0183] Decision variables: spatial location and development intensity of various types of land use;

[0184] Constraints: spatial layout constraints, resource and environmental constraints;

[0185] Objective function: maximize land use efficiency and minimize environmental impact.

[0186] S6: converting the collision avoidance rule into a mathematical inequality as a constraint condition, using the NSGA-III algorithm to solve the multi-objective optimization problem, generating multiple non-dominated solutions, and executing a spatial planning strategy based on the multiple non-dominated solutions.

[0187] In step S6, the embodiment of the present application uses NSGA-III (non-dominated sorting genetic algorithm III) to solve the multi-objective optimization problem. This algorithm is particularly suitable for dealing with optimization problems with multiple objectives and complex constraints. The specific implementation process includes:

[0188] Initialization: Generate an initial population that satisfies the constraints, where each individual represents a possible spatial layout solution.

[0189] Evolutionary operation: Generate new solutions through crossover and mutation operations while ensuring that the new solutions satisfy the collision avoidance rules.

[0190] Non-dominated sorting: Sort the solutions based on the Pareto dominance relationship and retain the non-dominated solutions.

[0191] Reference point fitness: Use the reference point method to maintain the diversity of the population and ensure that a uniformly distributed set of solutions is obtained.

[0192] Selection operation: The next generation of individuals is selected by comprehensively considering the non-dominated level and the distance to the reference point.

[0193] After multiple iterations, the algorithm eventually generates a set of non-dominated solutions, each of which represents a feasible planning scheme. These schemes reflect different trade-offs between land use efficiency and environmental protection. For example, some schemes may prioritize the protection of ecological corridors but have a lower floor area ratio, while others may achieve a higher development intensity but require more environmental compensation measures.

[0194] like Figure 7 As shown, the present application also provides a national land space planning and mapping system based on dynamic remote sensing technology, including:

[0195] The acquisition module 10 is used to construct a multi-source remote sensing data acquisition network based on a linear parameter variation LPV model, wherein the LPV model includes a state matrix describing changes in ground objects and a time-varying parameter vector describing seasonal changes and climate conditions. Through the multi-source remote sensing data acquisition network, time-series remote sensing data including satellite remote sensing, UAV remote sensing and ground observation station remote sensing are acquired;

[0196] An extraction module 20 is used to perform spatiotemporal registration and noise filtering on the time series remote sensing data according to the state matrix and the time-varying parameter vector to obtain standardized dynamic remote sensing data, and extract spatiotemporal features using a deep neural network including a spatiotemporal convolution layer and an attention mechanism based on the standardized dynamic remote sensing data, including multispectral data, SAR data and DEM data;

[0197] An identification storage module 30 is used to identify key change areas through a dual attention module according to the extracted spatiotemporal features, store the spatiotemporal features and key change area information into a hierarchical dynamic feature library, and obtain dynamic change features of ground objects;

[0198] A construction module 40 is used to construct an uncertainty distribution model using a Gaussian mixture model and a Beta distribution model according to the dynamic change characteristics of the ground object, and to design a probability constraint according to the uncertainty distribution model;

[0199] A generation module 50 is used to calculate the minimum safe distance based on the uncertainty distribution model of the probability constraint in combination with the uncertainty ellipsoid, generate an adaptive safe distance rule, obtain a collision avoidance rule through verification, and establish a multi-objective optimization function of land use efficiency index and environmental protection index according to the collision avoidance rule;

[0200] The execution module 60 is also used to convert the collision avoidance rules into mathematical inequalities as constraints, use the NSGA-III algorithm to solve the multi-objective optimization problem, generate multiple non-dominated solutions, and execute the space planning strategy according to the multiple non-dominated solutions.

[0201] The basic principles of the present application are described above in conjunction with specific embodiments. However, it should be noted that the advantages, strengths, effects, etc. mentioned in the present application are only examples and not limitations, and it cannot be considered that these advantages, strengths, effects, etc. are required by each embodiment of the present application. In addition, the specific details disclosed above are only for the purpose of illustration and ease of understanding, not for limitation, and the above details do not limit the present application to being implemented by adopting the above specific details.

[0202] It should also be noted that in the system and method of the present application, each component or each step can be decomposed and / or recombined. Such decomposition and / or recombination should be regarded as equivalent solutions of the present application.

Claims

1. A national land space planning surveying and mapping method based on dynamic remote sensing technology, characterized in that: include: Construct a multi-source remote sensing data acquisition network based on a linear parameter variation LPV model, wherein the LPV model includes a state matrix describing changes in ground objects and a time-varying parameter vector describing seasonal changes and climate conditions, and acquires time-series remote sensing data including satellite remote sensing, UAV remote sensing and ground observation station remote sensing through the multi-source remote sensing data acquisition network; Performing spatiotemporal registration and noise filtering on the time series remote sensing data according to the state matrix and the time-varying parameter vector to obtain standardized dynamic remote sensing data, wherein the standardized dynamic remote sensing data includes multispectral data, SAR data and DEM data, and extracting spatiotemporal features using a deep neural network including a spatiotemporal convolution layer and an attention mechanism according to the standardized dynamic remote sensing data; According to the extracted spatiotemporal features, the key change areas are identified through a dual attention module, and the spatiotemporal features and the key change area information are stored in a hierarchical dynamic feature library to obtain the dynamic change features of the ground objects; According to the dynamic change characteristics of the ground object, a Gaussian mixture model and a Beta distribution model are used to construct an uncertainty distribution model, and probability constraints are designed according to the uncertainty distribution model; Based on the uncertainty distribution model of the probability constraint, the minimum safety distance is calculated in combination with the uncertainty ellipsoid, an adaptive safety distance rule is generated, and a collision avoidance rule is obtained through verification. According to the collision avoidance rule, a multi-objective optimization function of land use efficiency index and environmental protection index is established; The collision avoidance rules are converted into mathematical inequalities as constraints, and the NSGA-III algorithm is used to solve the multi-objective optimization problem, generate multiple non-dominated solutions, and execute a spatial planning strategy based on the multiple non-dominated solutions.

2. The method according to claim 1, characterized in that The method of performing spatiotemporal registration and noise filtering on the time series remote sensing data according to the state matrix and the time-varying parameter vector to obtain standardized dynamic remote sensing data includes: constructing a state space model according to the time series remote sensing data, wherein a system matrix of the state space model depends on a time-varying parameter vector; Based on the state space model, a space-time registration algorithm is used to calculate the transformation matrix between different data sources to obtain the registered remote sensing data; For the registered remote sensing data, using a Kalman filter to perform adaptive noise suppression to obtain filtered data; According to the filtered data, updating the state matrix describing the change of the ground object in the LPV model; The filtered data is standardized based on the updated state matrix to obtain standardized dynamic remote sensing data.

3. The method according to claim 1, characterized in that The method of extracting spatiotemporal features based on the standardized dynamic remote sensing data using a deep neural network including a spatiotemporal convolution layer and an attention mechanism includes: Receiving the standardized dynamic remote sensing data, performing normalization preprocessing on the multispectral data, the SAR data and the DEM data respectively through a multi-branch input layer to obtain preprocessed multi-source data; For the preprocessed multi-source data, a 3D convolution layer is constructed, the convolution kernel size is set to 3, and the parameters of the step size are set to 1 to obtain preliminary spatiotemporal features; For the preliminary spatiotemporal features, multi-scale features are obtained by adding a residual block including two convolutional layers and a batch normalization layer; For the multi-scale features, a temporal attention module is set to calculate the attention weights of features at different time points to obtain a weighted spatiotemporal feature map.

4. The method according to claim 3, characterized in that: The method of identifying key change areas through a dual attention module based on the extracted spatiotemporal features includes: For the weighted spatiotemporal feature map, a correlation matrix between positions of the spatiotemporal feature map is calculated, and a softmax function is used to generate spatial weights to obtain a spatial attention map; For the weighted spatiotemporal feature map, global average pooling and a fully connected layer are used to calculate the dependency between channels to obtain a channel attention map; Combining the spatial attention map and the channel attention map, weighting the weighted spatiotemporal feature map in both the spatial dimension and the channel dimension to obtain a key change area; Based on the key change area, the spatial attention map and the channel attention map are fused to obtain a change significance map.

5. The method according to claim 4, characterized in that The spatiotemporal features and key change area information are stored in a hierarchical dynamic feature library to obtain the dynamic change features of the ground objects, including: Build a hierarchical dynamic feature library including basic feature layer, change pattern layer and knowledge graph layer; Based on the weighted spatiotemporal feature map and the change significance map, storing the original feature vector in the basic feature layer and establishing a feature index; According to the data of the basic feature layer, a density clustering algorithm is used to identify the change pattern, and the change pattern of the ground object is stored in the change pattern layer; Based on the analysis results of the change pattern layer, storing information representing the causal relationship of the change in the knowledge graph layer; With respect to the multi-layer feature structure of the hierarchical dynamic feature library, the feature distribution is dynamically maintained through an incremental update mechanism based on similarity calculation to obtain the dynamic change features of the ground objects.

6. The method according to claim 5, characterized in that According to the dynamic change characteristics of the ground object, the Gaussian mixture model and the Beta distribution model are used to construct an uncertainty distribution model, including: For the dynamic change characteristics of the ground objects, the EM algorithm is used to estimate the parameters of the Gaussian mixture model describing the distribution of the ground object positions, so as to obtain a probability distribution model of the ground object positions; According to the probability distribution model, a multidimensional Beta distribution model describing the shape parameters of the land object is established, and a membership function is introduced to represent the boundary fuzziness, thereby forming an uncertainty description of the land object characteristics; Based on the uncertainty description of the features of the land object, a dynamic Bayesian network is constructed to obtain an uncertainty distribution model.

7. The method according to claim 6, characterized in that The designing of probability constraints according to the uncertainty distribution model comprises: Based on the prediction results of the uncertainty distribution model, a collision risk probability threshold ε is set for each ground object, and the value range of ε is 0.01 to 0.05, so as to establish basic constraint conditions; For the basic constraints, a Copula function is used to construct a joint probability constraint matrix among multiple ground objects; According to the joint probability constraint matrix, the uncertainty range of the Wasserstein distance evaluation distribution is introduced to construct a fuzzy set of probability distribution; Based on the fuzzy set and in accordance with the risk assessment result of real-time monitoring, the constraint parameters of the basic constraint conditions are dynamically adjusted to obtain a probabilistic constraint.

8. The method according to claim 7, characterized in that The uncertainty distribution model based on the probability constraint is combined with the uncertainty ellipsoid to calculate the minimum safety distance, generate an adaptive safety distance rule, and obtain a collision avoidance rule, including: Based on the probability constraints, construct a confidence ellipsoid around the object and establish a safety envelope of the object; Calculate the minimum separation distance between adjacent ground object confidence ellipsoids and determine the initial safety distance; For the initial safety distance, a fuzzy controller is used to achieve smooth adjustment of the safety distance, thereby obtaining an adaptive safety distance rule; The adaptive safety spacing rule is verified by generating multiple sets of scenarios using the Monte Carlo method to obtain a collision avoidance rule.

9. The method according to claim 1, characterized in that: The multi-objective optimization function of land use efficiency index and environmental protection index is established according to the collision avoidance rule, including: converting the collision avoidance rules into spatial layout constraints; Based on the spatial layout constraints, a land use efficiency objective function with volume ratio and building density as parameters is constructed; According to the spatial layout constraints, an environmental impact objective function including ecological connectivity and environmental carrying capacity is set; Combining the constraints and objective functions, a multi-objective optimization model is constructed.

10. A national land space planning and mapping system based on dynamic remote sensing technology, characterized in that: include: An acquisition module is used to construct a multi-source remote sensing data acquisition network based on a linear parameter variation LPV model, wherein the LPV model includes a state matrix describing changes in ground objects and a time-varying parameter vector describing seasonal changes and climate conditions. Through the multi-source remote sensing data acquisition network, time-series remote sensing data including satellite remote sensing, UAV remote sensing and ground observation station remote sensing are acquired; An extraction module, for performing spatiotemporal registration and noise filtering on the time series remote sensing data according to the state matrix and the time-varying parameter vector to obtain standardized dynamic remote sensing data, and extracting spatiotemporal features using a deep neural network including a spatiotemporal convolution layer and an attention mechanism based on the standardized dynamic remote sensing data, including multispectral data, SAR data and DEM data; An identification and storage module is used to identify key change areas through a dual attention module according to the extracted spatiotemporal features, store the spatiotemporal features and key change area information into a hierarchical dynamic feature library, and obtain dynamic change features of ground objects; A construction module is used to construct an uncertainty distribution model using a Gaussian mixture model and a Beta distribution model according to the dynamic change characteristics of the ground object, and to design a probability constraint according to the uncertainty distribution model; A generation module is used to calculate the minimum safe distance based on the uncertainty distribution model of the probability constraint in combination with the uncertainty ellipsoid, generate an adaptive safe distance rule, obtain a collision avoidance rule through verification, and establish a multi-objective optimization function of land use efficiency index and environmental protection index according to the collision avoidance rule; The execution module is used to convert the collision avoidance rules into mathematical inequalities as constraints, adopt the NSGA-III algorithm to solve the multi-objective optimization problem, generate multiple non-dominated solutions, and execute the space planning strategy according to the multiple non-dominated solutions.

Images