National space planning suitability evaluation system based on spatiotemporal big data

Abstract

The present application belongs to the technical field of space planning, and discloses a land space planning suitability evaluation system based on space-time big data; by constructing a projection zone distribution topology graph and an anisotropic distortion tensor field, coordinate distortion at the junction of projection zones is accurately identified and corrected. By using a spatial autocorrelation analysis method, systematic deviation areas in the residual system in the process of digitizing historical maps are identified, and a residual propagation influence domain model is established, so that accurate coordinate correction of each position point is realized. Through an adaptive fusion mechanism of compensation components, the fusion accuracy of multi-source space-time data is improved. A coordinate correction confidence evaluation mechanism is introduced, data uncertainty information is included in the suitability evaluation process, and more scientific and reliable evaluation results and confidence intervals are provided for decision makers.

Description

Technical Field

[0001] This invention relates to the field of spatial planning technology, and more specifically, to a land spatial planning suitability evaluation system based on spatiotemporal big data. Background Technology

[0002] Existing land spatial planning suitability evaluation systems face the problem of coordinate transformation accuracy loss when processing multi-source spatiotemporal big data, which restricts the accuracy and reliability of evaluation results. In practical applications, the transformation between different projection coordinate systems often produces significant accuracy attenuation in boundary areas. In particular, the deformation distortion of the UTM projection system in the edge areas can reach tens of centimeters or even meters, which causes a series of chain errors in high-precision urban planning and the delineation of ecologically sensitive areas. At the same time, the digitization of historical maps inevitably introduces control point residuals, which manifest as systematic offsets in fine-scale land suitability evaluations, leading to erroneous judgments in historical evolution analysis. In the Gauss-Kruger projection system widely used in land surveying, the problem of mixed use of 3-degree and 6-degree zone data is particularly prominent. In key areas such as inter-provincial boundaries, large watersheds, and ecological corridors, there is a significant coordinate jump phenomenon near the zone boundaries, resulting in inconsistent geometric representations of the same geographic entity in different data sources, and its deformation gradient exhibits complex anisotropic characteristics along the boundary line. This coordinate system difference is particularly pronounced in mountainous and hilly areas. When the terrain undulations intersect with the projection zone boundary, the elevation factor further amplifies the coordinate distortion effect, artificially disrupting the spatial continuity of key factors such as soil and water conservation and ecological suitability. Traditional methods often employ global transformation parameters or simple linear models for coordinate transformation, failing to effectively capture and correct these precision losses due to significant spatial heterogeneity. This results in uncertainty blind spots in suitability evaluation results in key areas of planning decisions, affecting the scientific nature and implementation effectiveness of national land spatial planning.

[0003] In view of this, the present invention proposes a land spatial planning suitability evaluation system based on spatiotemporal big data to solve the above problems. Summary of the Invention

[0004] To overcome the aforementioned deficiencies of the existing technology and to achieve the above objectives, the present invention provides the following technical solution: a land spatial planning suitability evaluation system based on spatiotemporal big data, comprising: The spatial data acquisition module is used to acquire multi-source spatial data of the evaluation area and the corresponding projection coordinate system metadata of each data source. The topology construction module is used to parse the projection zone parameters of each data source based on the projection coordinate system metadata, and construct the projection zone distribution topology map of the evaluation area. The sensitive region identification module is used to extract the boundary line of the projection zone based on the topological map of the projection zone distribution and calculate the scale factor gradient on both sides of the boundary line to identify the coordinate distortion sensitive region. The tensor field construction module is used to construct anisotropic distortion tensor fields in coordinate distortion sensitive regions. The anisotropic distortion tensor fields characterize the differences in coordinate distortion along different directions. The offset recognition module is used to obtain the labeled coordinates and measured coordinates of the digitized control points on the historical map, calculate the residual vector of each control point and perform spatial autocorrelation analysis to identify the systematic offset area of ​​the residual. The influence domain modeling module is used to construct an influence domain model for residual propagation based on the distribution characteristics of the residual systematic offset region. The compensation calculation module is used to calculate the comprehensive coordinate compensation of each location point in multi-source spatial data based on the anisotropic distortion tensor field and residual propagation influence domain model. The baseline layer generation module is used to correct the accuracy of multi-source spatial data based on the comprehensive coordinate compensation amount and generate a suitability evaluation factor layer under a unified spatial baseline. The suitability assessment module is used to calculate the land spatial planning suitability index of each spatial unit in the assessment area based on the suitability assessment factor layer and coordinate correction confidence. The modules are connected via wired and / or wireless means to enable data transmission between them.

[0005] The technical effects and advantages of the land spatial planning suitability evaluation system based on spatiotemporal big data of this invention are as follows: This invention improves the accuracy and reliability of suitability assessment for land spatial planning, effectively overcoming the problem of coordinate inconsistency in the process of multi-source spatiotemporal data fusion. By accurately processing coordinate distortion in the boundary areas of projection zones, it eliminates the accuracy cliff effect in boundary areas in traditional methods, achieving seamless splicing and fusion of spatial data. For systematic residuals generated by the digitization of historical maps, this invention can intelligently identify and correct their spatial propagation characteristics, improving the temporal continuity and reliability of land use change analysis. The anisotropy correction method makes the geographical representation of planning units more accurate, especially in complex areas such as mountainous areas, watersheds, and administrative boundary intersections, ensuring the spatial consistency of suitability assessment factors. The data quality assessment mechanism introduced in this invention presents uncertainty information intuitively to decision-makers, enhancing the risk control capabilities of planning schemes. Attached Figure Description

[0006] Figure 1 This is a schematic diagram of the land spatial planning suitability evaluation system based on spatiotemporal big data of the present invention. Detailed Implementation

[0007] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0008] This application provides a land spatial planning suitability evaluation system based on spatiotemporal big data. The executing entities of this system include, but are not limited to, those that run on the system: land spatial planning platforms, geographic information systems, spatial decision support systems, intelligent planning assistance systems, etc., which can be regarded as general computing nodes of this application. The suitability evaluation system includes, but is not limited to, at least one of the following: cloud spatial analysis engine, distributed geographic computing platform, multi-source data fusion system.

[0009] This invention provides a land spatial planning suitability evaluation system based on spatiotemporal big data. Through spatial data acquisition, analytical topology construction, sensitive area identification, tensor field construction, offset identification, influence domain modeling, compensation calculation, benchmark layer generation, and suitability evaluation, it achieves accurate assessment of the suitability of multi-source spatiotemporal data for land spatial planning after projection transformation and historical error correction. It features high accuracy and reliability, effectively resolving issues such as coordinate system differences between multi-source data, projection zone boundary deformation, and historical map digitization errors, providing scientific evaluation results under a unified benchmark.

[0010] Please see Figure 1 In this embodiment of the invention, the land spatial planning suitability evaluation system based on spatiotemporal big data includes: The spatial data acquisition module is used to acquire multi-source spatial data of the evaluation area and the corresponding projection coordinate system metadata of each data source. Multi-source spatial data includes various spatiotemporal information such as remote sensing imagery, vector maps, historical maps, statistical data, and field survey data, acquired in real time through a multi-channel data interface. The projection coordinate system metadata records key information such as the coordinate system type, central meridian parameter, projection method, and transformation parameters of the data source. Remote sensing imagery provides current land use information for the area; vector maps record basic geographic elements such as administrative boundaries and transportation facilities; historical maps record the historical evolution of the area; statistical data includes socioeconomic indicators such as population, economy, and resources; and field survey data provides the latest ground-based verification information. This data provides comprehensive raw material for subsequent evaluation, ensuring the completeness and accuracy of the evaluation.

[0011] All multi-source spatial data used in this system originates from publicly available and legal channels, including authorized public databases. The data collection and usage process strictly adheres to relevant laws and regulations, ensuring the legality, compliance, and privacy protection requirements of data acquisition, and does not involve any illegally collected or unauthorized data sources. During data processing, only necessary coordinate unification and accuracy correction are performed on publicly available data; the original public attributes and usage permissions of the data remain unchanged.

[0012] The topology construction module is used to analyze the projection zone parameters of each data source based on the projection coordinate system metadata, and construct a projection zone distribution topology map of the evaluation area. The projection zone distribution topology map includes the spatial distribution, boundary relationships, and data coverage of different projection zone types such as 3-degree zones and 6-degree zones. Projection zone parameters include key information such as the position of the central meridian, bandwidth type, and coordinate origin setting. By analyzing these parameters, the projection zone location and type of each data source within the evaluation area are identified. The topology map records the adjacency relationships and regional coverage between projection zones in a graph structure, providing a spatial reference framework for subsequent analysis.

[0013] The sensitive region identification module is used to extract the boundary lines of projection zones based on the topological map of the projection zone distribution and calculate the scale factor gradients on both sides of the boundary lines to identify coordinate distortion-sensitive regions. Coordinate distortion-sensitive regions refer to the areas of precision discontinuity at the boundary lines of projection zones caused by coordinate system transformations; these regions are a key challenge in multi-source data fusion. The module samples and detects the boundary lines, calculates the coordinate jump vectors and scale factor gradients between adjacent projection zones, and accurately locates sensitive segments with severe coordinate distortion, providing target regions for subsequent tensor field construction.

[0014] The tensor field construction module is used to construct an anisotropic distortion tensor field within the coordinate distortion-sensitive region. This anisotropic distortion tensor field characterizes the differences in coordinate distortion along different directions. The distortion tensor field is a high-order mathematical model describing spatial deformation characteristics, expressing the differences in coordinate deformation in different directions through a second-order tensor. The module first establishes a local coordinate system in the sensitive region, calculates the distortion rates in the tangential and normal directions, constructs the distortion tensor, and performs eigenvalue decomposition to obtain the principal distortion directions and principal values, forming a complete anisotropic distortion tensor field, providing a precise mathematical basis for coordinate compensation.

[0015] The offset identification module acquires the labeled and measured coordinates of digitized control points on historical maps, calculates the residual vectors for each control point, performs spatial autocorrelation analysis, and identifies regions of systematic residual offset. These systematic offset regions reflect systematic errors generated during the digitization of historical maps, which exhibit spatial clustering and directional consistency. By comparing the labeled and measured coordinates of control points, the module calculates the residual vectors and uses spatial autocorrelation analysis to identify clusters of residual points with significant spatial correlation, determining the spatial extent of the systematic offset regions and providing a data foundation for subsequent influence domain modeling.

[0016] The influence domain modeling module is used to construct an influence domain model for residual propagation based on the distribution characteristics of the systematic offset regions of the residuals. The residual propagation influence domain model describes the propagation and attenuation law of systematic errors in space; the influence gradually weakens with increasing distance from the propagation source point. The module first determines the geometric center of each offset region as the propagation source point, calculates its residual strength and attenuation coefficient, constructs a spatial attenuation function, and superimposes the influence of multiple propagation source points to generate a complete influence domain model, providing a theoretical basis for accurately calculating the residual compensation amount at each location.

[0017] The compensation calculation module is used to calculate the comprehensive coordinate compensation for each location point in multi-source spatial data based on the anisotropic distortion tensor field and residual propagation influence domain model. The comprehensive coordinate compensation integrates spatial error correction information from two dimensions: projection distortion and historical residuals. The module first calculates the distortion compensation component and the residual compensation component for the point to be compensated, then analyzes the synergy and interference relationship between the two compensation components, constructs adaptive fusion coefficients, and achieves the organic synthesis of the compensation components to generate the final comprehensive coordinate compensation, providing precise adjustment parameters for data accuracy correction.

[0018] The baseline layer generation module is used to correct the accuracy of multi-source spatial data based on comprehensive coordinate compensation, generating a suitability evaluation factor layer under a unified spatial baseline. This suitability evaluation factor layer is a high-precision spatial dataset in a unified coordinate system, eliminating the effects of projection distortion and historical errors. The module first corrects the coordinates of each data element, then confirms the correction effect through spatial consistency checks, selects a unified coordinate system for data transformation and resampling, and finally fuses the data based on spatial location confidence to generate a complete evaluation factor layer. It also records coordinate correction confidence information, providing data support for the final suitability evaluation.

[0019] The suitability assessment module calculates the land use suitability index for each spatial unit within the assessment area based on a suitability assessment factor layer and coordinate-corrected confidence levels. The land use suitability index is a comprehensive evaluation indicator reflecting the suitability of land use in a region. The module first extracts the assessment factor values ​​for each unit from the assessment factor layer, calculates reliability weights based on coordinate-corrected confidence levels, and labels the uncertainty of factor values ​​located in sensitive areas. Through weighted comprehensive calculation, it derives the suitability index and confidence interval, providing a scientific basis for land use planning decisions.

[0020] The modules are connected via wired and / or wireless means to enable data transmission between them.

[0021] In this embodiment of the invention, the detailed implementation steps for constructing the projection zone distribution topology map of the evaluation region include: The process involves analyzing the central meridian and projection band width of each data source to determine its projection band type, which includes 3-degree and 6-degree bands. The analysis first extracts projection parameter information from the data source's metadata, including key parameters such as the central meridian longitude, projection method, and bandwidth settings. For historical data lacking metadata, projection parameters are estimated using a control point inversion method; that is, feature points with known geographic coordinates are selected, and the most likely combination of projection parameters is calculated using coordinate inversion. The projection band type is determined using an empirical threshold method, based on bandwidth parameters or the interval between adjacent central meridians to determine whether it is a 3-degree or 6-degree band. A 3-degree band has a central meridian interval of 3 degrees and is generally used for large-scale mapping; a 6-degree band has an interval of 6 degrees and is commonly used for small- to medium-scale regional coverage. This step provides basic information about the data source's projection type for subsequent analysis and is a prerequisite for constructing topological relationships.

[0022] Based on the spatial coverage and projection band type of each data source, an attribution table is generated between the data source and the projection band. This attribution table is the core data structure linking the data source to the spatial reference system. The generation process first calculates the spatial extent (minimum bounding rectangle) of each data source, then performs spatial overlay analysis with the coverage of each projection band to calculate the overlay ratio and determine the primary and secondary coverage bands. For cross-band data, the primary band attribution is determined based on the coverage area ratio, while secondary band information is recorded. The attribution table is stored in an association matrix format, with rows representing data source IDs, columns representing projection band numbers, and cell values ​​representing coverage levels. This structured attribution table clearly defines the spatial reference framework for each data source, providing the basic mapping relationships for subsequent topology construction.

[0023] An adjacency graph of projection zones is constructed, using projection zones as nodes and the boundaries between adjacent projection zones as edges. The projection zone adjacency graph is a graph theory model that expresses the spatial relationships between projection zones, intuitively demonstrating the organizational structure of regional projection zones. The construction process employs graph theory methods, representing each projection zone as a node in the graph and assigning attributes such as zone number, central meridian, and bandwidth to each node. Then, based on the geographic adjacency relationships, connecting edges are established between adjacent projection zones, forming a complete adjacency graph structure. The weight of the edges can be set as the length of the boundary line, reflecting the importance of the boundary relationship. For special cases, such as 3-degree and 6-degree zones covering the same area, special types of associated edges are established and marked as type overlap relationships. This graph structure clearly expresses the spatial organization of projection zones, providing a topological basis for identifying sensitive boundary areas.

[0024] The coverage areas of each data source are marked in the adjacency graph of the projection zone, generating a topology map of the projection zone distribution. The topology map is an enhanced model that integrates data source distribution information on top of the adjacency graph, comprehensively displaying the spatial relationship between data and projection zones. The marking process first maps the data source information in the attribution table to the corresponding nodes in the adjacency graph, adding a data source coverage list attribute to the nodes; then, based on the actual coverage range of the data sources, precise coverage area boundaries are marked on the graph, forming an enhanced topology map containing data distribution information. For important boundary areas, detailed descriptions are added, recording the precise geometric shape and location characteristics of the boundary lines. The final generated topology map not only includes the adjacency relationships of projection zones but also integrates the spatial distribution of data sources, intuitively displaying the organization structure of the projection system and data within the evaluation area, providing a comprehensive spatial reference framework for subsequent sensitive area identification and zone conflict analysis.

[0025] Overlapping regions simultaneously covered by both 3-degree and 6-degree projection zone data are identified in the topology map of the projection zone distribution and marked as mixed-use conflict zones. These mixed-use conflict zones are a key challenge in multi-source data fusion, as data from different projection zone systems exhibit inconsistencies in coordinate system transformation within these areas. The identification process first extracts regions simultaneously covered by both 3-degree and 6-degree projection zone data from the topology map, forming preliminary candidate overlap areas. Then, spatial overlay analysis is used to accurately calculate the boundary range of the overlapping regions. Finally, based on the degree of overlap and data importance, the conflict level is assessed, and these areas are marked as mixed-use conflict zones of different priorities. For high-priority conflict zones, detailed information on the data sources within the region and a description of the conflict characteristics are additionally recorded to provide a basis for subsequent precise processing. This identification and classification of conflict zones provides a targeted strategy for the system to handle the mixing of data from different projection systems and is a crucial step in ensuring the effective fusion of multi-source data.

[0026] In this embodiment of the invention, the detailed implementation steps for identifying coordinate distortion-sensitive regions include: Detection point pairs are deployed along the boundary line of the projection zone at preset sampling intervals, with each pair of detection points located symmetrically on both sides of the boundary line. These detection point pairs are the basic sampling units for analyzing the distortion characteristics of the boundary line, revealing the discontinuity of the projection transformation by comparing the coordinate differences on both sides. The deployment process first extracts the precise geometry of the boundary line of the projection zone and parameterizes it as a continuous curve. Then, an adaptive sampling strategy is used to set the interval, increasing the sampling density in areas with large curvature changes and appropriately reducing the number of sampling points in areas with gentle changes, ensuring the representativeness and efficiency of the sampling. Finally, at each sampling location, symmetrical detection point pairs are deployed on both sides of the boundary line's normal direction, recording the geographical coordinates and the projection zone information of the point pairs. The sampling interval is usually dynamically determined based on the total length and complexity of the boundary line, generally between 100 meters and 1 kilometer, ensuring that local variation features are captured without incurring excessive computational burden. This structured detection point pair deployment method provides a systematic sampling framework for subsequent coordinate jump analysis.

[0027] For each pair of detection points, coordinate forward calculation is performed using the parameters of their respective projection zones to obtain the planar coordinates in the two projection zone coordinate systems. Coordinate forward calculation is the mathematical process of converting geographic coordinates (latitude and longitude) into projected planar coordinates, and it is a fundamental step in analyzing coordinate differences. The calculation process first determines the projection zone parameters of each detection point, including key information such as the central meridian, projection type, and zone number; then, the standard Gauss-Kruger projection algorithm is applied to convert the geographic coordinates of the detection points into planar coordinates (x, y) under the corresponding projection zone; finally, the coordinate results under the two projection zones are saved to form the basic data for coordinate comparison. The core calculation of Gauss-Kruger projection involves complex ellipsoidal mathematics, but modern GIS libraries have provided efficient implementations. The forward calculation results directly reflect the differences in the representation of the same geographical location by different projection zones, providing direct evidence for quantifying coordinate jump phenomena.

[0028] The difference in planar coordinates of the same geographical location under two projection zone coordinate systems is defined as the coordinate jump vector. The coordinate jump vector is a core indicator for measuring discontinuities at the boundary of projection zones, directly reflecting the coordinate changes during cross-zone transformation. The definition process involves performing a vector difference operation on the planar coordinates of the same detection point under two projection zones to obtain a two-dimensional jump vector representing direction and magnitude. The formula is as follows: ; in, For testing points The coordinate jump vector, and These are the planar coordinates of the point in the first and second projection zone coordinate systems, respectively.

[0029] The magnitude of the jump vector reflects the severity of the coordinate jump, while the direction reveals the dominant direction of the deformation. The jump vector is typically largest near the boundary and gradually decreases with distance from the boundary. This vectorized representation of the jump not only quantifies the degree of coordinate discontinuity but also preserves the directional information of the deformation, providing complete change characteristics for subsequent gradient analysis.

[0030] The rate of change of coordinate jump vectors between adjacent detection point pairs is calculated as the scale factor gradient. The scale factor gradient is a differential index describing the rate of change of coordinate deformation space, reflecting the local acceleration or deceleration trend of deformation. The calculation process first constructs a one-dimensional sequence of detection points along the boundary line direction to ensure the order of the points; then, the difference in jump vectors between adjacent detection point pairs is calculated and divided by the geographical distance between the point pairs to obtain the normalized rate of change, as shown in the formula: ; in, For the first and Detect the scale factor gradient between point pairs. For the first The coordinate jump vector of each detection point This represents the geographical distance between adjacent detection points.

[0031] Regions with high gradient values ​​indicate drastic spatial changes in coordinate deformation, and are typically sensitive locations requiring special attention. This differential gradient analysis method effectively identifies regions of accelerated deformation, providing quantitative evidence for accurately locating sensitive areas.

[0032] Sections with scale factor gradients exceeding a preset gradient threshold are marked as coordinate distortion sensitive regions. The final determination of sensitive regions is based on a gradient threshold screening process aimed at identifying the most significant deformation-prone key sections. This process first analyzes the gradient distribution characteristics of the entire boundary line, calculating statistical indicators such as the mean and standard deviation. Then, an appropriate threshold is set, typically the mean plus 2-3 times the standard deviation, or the quantile method is used to select the top 10%-20% of high gradient values. Finally, sections with gradients continuously exceeding the threshold are identified and marked as coordinate distortion sensitive regions, with their spatial range and severity recorded. For particularly complex regions, the threshold can be lowered or auxiliary judgment indicators added to ensure complete capture of sensitive regions. This statistical threshold-based sensitive region delineation method considers both the global deformation level and highlights local anomaly characteristics, providing clear target regions for subsequent tensor field construction and early warning information for key and challenging areas in multi-source data fusion.

[0033] In this embodiment of the invention, the detailed implementation steps for constructing the anisotropic distortion tensor field include: A local coordinate system is established within the region sensitive to coordinate distortion. The X-axis of the local coordinate system is tangential to the boundary line of the projection zone, and the Y-axis is normal to the boundary line. The local coordinate system serves as a reference frame describing the deformation around the boundary line, simplifying deformation analysis through appropriate axial setting. The establishment process first constructs a parametric curve along the boundary line, calculating the tangent direction at each point as the X-axis direction; then, the normal direction is determined as the Y-axis, forming an orthogonal two-dimensional local coordinate system; finally, a transformation relationship from global geographic coordinates to the local coordinate system is established, achieving seamless mapping between the coordinate systems. This local coordinate system, based on the boundary line, naturally decomposes deformation analysis into two main directions: along the boundary line and perpendicular to the boundary line. This aligns with the physical characteristics of projection deformation and provides an ideal reference frame for subsequent directional analysis.

[0034] Calculate the components of the coordinate jump vector along the X and Y axes in the local coordinate system. Calculating the direction components is a crucial step in decomposing the two-dimensional deformation into principal directions, facilitating the separation of distortion effects in different directions. The calculation process first transforms the jump vector of each detection point to the local coordinate system, and then obtains the scalar components along the X and Y axes through vector projection, using the following formula: ; ; in, and These are the components of the jump vector along the X-axis (tangential) and Y-axis (normal), respectively. and Let be the unit basis vectors of the local coordinate system. This represents the vector dot product operation.

[0035] This orthogonal decomposition method simplifies complex two-dimensional deformation into one-dimensional changes in two directions, preserving the complete information of the deformation while making the analysis more intuitive and targeted. The decomposed component data provides the basis for calculating the directional distortion rate.

[0036] The normal distortion rate is defined as the ratio of the Y-axis component to the distance from the boundary line of the projection zone. The normal distortion rate is a normalized index describing deformation perpendicular to the boundary line, reflecting the rate of change of distance. The calculation process first determines the precise distance from each detection point to the boundary line; then, it calculates the ratio of the Y-axis jump component to the distance to obtain the normalized distortion rate; finally, it interpolates the sampling points across the entire sensitive area to form a continuous normal distortion rate field. This index directly reflects the discontinuity of the coordinate system perpendicular to the boundary line, typically reaching its maximum near the boundary line and decaying with increasing distance. The spatial distribution of the normal distortion rate reveals the decay pattern of the boundary line's influence, providing important directional characteristics for subsequent tensor field construction.

[0037] The rate of change of the X-axis component along the boundary line is used as the tangential distortion rate. The tangential distortion rate characterizes the deformation gradient along the boundary line, reflecting the deformation differences between different segments of the boundary line. The calculation process first constructs a one-dimensional sequence of points along the boundary line to ensure spatial continuity; then, the difference in the X-axis component between adjacent points is calculated and divided by the distance along the boundary line to obtain the normalized rate of change; finally, interpolation along the line forms the complete tangential distortion rate distribution. Regions with high tangential distortion rates typically correspond to locations with complex boundary line shapes or significant topographic relief, requiring special attention. This gradient analysis method along the line reveals the spatial variability of the boundary effect, providing another key directional component for the tensor field.

[0038] A second-order distortion tensor is constructed based on the normal and tangential distortion rates. The second-order distortion tensor is a mathematical tool for comprehensively expressing directional deformation, capturing the interrelationships between deformations in different directions through matrix form. The construction process first uses the normal and tangential distortion rates as the main diagonal elements of the tensor; then, the off-diagonal elements are estimated based on the cross gradient to reflect the coupling effect of deformations in different directions; finally, the symmetry and positive definiteness of the tensor are ensured to meet the basic requirements of physical deformation. For each grid point within the sensitive region, a corresponding local tensor is constructed, forming a discrete tensor field sampling. This tensor expression based on physical deformation theory not only preserves directional information but also describes the complex coupling characteristics of deformation, making it a high-level mathematical model for accurately characterizing coordinate distortion.

[0039] Eigenvalue decomposition is performed on the second-order distortion tensor to obtain the principal distortion directions and principal distortion values, generating an anisotropic distortion tensor field. Eigenvalue decomposition is a key mathematical operation for extracting the essential features of the tensor, revealing the principal axis system and intensity of the deformation. The decomposition process uses a standard matrix eigenvalue algorithm to decompose the second-order tensor at each position into eigenvalues ​​and eigenvectors, as shown in the formula: ; in, It is a second-order distortion tensor. and These are eigenvalues ​​(distortion principal values). and This is the corresponding feature vector (the principal direction of distortion). and This is the corresponding transpose.

[0040] Eigenvalues ​​reflect the distortion intensity along the principal direction, while eigenvectors indicate the directions of maximum and minimum distortion. By integrating the eigenvalue decomposition results of all grid points, a complete anisotropic distortion tensor field is constructed, comprehensively describing the directional characteristics and spatial distribution of coordinate distortion within the sensitive region. This high-level mathematical expression based on tensors provides a theoretical foundation for subsequent precise coordinate correction, enabling targeted compensation for deformation effects in different directions and significantly improving the accuracy of multi-source data fusion.

[0041] In this embodiment of the invention, the detailed implementation steps for calculating the residual vector of each control point and performing spatial autocorrelation analysis include: The difference between the labeled coordinates and the measured coordinates of each control point is calculated to obtain the control point residual vector. The residual vector is a fundamental indicator for evaluating the accuracy of historical map digitization, directly reflecting the direction and magnitude of positional deviation. The calculation process first collects control point information from the historical map digitization process, including labeled coordinates (coordinates recorded during digitization) and measured coordinates (reference coordinates for high-precision measurements). Then, coordinate system unification processing is performed to ensure that the two sets of coordinates are within the same reference frame. Finally, the coordinate difference is calculated to obtain the residual vector representing the magnitude and direction of the deviation. Control points are typically selected from clearly identifiable and still existing geographical features on historical maps, such as the corners of important buildings, road intersections, or topographical features. The statistical distribution of the residual vector reflects the overall accuracy level of the digitization process, providing fundamental data for identifying systematic errors.

[0042] A spatial neighborhood of control points is constructed using a preset search radius, and the residual vectors of other control points within each control point's neighborhood are statistically analyzed. Spatial neighborhood analysis is a key method for discovering local patterns, identifying spatial correlations by comparing the residual characteristics of neighboring points. The analysis process first determines an appropriate search radius for each control point, typically based on the average density of control points and the spatial scale of the study area, ensuring that the neighborhood contains sufficient sample points without being excessively large. Then, a spatial index structure is constructed to quickly identify the neighborhood members of each point. Finally, the direction and magnitude characteristics of all residual vectors within the neighborhood are collected and statistically analyzed, providing a local sample set for subsequent correlation analysis. The search radius is set using an adaptive strategy, appropriately increasing in sparse control point areas and appropriately decreasing in dense areas to ensure statistical reliability and representativeness. This spatial neighborhood-based local analysis method effectively captures spatial variation patterns of residuals and is a fundamental step in identifying systematic shifts.

[0043] Calculate the ratio of the mean angle between the residual vector of each control point and the residual vectors of its neighboring control points to their magnitude. The mean angle and magnitude ratio are core indicators for quantifying the spatial correlation of residuals, reflecting the consistency of residuals in a local area. The calculation process first calculates the angle between the residual vector of each control point and the residual vectors of all points in its neighborhood, using the vector dot product formula: ; in, Control points and The angle between the residual vectors, and For the corresponding residual vector, Represents the vector dot product. This represents the magnitude of the vector.

[0044] Then, the mean of all included angles is calculated to obtain the orientation consistency index; simultaneously, the ratio of the residual vector magnitude to the average magnitude of the neighborhood is calculated to assess the relative consistency of magnitude. These two indices together reflect the degree of residual similarity between the control point and its surrounding environment, and are key evidence for identifying systematic patterns. Regions with high orientation consistency usually indicate the presence of systematic shifts rather than random errors, and require special attention.

[0045] A set of control points with a mean directional angle less than a preset angle threshold and a modulus ratio within a preset range is identified as a control point cluster with systematic offset characteristics. Control point cluster identification is a crucial step in extracting systematic patterns from discrete points, aiming to distinguish between random errors and systematic offsets. The identification process first sets appropriate thresholds, typically 30° to 45° for angles to reflect directional consistency requirements; the modulus ratio range is usually 0.5-2.0, allowing some variation but excluding extreme outliers. Then, control points meeting these criteria are selected to form a preliminary cluster set. Finally, connectivity analysis is performed to combine spatially connected points into complete control point clusters, while isolated points or small-scale point groups are removed to reduce noise impact. This cluster identification method based on multi-index thresholds effectively separates point sets with systematic characteristics, providing a reliable point set basis for determining the residual region.

[0046] Convex hull expansion is performed on control point clusters to generate the systematic offset region of the residuals. Convex hull expansion is a transformation process from discrete point clusters to a continuous region, aiming to determine the spatial extent of the systematic offset. The expansion process first calculates the standard convex hull for each control point cluster, obtaining the smallest convex polygon containing all points; then, based on the density and distribution characteristics of the points, the convex hull is appropriately buffered and expanded to ensure coverage of the potential influence area; finally, the expansion boundary is smoothed to eliminate jagged edges and sharp corners, forming a natural region boundary. For particularly complex point distributions, an alpha shape can be considered instead of the standard convex hull to better adapt to non-convex distribution patterns. This point-to-surface spatial expansion method transforms discrete control point observations into a continuous influence region, providing a clear spatial definition for subsequent residual influence domain modeling and a basis for regional division for error correction during the evaluation process.

[0047] In this embodiment of the invention, the detailed implementation steps for constructing the residual propagation influence domain model include: The geometric center of each systematic offset region in the residuals is designated as the propagation source point. The propagation source point is the starting location for the diffusion of residual influences and represents the core region of the systematic offset. The determination process first calculates the geometric center of each offset region, using either the centroid method or the weighted centroid method. Then, the representativeness of the center point is evaluated, ensuring it is located in a densely populated area of ​​control points rather than at the edges or in gaps. Finally, the center points that meet the conditions are set as propagation source points, and their precise location and region affiliation information are recorded. For offset regions with particularly irregular shapes, density-weighted centers are considered instead of simple geometric centers to better reflect the actual distribution of control points. As the theoretical origin of the residual influences, the propagation source point directly affects the accuracy of subsequent propagation model construction; its accurate location is a crucial guarantee for model effectiveness.

[0048] The average residual vector of the control point cluster corresponding to each propagation source point is calculated as the residual intensity of that source point. Residual intensity is a core indicator of the influence of the propagation source point, reflecting the strength and direction of the systematic migration. The calculation process first collects all residual vectors within the control point cluster to which the propagation source point belongs; then, these vectors are statistically analyzed to calculate the average direction and magnitude; finally, the average result is set as the residual intensity vector of the source point. For control points with unevenly distributed heights, a weighted averaging strategy is considered, assigning different weights according to the reliability or importance of the points. The residual intensity vector contains both directional information, indicating the dominant direction of the migration, and magnitude information, reflecting the significance of the migration. This intensity definition based on actual observations ensures that the initial conditions of the propagation model conform to the actual residual distribution characteristics, providing a reliable intensity benchmark for subsequent spatial influence calculations.

[0049] The residual propagation attenuation coefficient is calculated based on the residual intensity and the control point density in the region where the propagation source point is located. The attenuation coefficient is a key parameter describing the degree to which the residual influence weakens with distance, determining the spatial extent of the influence domain. The calculation process first analyzes the control point distribution characteristics of each source point region, including density, evenness, and spatial autocorrelation; then, based on these characteristics and the magnitude of the residual intensity, the most suitable attenuation coefficient is estimated using empirical formulas or machine learning methods; finally, the coefficient is normalized to ensure comparability between different regions. The physical meaning of the attenuation coefficient is the magnitude of spatial damping during the propagation of the residual influence, reflecting the degree to which the geographical environment within the region suppresses error propagation. High-density areas typically experience faster attenuation, while open areas experience slower attenuation. This data-characteristic-based attenuation coefficient calculation method enables the propagation model to adaptively adjust to the characteristics of different regions, improving the model's applicability and accuracy.

[0050] Using the residual intensity as the initial value and the residual propagation attenuation coefficient as the attenuation rate, a spatial attenuation function for the residual influence at each propagation source point is constructed. The spatial attenuation function is a mathematical model describing how the residual influence changes with distance and is a core component of the influence domain calculation. The construction process employs an improved exponential attenuation model, considering directionality and regional characteristics; the formula is: ; in, The distance from the source point is , direction angle is The residual influence vector at the location, The residual intensity vector of the source point. The basic attenuation coefficient, It is the directional coefficient. The residual principal direction angle is denoted as .

[0051] This model not only considers the fundamental law of distance decay but also introduces a directional adjustment term, allowing the influence to propagate further along the principal direction of the residuals and restricting its propagation in the vertical direction, thus conforming to the observed residual distribution characteristics. For each propagation source point, a dedicated decay function is constructed to reflect its unique influence propagation mode, providing a functional basis for subsequent influence superposition.

[0052] A residual propagation influence domain model is generated by superimposing the spatial decay functions of the residual influence from multiple propagation source points. This influence domain model is a comprehensive expression integrating the influences of multiple source points, describing the overall residual influence at each location within the study area. The superposition process employs a vector weighted summation method, calculating the contribution of influence from all source points for each spatial location and adaptively adjusting the weights based on the consistency of the influences. When the influence directions of multiple source points are consistent, the effects are enhanced by superposition; when the directions are opposite, the effects cancel each other out. The final generated influence domain model is a continuous vector field describing the distribution trend and intensity variation of the residuals throughout the region. This influence domain model, based on physical propagation theory, reasonably reproduces the spatial distribution pattern of actual residuals, providing a theoretical foundation for subsequent coordinate compensation calculations and a scientific tool for understanding the spatial patterns of historical map digitization errors.

[0053] In this embodiment of the invention, the detailed implementation steps for calculating the comprehensive coordinate compensation of each location point in multi-source spatial data include: Obtain the geographic coordinates of the point to be compensated, and calculate the normalized distance from this point to the boundary of the coordinate distortion-sensitive area and the spatial distance to each residual propagation source point. Calculating location parameters is a prerequisite step in determining the compensation amount, assessing the point's impact sensitivity through spatial relationships. The calculation process first determines the precise geographic location of the point to be compensated; then, it calculates its distance to the nearest sensitive area boundary and normalizes it using the characteristic scale of the sensitive area to ensure comparability between different areas; simultaneously, it calculates the Euclidean distance from the point to each residual propagation source point, serving as the basic parameter for residual impact assessment. Distance calculations use either spherical distance or projected plane distance, selecting the appropriate method based on the scale of the study area. Normalization makes the distance parameters more universal, facilitating standardized calculations for subsequent impact assessments. These spatial location parameters provide a fundamental basis for determining the point's impact status and are key inputs for subsequent accurate compensation calculations.

[0054] Based on the normalized distance, the local distortion tensor corresponding to the point to be compensated is obtained by interpolation from the anisotropic distortion tensor field. The local distortion tensor is a mathematical expression describing the projection deformation characteristics at the point location, obtained from the discrete tensor field through spatial interpolation. The interpolation process employs a tensor field-specific interpolation algorithm to ensure that the physical properties of the tensor are preserved during the interpolation process. For points close to the sensitive area, local polynomial interpolation is used to accurately capture the changing gradients near the boundary; for points far from the area, distance-weighted interpolation is used for a smooth transition. The interpolation result is a 2×2 symmetric matrix describing the directional deformation characteristics at the point location, containing the principal distortion value and principal direction information. This tensor interpolation method based on physical properties ensures the continuity and physical rationality of the compensation calculation, providing local deformation characteristics for accurate correction.

[0055] The distortion compensation component is obtained by performing a tensor product operation between the position vector of the point to be compensated in the original coordinate system and the local distortion tensor. Tensor product operation is a key step in transforming theoretical deformation into practical compensation, calculating the coordinate correction amount through mathematical transformation. The calculation process first constructs the point's position vector, typically using a reference point at the boundary of the sensitive area as the origin; then, the tensor product operation is applied, multiplying the position vector by the local distortion tensor to obtain the deformation effect vector; finally, based on the actual position of the point and the characteristics of the sensitive area, the compensation direction and scale are adjusted to generate the final distortion compensation component. The tensor product operation preserves the directional characteristics of the deformation, enabling precise correction of deformation differences in different directions. This compensation calculation method based on tensor theory achieves high-precision correction of anisotropic deformation, significantly improving the accuracy of cross-projection zonation data fusion.

[0056] Based on the spatial distance between the location to be compensated and each residual propagation source point, the influence weight of each propagation source point on that location is calculated using the residual propagation influence domain model. The influence weight is a quantitative indicator for assessing the importance of residual sources, determining the contribution ratio of each source point in the overall compensation. The calculation process is based on the previously constructed residual propagation influence domain model, considering distance attenuation and directional characteristics, to calculate the influence intensity of each source point on the target location; then these intensities are normalized and converted into weight coefficients, ensuring that the sum of the weights of all source points is 1. The weight calculation considers both spatial distance factors and the residual intensity and propagation characteristics of the source points, comprehensively reflecting the spatial differentiation law of residual influence. This weight calculation method based on the propagation model enables residual compensation to accurately reflect the comprehensive influence of multiple source points, providing a scientific basis for reasonable error correction.

[0057] The residual intensity vectors of each propagation source point are weighted and summed according to their influence weights to obtain the residual compensation component. The residual compensation component is a key vector for correcting errors in historical map digitization, reflecting the comprehensive effect of multi-source influences through weighted combination. The calculation process first obtains the residual intensity vector of each propagation source point; then, based on the influence weights calculated in the previous step, these vectors are weighted and summed to obtain a composite vector representing direction and magnitude; finally, necessary scale adjustments are made according to the location characteristics of the points and the data source to generate the final residual compensation component. The weighted summation process considers both the direction of the vectors, preserving the directional characteristics of the residuals, and the magnitude of the intensity, reflecting the significance of the influence. This vector synthesis-based compensation calculation method can reasonably reproduce complex spatial distributions of residuals, providing an effective tool for improving the accuracy of historical data.

[0058] The magnitudes of the vector dot product and cross product of the distortion compensation component and the residual compensation component are calculated. The dot product determines the degree of synergistic enhancement between the two compensation components, and the cross product magnitude determines the degree of orthogonal interference. Vector relationship analysis is a crucial step in evaluating the interaction between the two compensation mechanisms, quantifying synergistic and interference effects through mathematical operations. The calculation process first determines the two compensation component vectors (distortion and residual); then, their dot product is calculated to measure directional consistency and enhancement effect; simultaneously, the magnitude of the cross product is calculated to assess the interference degree of the perpendicular component. The formulas for calculating the dot product and cross product are: ; ; in, The dot product value, The cross product modulus value. and These are the distortion and residual compensation components, respectively.

[0059] A large and positive dot product value indicates that the two compensations are highly consistent and mutually reinforcing; a dot product close to zero indicates that they are independent; a negative dot product indicates that they are in opposite directions and cancel each other out. A large cross product magnitude indicates that the vertical component is significant and there is strong orthogonal interference. This relational analysis method based on vector algebra reveals the interaction patterns of different compensation mechanisms and provides a mathematical basis for reasonable integration.

[0060] An adaptive fusion coefficient for the compensation components is constructed based on the ratio of the degree of synergistic enhancement to the degree of orthogonal interference. The fusion coefficient is a key parameter for adjusting the contributions of different compensation components, and its optimal combination is achieved through adaptive calculation. The construction process first calculates the ratio of the degree of synergistic enhancement to the degree of orthogonal interference, which reflects the compatibility of the compensation components. Then, a nonlinear mapping function is designed based on this ratio to convert the original ratio into a fusion coefficient between 0 and 1. Finally, the coefficient is locally adjusted according to the location characteristics and data attributes of the points to ensure the continuity and rationality of the compensation results. The fusion coefficient design considers the physical meaning and interaction characteristics of the compensation components, favoring the enhancement combination effect when the components are consistent and selecting the dominant component when components conflict, thus avoiding erroneous superposition. This adaptive fusion method based on vector relationships achieves intelligent integration of multiple error sources, improving the accuracy and robustness of the compensation.

[0061] The distortion compensation component and the residual compensation component are vector-synthesized using adaptive fusion coefficients to generate a comprehensive coordinate compensation. Vector synthesis is a crucial step in generating the final compensation result, integrating correction effects from different sources through mathematical operations. The synthesis process employs a weighted vector combination method, adjusting the relative contributions of the two compensation components according to the fusion coefficients. The formula is as follows: ; in, For the comprehensive coordinate compensation amount, The fusion coefficient is... As an orthogonal adjustment parameter, it is usually set to a small value (0.1-0.3) or dynamically determined according to the regional characteristics.

[0062] The first two terms of the formula achieve the main weighted combination, while the third term is a special adjustment for the case of strong orthogonal interference, introducing an appropriate vertical component to balance the interactive effect of the two compensation mechanisms. This advanced vector synthesis method not only considers the linear combination of components but also incorporates a mechanism for handling nonlinear interactive effects, achieving precise correction of complex error structures and providing a mathematical foundation for high-precision fusion of multi-source spatiotemporal data.

[0063] In this embodiment of the invention, the detailed implementation steps for generating a suitability evaluation factor layer under a unified spatial benchmark include: The original coordinates of each element in multi-source spatial data are vector-superimposed with the comprehensive coordinate compensation values ​​of their corresponding locations to obtain the accurate element coordinates. Vector superposition is the core operation for coordinate correction, improving positional accuracy through mathematical calculations. The correction process first determines the original coordinates of all vertices or raster center points of the data elements; then, it calculates the corresponding comprehensive coordinate compensation value for each point; finally, it vector-sums the original coordinates with the compensation values ​​to obtain the corrected coordinates. For vector data, the topological relationships of the elements are maintained, and only the coordinate positions are adjusted; for raster data, a reverse mapping method is used for resampling to ensure data integrity. This point-level compensation-based correction method achieves precise adjustment of the data geometry while preserving the attribute information and internal structure of the original data, ensuring data integrity and consistency.

[0064] Spatial consistency check points are selected within the overlapping area of ​​the coordinate distortion-sensitive region and the residual systematic offset region. The residual coordinate bias of different data sources at these check points is calculated after correction. Spatial consistency checking is a crucial step in evaluating the correction effect, using residual bias to measure the accuracy of data fusion. The checking process first selects representative check points in key and challenging areas, typically choosing corresponding feature points from multi-source data that clearly correspond to the same ground features. Then, the coordinates of each data source at these points after correction are obtained. Finally, the coordinate differences between different data sources are calculated as a residual bias indicator. The point selection strategy focuses on the overlapping parts of the sensitive and offset regions, which are the most difficult to correct and require the highest accuracy. Residual bias statistics provide a quantitative assessment of the correction effect, serving as an important basis for judging whether the algorithm has achieved its expected goals and providing a clear indicator for iterative optimization.

[0065] The process involves determining whether the residual coordinate deviation exceeds a preset convergence threshold. If it does, the comprehensive coordinate compensation is iteratively corrected until the residual coordinate deviation converges below the preset threshold. Iterative correction is an adaptive process to improve accuracy, achieving high-precision fusion through iterative optimization. The judgment process first analyzes the statistical distribution of the residual deviation, calculating key indicators such as the mean, maximum, and standard deviation. Then, it compares this with a preset convergence threshold, typically set as the minimum accuracy requirement for planning. If the deviation exceeds the threshold, the iterative correction phase begins. The iterative process employs a gradient descent strategy, adjusting the comprehensive compensation based on the direction and magnitude of the residual deviation. The residual deviation is reassessed after each iteration until the convergence condition is met. This adaptive iterative method can handle complex nonlinear error structures, achieving high-precision matching through gradual approximation. It is an effective means of handling difficult areas, ensuring the quality and reliability of data fusion.

[0066] Based on the topological map of the projection zone distribution, the central meridian and projection parameters of the target unified coordinate system are determined, and the accuracy-corrected data sources are transformed to the target unified coordinate system. Unified coordinate system transformation is a crucial step in achieving consistent data representation, establishing a common spatial reference through projection transformation. The process first analyzes the geographical features and projection zone distribution of the evaluation area, selecting the most suitable central meridian and projection parameters, typically using the parameters of the projection zone where the main body of the area is located or specially designed optimized parameters. Then, a coordinate transformation model is constructed to transform the corrected data sources from the original coordinate system to the target unified system. Finally, a transformation quality assessment is conducted to ensure that the transformation process does not introduce new errors. The selection of the target coordinate system considers the region's geographical location, shape characteristics, and the distribution of major data sources, aiming to minimize overall deformation and data distortion. This unified coordinate selection and transformation method based on regional characteristics provides a consistent geometric basis for subsequent spatial analysis and evaluation, and is a necessary step in multi-source data fusion.

[0067] Spatial resampling is performed on the transformed data sources to unify raster resolution and vector feature node density. Spatial resampling is a technical process to achieve data structure unification, adjusting the spatial precision of the data through sampling transformation. The resampling process employs specific strategies for different data types: for raster data, the target resolution is first determined, typically selecting the optimal scale that preserves key information, and then resampling transformation is performed using methods such as bilinear or cubic convolution; for vector data, an appropriate node density standard is set, and simplification or encryption is performed using algorithms such as Douglas-Peucker to adjust to a uniform complexity level. The setting of resampling parameters considers a balance between the original data precision, planned application requirements, and computational efficiency, ensuring the preservation of key information while avoiding excessive computational burden. This unified data structure processing provides a consistent foundation of data for subsequent fusion analysis, improving computational efficiency and result reliability.

[0068] Based on the correction magnitude of each data source within the coordinate distortion-sensitive region and the residual systematic offset region, the spatial location reliability of each data source is calculated. Spatial location reliability is a key indicator for evaluating data quality, quantifying the reliability of the data through correction details. The calculation process first analyzes the distribution of correction vectors for each data source within the sensitive and offset regions, including characteristics such as directional consistency, amplitude uniformity, and spatial continuity. Then, a reliability assessment model is constructed based on the correction characteristics to quantify the location accuracy of the data sources in different regions. Finally, a reliability distribution map covering the entire evaluation area is generated, reflecting spatial heterogeneity. The reliability calculation considers the quality of the original data, the difficulty of correction, and the stability of the correction results, providing an objective basis for the weight allocation of data from different sources. This reliability assessment method based on correction characteristics achieves an explicit spatial expression of data quality, providing a scientific foundation for subsequent weighted fusion and also providing important input for the uncertainty analysis of the final results.

[0069] Weighted fusion of resampled data sources based on spatial location credibility generates a suitability assessment factor layer. Weighted fusion is a key technology for generating high-quality assessment factors, leveraging credibility weighting to highlight the contribution of high-quality data. The fusion process first determines the data source combination for each assessment factor and selects an appropriate fusion strategy based on the factor type. Then, a weighting scheme is designed based on spatial location credibility to ensure that high-credibility data dominates the results. Finally, spatial analysis algorithms such as weighted averaging, priority coverage, or Dempster-Shafer evidence theory are applied to generate a comprehensive assessment factor layer. The selection of fusion algorithms considers data characteristics and factor requirements; weighted averaging is preferred for continuous factors, while probability-based or evidence-based fusion methods are used for discrete factors. This intelligent fusion technology, which considers data quality, fully utilizes the complementary advantages of multi-source data, significantly improving the accuracy and reliability of assessment factors and providing a high-quality data foundation for the final suitability assessment.

[0070] The residual coordinate deviation and spatial location confidence of each spatial unit are recorded as the coordinate correction confidence of that unit. The coordinate correction confidence is a comprehensive indicator of data processing quality, quantifying the reliability of the results through multi-dimensional evaluation. The recording process first divides the evaluation area into a spatial unit grid of appropriate scale; then, it calculates the residual deviation statistics and average location confidence for each unit; finally, these two indicators are converted into a unified confidence index through weighted combination or function mapping and stored as the unit's metadata attribute. The confidence index intuitively reflects the quality and reliability of data fusion; high confidence areas indicate guaranteed data accuracy, while low confidence areas suggest potentially significant uncertainty. This explicit quality marking mechanism provides an important reference for uncertainty control in subsequent suitability assessments, enabling the evaluation results to consider spatial differences in data quality and enhancing the scientific rigor and credibility of the evaluation.

[0071] In this embodiment of the invention, the method for identifying band-type mixed conflict zones includes: This method extracts geographical features with the same name from areas simultaneously covered by both 3-degree and 6-degree zone data. These features are key comparative objects for identifying coordinate system conflicts, revealing differences by comparing geometric characteristics under different systems. The extraction process first identifies the overlapping area between the 3-degree and 6-degree zone data, which is a potential conflict zone. Then, it searches within the overlapping area for geographical features that clearly correspond to each other in both types of data, such as prominent buildings, road intersections, or natural landform features. Finally, it establishes a one-to-one correspondence between features, forming a comparison dataset. Feature matching typically employs a combination of spatial relationships and attributes to ensure accuracy. For historical or heterogeneous data, expert assistance may be needed to confirm complex correspondences. This comparison method based on features with the same name, through concrete comparison at the entity level, can intuitively reveal the actual impact of coordinate system differences, providing reliable basic data for conflict zone identification.

[0072] Geometric feature parameters of corresponding ground features were calculated in both 3-degree and 6-degree coordinate systems. These parameters included area, perimeter, and centroid coordinates. Geometric feature parameters are key indicators for quantitatively assessing the impact of coordinate system differences, comprehensively reflecting deformation effects through multi-dimensional comparison. The calculation process first ensured that the representation of corresponding features was consistent in both coordinate systems, guaranteeing the effectiveness of parameter comparison. Then, standard GIS spatial analysis functions were applied to accurately calculate the geometric parameter set of each feature in both coordinate systems. Finally, a parameter comparison table was generated, clearly recording the parameter differences for each pair of corresponding features. Area and perimeter reflect scale variation effects, while centroid coordinates reveal the overall displacement trend. These parameters together constitute a comprehensive set of geometric deformation indicators. The calculation method emphasizes precision control, employing high-precision algorithms for complex features to ensure reliable results. This multi-parameter geometric feature analysis method deeply reveals the actual impact of different coordinate systems on ground feature representation, providing a quantitative basis for conflict severity assessment.

[0073] Calculate the relative deviation of geometric feature parameters of the same geographic feature in two coordinate systems. Relative deviation is a standardized indicator that quantifies the influence of the coordinate system, expressed as a percentage for easy comparison and threshold judgment. The calculation process uses a standardized formula to convert absolute differences into relative change ratios. The formula is: ; in, For parameters The percentage of relative deviation and These are the parameter values ​​of the feature in the 3-degree zone and the 6-degree zone coordinate system, respectively.

[0074] The aforementioned relative ratio formula is used for area and perimeter, while for centroid coordinates, the ratio of displacement distance to feature scale is calculated. Relative deviation eliminates the influence of differences in feature size, allowing features of different types and sizes to be compared under a unified standard, facilitating the setting of universal thresholds and the identification of significant deviations. This standardized analysis method based on relative change objectively reflects the actual impact of coordinate system differences, providing a fair evaluation standard for identifying conflicting features.

[0075] Land features with relative deviations exceeding a preset deviation ratio are marked as conflict features. Conflict feature identification involves screening significantly affected objects from all features with the same name, using threshold judgments to determine targets requiring special attention. The marking process first sets appropriate deviation thresholds based on planning accuracy requirements and data quality standards; typically, relative deviation thresholds for area and perimeter are set at 2%-5%, and for centroid displacement at 1%-3% of the feature's characteristic scale. Then, the relative deviation of each feature is compared to the thresholds to identify those exceeding the standards. Finally, spatial analysis is performed on the identification results to assess the distribution patterns and clustering characteristics of conflict features. The threshold settings consider the accuracy requirements of land spatial planning and the technical difficulties of multi-source data fusion, ensuring the reliability of the identification while meeting planning needs. This threshold-based conflict identification method enables quantitative assessment and precise location of the impact of mixed-use zones, providing a core basis for subsequent conflict zone delineation.

[0076] The boundaries of mixed-use conflict zones are determined by the spatial distribution of conflict elements. Boundary determination is a spatial generalization process from discrete conflict elements to continuous conflict areas, defining the problem area through appropriate spatial analysis methods. The process first analyzes the spatial distribution characteristics of conflict elements, including density distribution, directionality, and clustering. Then, based on these characteristics, a suitable spatial expansion method is selected: convex hull or alpha shape methods are used for areas with significant clustering, while density thresholding or kernel density estimation methods are used for dispersed areas. Finally, the preliminary boundaries are smoothed and optimized to obtain a natural and coherent conflict zone boundary. Boundary determination not only considers the locational distribution of conflict elements but also integrates background information such as topography and administrative divisions, ensuring the practicality and operability of the delineated area. This point-to-area spatial integration method transforms discrete observations into a systematic regional division, providing a clear spatial definition for the management and handling of mixed-use conflict zones, and also providing regional guidance for data selection and processing strategies in evaluation work.

[0077] In this embodiment of the invention, the detailed implementation steps for calculating the suitability index of land spatial planning for each spatial unit in the evaluation area include: Evaluation factor values ​​for each spatial unit are extracted from the suitability evaluation factor layer. Evaluation factor extraction is a data preparation step in suitability analysis, obtaining characteristic parameters for each unit through spatial queries. The extraction process first determines the basic spatial unit for evaluation, which can be a regular grid or functional units such as administrative villages / communities. Then, factor values ​​for each unit are extracted from the previously generated evaluation factor layer through spatial overlay analysis. Finally, missing or outlier values ​​are processed to ensure data integrity and rationality. For vector factors, the statistical characteristics of the unit are calculated using spatial relationships; for raster factors, representative values ​​for the unit are obtained through regional statistics. Data extraction emphasizes spatial precision matching to ensure that factor values ​​accurately correspond to spatial units. This systematic factor extraction method provides structured input data for subsequent suitability calculations and is a fundamental step in achieving accurate evaluation.

[0078] The reliability weights of evaluation factor values ​​are calculated based on the coordinate-corrected confidence scores of each spatial unit. Reliability weights are a key mechanism for incorporating data quality into the evaluation process, reflecting the impact of uncertainty through weight adjustments. The calculation process first obtains the coordinate-corrected confidence score for each spatial unit, an indicator that integrates residual bias and location confidence information. Then, using an appropriate mapping function, the confidence scores are converted into standardized weight coefficients to ensure a reasonable weight distribution. Finally, the influence of the weights is adjusted according to the factor type and importance, forming the final reliability weight set. The mapping function typically uses an S-curve or piecewise linear function to ensure good discrimination in the medium-to-high confidence range. This data quality-based weighting mechanism enables the evaluation process to adaptively adjust factor contributions, reducing the negative impact of low-quality data and improving the reliability and robustness of the evaluation results. It is an effective method for handling data uncertainty.

[0079] Uncertainty labeling is applied to evaluation factor values ​​located in coordinate distortion-sensitive regions or residual systematic offset regions for each spatial unit. Uncertainty labeling serves as a warning mechanism to highlight risk areas, clearly indicating potential problems through explicit marking. The labeling process first identifies the overlap between spatial units and sensitive and offset regions to determine the degree of impact. Then, based on the overlap and severity of the impact, appropriate uncertainty labels are added to the relevant factor values, including information such as error range, confidence interval, or quality level. Finally, the labeling results are integrated into the factor data structure as auxiliary information in the evaluation. The labeling strategy is designed differently based on regional characteristics and factor types, employing more detailed uncertainty descriptions for key factors and high-risk areas. This explicit quality labeling method makes the evaluation process transparent in handling data uncertainty, providing risk warnings for decision-makers and background explanations for interpreting evaluation results, thus enhancing the scientific rigor and accountability of the evaluation work.

[0080] Based on reliability weights and uncertainty labeling, the evaluation factor values ​​are weighted and synthesized. Weighted synthesis is the core calculation process for generating the suitability index, integrating the influence of multiple factors through mathematical operations. The synthesis process first determines the evaluation model, selecting an appropriate mathematical framework based on planning objectives and regional characteristics, such as weighted linear combination, analytic hierarchy process (AHP), or fuzzy comprehensive evaluation. Then, based on the previously calculated reliability weights, the contribution coefficients of each factor are designed. Finally, the selected model is applied, combined with uncertainty information, to calculate and obtain preliminary comprehensive evaluation results. Model selection and parameter setting fully consider planning needs and regional characteristics, ensuring scientific rationality through a combination of expert knowledge and data-driven approaches. This weighted synthesis method, which considers data quality, achieves the organic integration of multiple factors while controlling the transmission effect of uncertainty, providing a reasonable numerical basis for the final suitability index, and is a core technical link in the evaluation process.

[0081] The suitability index for each spatial unit is calculated based on the weighted composite results, and the confidence interval of the suitability index is output. The suitability index is the final output of the evaluation, quantifying the suitability of a region through comprehensive indicators and providing a scientific basis for planning decisions. The calculation process is based on the previous weighted composite results, applying standardization and grading techniques to convert the original evaluation values ​​into a standardized suitability index for easier understanding and comparison. Simultaneously, considering uncertainty information, the confidence interval of the suitability index is calculated using Monte Carlo simulation or interval analysis methods, reflecting the reliable range of the results. Standardization employs quantile methods or natural breakpoint methods, and the grades are typically divided into three levels (high, medium, and low suitability) or a more detailed five-level division. The width of the confidence interval directly reflects the degree of certainty of the evaluation results, providing an important reference for decision-making risk assessment. This combination of point values ​​and intervals comprehensively demonstrates the results and reliability of the suitability evaluation, providing clear decision-making recommendations while retaining uncertainty information, thus providing comprehensive technical support for scientific decision-making in territorial spatial planning.

[0082] This invention achieves high-precision suitability evaluation of land spatial planning based on spatiotemporal big data through spatial data acquisition, analytical topology construction, sensitive area identification, tensor field construction, offset identification, influence domain modeling, compensation calculation, benchmark layer generation, and suitability evaluation. The tensor field analysis and residual propagation model of this invention solve the problems of coordinate system differences and historical errors in multi-source data, improving the accuracy and reliability of suitability evaluation.

[0083] The above are merely preferred embodiments of the present invention and are not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

[0084] It should be noted that all formulas in this manual are calculated by removing dimensions and taking their numerical values. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters and thresholds in the formulas are set by those skilled in the art according to the actual situation.

[0085] Although embodiments of the invention have been shown and described, those skilled in the art will understand that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the claims and their equivalents.

Claims

1. A land spatial planning suitability evaluation system based on spatiotemporal big data, characterized in that, include: The spatial data acquisition module is used to acquire multi-source spatial data of the evaluation area and the corresponding projection coordinate system metadata of each data source. The topology construction module is used to parse the projection zone parameters of each data source based on the projection coordinate system metadata, and construct the projection zone distribution topology map of the evaluation area. The sensitive region identification module is used to extract the boundary line of the projection zone and calculate the scale factor gradient on both sides of the boundary line based on the topology map of the projection zone distribution, and identify the coordinate distortion sensitive region. The tensor field construction module is used to construct an anisotropic distortion tensor field within the coordinate distortion sensitive region. The offset recognition module is used to obtain the labeled coordinates and measured coordinates of the digitized control points on the historical map, calculate the residual vector of each control point and perform spatial autocorrelation analysis to identify the systematic offset area of ​​the residual. The influence domain modeling module is used to construct a residual propagation influence domain model based on the distribution characteristics of the residual systematic offset region. The compensation calculation module is used to calculate the comprehensive coordinate compensation of each location point in the multi-source spatial data based on the anisotropic distortion tensor field and the residual propagation influence domain model. The baseline layer generation module is used to perform accuracy correction on the multi-source spatial data based on the comprehensive coordinate compensation amount, and generate a suitability evaluation factor layer under a unified spatial baseline. The suitability assessment module is used to calculate the land and space planning suitability index of each spatial unit in the assessment area based on the suitability assessment factor layer and the coordinate correction confidence level.

2. The system according to claim 1, characterized in that, The constructed evaluation region projection zone distribution topology map includes: The central meridian and projection zone width of each data source are analyzed to determine the projection zone type of each data source. The projection zone types include 3-degree zone type and 6-degree zone type. Based on the spatial coverage and projection band type of each data source, generate a table showing the attribution relationship between data sources and projection bands; Construct an adjacency graph of projection zones using projection zones as nodes and the boundary relationships between adjacent projection zones as edges; The coverage areas of each data source are marked in the adjacency map of the projection zone, and the topology map of the projection zone distribution is generated. Identify overlapping areas in the projection zone distribution topology map that are simultaneously covered by 3-degree zone data and 6-degree zone data, and mark them as zone-type mixed conflict areas.

3. The system according to claim 1, characterized in that, The region sensitive to coordinate distortion includes: Detection point pairs are arranged along the boundary line of the projection zone at a preset sampling interval, with each pair of detection points located symmetrically on both sides of the boundary line. For each pair of detection points, coordinate forward calculation is performed using the parameters of the respective projection zone to obtain the planar coordinates in the coordinate systems of the two projection zones; The difference in planar coordinates of the same geographical location in two projection zone coordinate systems is defined as the coordinate jump vector. Calculate the rate of change of the coordinate jump vector between adjacent detection point pairs, and use it as the scale factor gradient; The segment where the scale factor gradient is greater than a preset gradient threshold is marked as the coordinate distortion sensitive region.

4. The system according to claim 1, characterized in that, The construction of the anisotropic distortion tensor field includes: A local coordinate system is established within the coordinate distortion sensitive area, wherein the X-axis of the local coordinate system is tangential to the boundary line of the projection zone, and the Y-axis is normal to the boundary line of the projection zone. Calculate the components of the coordinate jump vector in the X and Y axes of the local coordinate system, respectively; The ratio of the Y-axis component to the distance from the boundary of the projection zone is used as the normal distortion rate; The rate of change of the X-axis component along the boundary line is taken as the tangential distortion rate. Based on the normal distortion rate and the tangential distortion rate, a second-order distortion tensor is constructed; The second-order distortion tensor is decomposed into eigenvalues ​​to obtain the principal distortion directions and principal distortion values, thereby generating the anisotropic distortion tensor field.

5. The system according to claim 1, characterized in that, The calculation of the residual vector for each control point and the performance of spatial autocorrelation analysis include: Calculate the difference between the labeled coordinates and the measured coordinates of each control point to obtain the control point residual vector; Construct a spatial neighborhood of the control point with a preset search radius, and count the residual vector of other control points in the neighborhood of each control point. Calculate the ratio of the mean angle between the direction of the residual vector of each control point and the residual vector of its neighboring control points to the magnitude; The set of control points whose mean directional angle is less than a preset angle threshold and whose modulus ratio is within a preset ratio range is identified as a cluster of control points with systematic offset characteristics. The control point cluster is expanded by convex hull to generate the residual systematic offset region.

6. The system according to claim 1, characterized in that, The construction of the residual propagation influence domain model includes: The geometric center of each of the systematic offset regions of the residual system is taken as the propagation source point; Calculate the average residual vector of the control point cluster corresponding to each propagation source point, and use it as the residual intensity of that propagation source point; The residual propagation attenuation coefficient is calculated based on the residual strength and the control point density in the region where the propagation source is located. Using the residual intensity as the initial value and the residual propagation attenuation coefficient as the attenuation rate, a spatial attenuation function of the residual influence of each propagation source point is constructed. The residual propagation influence domain model is generated by superimposing the residual influence space decay functions of multiple propagation source points.

7. The system according to claim 6, characterized in that, The calculation of the comprehensive coordinate compensation amount for each location point in the multi-source spatial data includes: Obtain the geographic coordinates of the location point to be compensated, and calculate the normalized distance of the location point from the boundary of the coordinate distortion sensitive area and the spatial distance from each residual propagation source point. Based on the normalized distance, the local distortion tensor corresponding to the location point to be compensated is obtained by interpolation from the anisotropic distortion tensor field; The distortion compensation component is obtained by performing a tensor product operation between the position vector of the point to be compensated in the original coordinate system and the local distortion tensor. Based on the spatial distance between the location to be compensated and each residual propagation source point, the influence weight of each propagation source point on the location is calculated using the residual propagation influence domain model. The residual intensity vectors of each propagation source point are weighted and summed according to the influence weights to obtain the residual compensation components. Calculate the magnitude of the vector dot product and the vector cross product of the distortion compensation component and the residual compensation component. Determine the degree of synergistic enhancement of the two compensation components based on the vector dot product, and determine the degree of orthogonal interference of the two compensation components based on the magnitude of the vector cross product. Based on the ratio of the degree of synergistic enhancement to the degree of orthogonal interference, an adaptive fusion coefficient for the compensation component is constructed; The distortion compensation component and the residual compensation component are vector-synthesized using the adaptive fusion coefficient of the compensation component to generate the comprehensive coordinate compensation amount.

8. The system according to claim 1, characterized in that, The generation of the suitability evaluation factor layer under a unified spatial benchmark includes: The original coordinates of each element in the multi-source spatial data are vector-superimposed with the comprehensive coordinate compensation amount of the corresponding location point to obtain the element coordinates after accuracy correction. In the overlapping area of ​​the coordinate distortion sensitive area and the residual systematic offset area, a spatial consistency check point is selected, and the coordinate residual deviation of different data sources at the check point is calculated after correction. Determine whether the coordinate residual deviation is greater than a preset convergence threshold. If it is, iteratively correct the comprehensive coordinate compensation amount until the coordinate residual deviation converges to below the preset convergence threshold. Based on the topology map of the projection zone distribution, determine the central meridian and projection parameters of the target unified coordinate system, and convert each data source after accuracy correction to the target unified coordinate system; Spatial resampling is performed on the converted data sources to unify the raster resolution and vector feature node density; The spatial location reliability of each data source is calculated based on the correction magnitude of each data source in the coordinate distortion sensitive area and the residual systematic offset area. Based on the spatial location credibility, the resampled data sources are weighted and fused to generate the suitability evaluation factor layer; The residual coordinate deviation and spatial position confidence of each spatial unit are recorded as the coordinate correction confidence of that spatial unit.

9. The system according to claim 2, characterized in that, The identification method for the strip-type mixed conflict zone includes: Extract the same-named ground features within areas simultaneously covered by both 3-degree and 6-degree zone data; Calculate the geometric feature parameters of the corresponding land feature in the 3-degree zone coordinate system and the 6-degree zone coordinate system respectively. The geometric feature parameters include area, perimeter and centroid coordinates. Calculate the relative deviation of geometric characteristic parameters of the same ground feature in two coordinate systems; Land features whose relative deviation is greater than a preset deviation ratio are marked as conflicting features; The boundary of the strip-type mixed-use conflict zone is determined by the spatial distribution range of the conflict elements.

10. The system according to claim 1, characterized in that, The calculation of the land spatial planning suitability index for each spatial unit in the evaluation area includes: Extract the evaluation factor values ​​of each spatial unit from the suitability evaluation factor layer; The confidence level is adjusted based on the coordinates of each spatial unit, and the reliability weight of the evaluation factor value is calculated. Uncertainty labeling is applied to the evaluation factor values ​​of each spatial unit located in the coordinate distortion sensitive region or the residual systematic offset region; The evaluation factor values ​​are weighted and synthesized based on the reliability weight and the uncertainty label. The suitability index for land use planning of each spatial unit is calculated based on the weighted composite results, and the confidence interval of the suitability index is output.

Images