Urban land utilization simulation method and system based on digital twinning

Abstract

The invention relates to the technical field of land resource management, in particular to an urban land utilization simulation method and system based on digital twinning, and the method comprises the steps: recognizing commercial nodes based on a power load, generating a resident viscosity index in combination with terminal signaling, extracting a non-commercial grid center to construct an expansion vector, and fusing physical blocking to correct boundary displacement. And reconstructing a commercial overflow area and rendering a morphological evolution result. According to the method, the building energy consumption characteristics and the mobile terminal signaling data are fused, the actual function attributes and the crowd gathering intensity of the urban space are mapped in real time, the resident viscosity index is quantified to accurately identify the potential high-activity area around the commercial node, and the discrete crowd heat is converted into the geometric expansion guide vector; the specific direction and amplitude of outward permeation of the commercial function are determined, and a physical blocking detection mechanism is used for automatically correcting and simulating displacement parameters for water flow and road isolation belts, so that land form dynamic reconstruction driven by multi-dimensional data is realized.

Description

Technical Field

[0001] This invention relates to the field of land resource management technology, and in particular to a method and system for simulating urban land use based on digital twins. Background Technology

[0002] The field of land resource management technology involves the entire process of management activities, including investigation and evaluation, planning and utilization, protection and remediation, and rights and registration supervision of land as a natural and economic complex. This field mainly uses surveying and mapping science, geographic information systems, and remote sensing technology to systematically monitor and analyze the quantity, quality, and spatial distribution of land resources in a region.

[0003] Among them, the urban land use simulation method refers to the technical means of using computer software to predict the changing trends of various types of land use within a city. It usually uses the land use status map of historical years as a basic reference, and selects the distance of the plot from the main road or the city center, the distribution density of surrounding infrastructure, and the flatness of the land as the influencing factors. By calculating the probability of each plot changing its use, and combining the rules of the driving or restricting effect of the use attributes of adjacent plots on the central plot, the land use type of each plot at a future time point is determined one by one, thereby generating a pre-show map reflecting the evolution of urban spatial form.

[0004] Existing technologies rely on historical static land use maps and physical attributes such as road distances as the basis for extrapolation, ignoring the dynamic changes under high-frequency urban operation. This results in a time lag between the simulation basis and the actual operating conditions, failing to reflect the functional evolution trend driven by real-time population behavior and energy consumption patterns. Consequently, the change probability calculated based on physical proximity deviates from the actual economic activity trend, making it difficult to capture the commercial expansion path driven by population gathering. This results in the pre-simulation map being merely a mechanical extrapolation of historical forms, lacking consideration of the interaction between micro-level vitality flow and physical environmental constraints, thus reducing the accuracy of spatial form evolution prediction and its practical reference value. Summary of the Invention

[0005] The purpose of this invention is to address the shortcomings of existing technologies by proposing a digital twin-based urban land use simulation method and system.

[0006] To achieve the above objectives, the present invention adopts the following technical solution: a method for simulating urban land use based on digital twins, comprising the following steps: S1: Obtain the instantaneous power load sequence of the building, construct the load feature vector and calculate the difference with the commercial and residential load benchmarks, and mark the commercial operation nodes; S2: Collect terminal signaling data of the commercial operation node, calculate the movement rate and generate dwell points based on the timestamp and coordinates of the terminal signaling data, and calculate the dwell stickiness index by counting the number and duration of dwell points. S3: Identify the non-commercial grid surrounding the commercial operation node, extract the geometric center of the non-commercial grid with the residency stickiness index, and construct an expansion guide vector based on the geometric center and the commercial operation node; S4: Calculate the simulated displacement parameters of the boundary of the commercial operation node based on the expansion guide vector and the residence viscosity index; detect physical blocking information based on the water flow boundary and road median layer; correct the simulated displacement parameters; and generate a corrected displacement vector. S5: Reconstruct the boundary of the commercial operation node based on the modified displacement vector, generate a commercial spillover area covering the non-commercial grid, and render the simulation results of the commercial land morphology evolution.

[0007] The present invention is improved in that the commercial operation node includes a unique index identifier of the target building and a commercial function attribute label determined based on the load characteristic difference degree; the dwell stickiness index is specifically calculated as the regional population attraction based on the sum of the distribution number of dwelling points and the sum of dwell time within the spatial grid; the expansion guidance vector includes the vector starting coordinates anchored to the center of the commercial operation node and the guiding orientation attribute pointing to the geometric center of the non-commercial spatial grid; the corrected displacement vector includes the displacement direction that avoids the obstruction of water flow boundaries and road medians and the displacement modulus after truncation correction based on physical obstruction information; and the commercial land morphology evolution simulation results include the reconstructed vector polygon vertex coordinate sequence, the updated commercial spillover zone geometric layer, and the rendered dynamic evolution visual view of the land boundary.

[0008] The present invention is improved in that step S1 is specifically as follows: S101: Obtain urban geographic vector data and building energy operation logs of the target area, construct a digital twin model of urban land use based on the urban geographic vector data and map energy data, extract the instantaneous power load sequence of buildings in the target area within a preset monitoring period in the model, perform extreme value difference calculation and fluctuation frequency statistics on the instantaneous power load sequence, and generate a load feature vector. S102: Map the load feature vector to a multi-dimensional feature measurement space, call the preset commercial operation load benchmark and residential load benchmark as spatial reference points, calculate the Euclidean distance between the load feature vector and the commercial operation load benchmark and the residential load benchmark in the measurement space, and generate a load attribute difference measurement value. S103: Based on the load attribute difference measurement value, set the filtering conditions, traverse and retrieve building units that match the characteristics of commercial electricity consumption patterns from the buildings in the target area, extract the unique identifier and spatial location information of the target building unit in the urban land use digital twin model, and mark the attribute status of the target building unit as commercial use to generate a commercial operation node.

[0009] The present invention is improved in that step S2 is specifically as follows: S201: Collect mobile terminal location signaling data covering the area where the commercial operation node is located, extract the timestamp sequence and latitude and longitude coordinate sequence from the signaling data stream, calculate the latitude and longitude coordinate displacement distance corresponding to adjacent timestamps to obtain the movement rate, filter latitude and longitude coordinates whose movement rate exceeds a preset speed threshold, and generate a terminal dwell point sequence. S202: Establish a spatial grid based on the commercial operation nodes, map the terminal dwell point sequence to the corresponding spatial grid, traverse and count the cumulative distribution number of dwell points in each spatial grid, calculate the cumulative dwell time of dwell points based on the timestamp difference, associate the cumulative distribution number and cumulative dwell time with the corresponding spatial grid ID, and generate a grid dwell distribution statistics set. S203: Call the grid dwell distribution statistics set, perform weighted normalization calculation on the cumulative distribution quantity and cumulative dwell time, quantify the degree of crowd gathering and dwell stickiness in each spatial grid, and assign a corresponding dwell stickiness index to each spatial grid.

[0010] The present invention is improved in that step S3 is specifically as follows: S301: Identify non-commercial spatial grids located adjacent to the geometric boundary of the commercial operation node from the spatial grid, obtain the residence stickiness index corresponding to the non-commercial spatial grid, filter spatial grids whose residence stickiness index exceeds the preset commercial spillover activation threshold and calculate the geometric center coordinates of the spatial grids, and generate the center coordinates of the spillover activation grids. S302: Spatially pair the geometric center coordinates of the commercial operation node with the center coordinates of the spillover activation grid, establish an azimuth connection line from the commercial center to the center of the spillover area, calculate the azimuth angle and distance data between the two points, and determine the expansion trend guiding coordinate pair; S303: Based on the spatial direction and distance defined by the expansion trend guide coordinate pair, analyze the direction of commercial land penetrating into the surrounding non-commercial areas, and use it as the expansion guide vector.

[0011] The present invention is improved in that step S4 is specifically as follows: S401: Based on the expansion guide vector and the residence viscosity index, calculate the outward extension distance of the boundary of the commercial operation node, generate outward radiating virtual rays in combination with the initial position of the boundary, define the initial trajectory and amplitude of the boundary expansion, and generate simulated displacement parameters. S402: Control the virtual ray corresponding to the simulated displacement parameter to extend in the digital twin model of urban land use, detect whether the virtual ray intersects with the water flow boundary layer or the road isolation strip layer in space, if the intersection occurs, extract the coordinates of the intersection point of the virtual ray with the boundary of each type of layer, mark the physical restriction position of the blockage expansion, and generate physical blockage information; S403: Based on the intersection coordinates of the physical blocking information, the simulated displacement parameters are truncated to restrict the virtual ray from crossing the physical isolation zone. The displacement endpoint in the blocked direction is adjusted to the intersection position, the original displacement parameters in the unblocked direction are retained, the magnitude and shape of the boundary expansion are corrected, and a corrected displacement vector is generated.

[0012] The present invention is improved in that the process of calculating the outward extrapolation distance of the boundary of the business operation node based on the expansion guidance vector and the residence stickiness index is specifically as follows: Extract the geometric modulus of the expansion guide vector as the basic expansion step size, and at the same time obtain the preset global maximum viscosity reference value; Dividing the residence viscosity index by the global maximum viscosity reference value yields the dimensionless weighting coefficient characterizing the expansion intensity. The outward extrapolation distance is obtained by using a weighted product of the basic expansion step size with dimensionless weight coefficients. The emission orientation of the virtual ray is determined based on the direction cosine of the expansion guiding vector, and the termination point of the virtual ray is set by combining the outward extrapolation distance, thus obtaining simulated displacement parameters including the starting point coordinates, the ending point coordinates, and the emission orientation.

[0013] The present invention is improved in that step S5 is specifically as follows: S501: Extract the vector boundary coordinates of the commercial operation node in the digital twin model of urban land use, apply the data of the corrected displacement vector to the corresponding boundary node, calculate the new spatial coordinate position of the boundary node after correction, connect the new coordinate points in sequence to close and form the updated polygon boundary, and generate a reconstructed boundary node sequence. S502: Obtain the geographic spatial range enclosed by the reconstructed boundary node sequence, incorporate the non-commercial spatial grid covered by the geographic spatial range into commercial land, change the corresponding land use attribute fields, and establish the geometric shape and area of ​​the newly added commercial land to generate a commercial spillover zone. S503: Overlay the commercial spillover area onto the map layer of the digital twin visual interface, update the coverage of the rendered commercial land, and generate simulation results of the evolution of commercial land morphology.

[0014] The present invention is improved in that the process of applying the corrected displacement vector data to the corresponding boundary nodes and calculating the new spatial coordinate positions of the boundary nodes after correction is specifically as follows: Iterate through and extract the x and y coordinates of the vector boundary coordinates, and analyze and correct the horizontal and vertical components of the displacement vector in the plane coordinate system. Perform a linear superposition operation of the horizontal coordinate value and the horizontal component, and a linear superposition operation of the vertical coordinate value and the vertical component to obtain the target displacement coordinates corresponding to the boundary node. The target displacement coordinates are reorganized according to the original connection order of the boundary nodes, and polygon topology rule detection is performed to remove abnormal coordinate points that cause the edge lines to self-intersect, thus establishing the new spatial coordinate position.

[0015] A digital twin-based urban land use simulation system, the system comprising: The load characteristic analysis module obtains the instantaneous power load sequence of the building, constructs the load characteristic vector, calculates the difference with the commercial and residential load benchmarks, and marks the commercial operation nodes; The dwell stickiness calculation module collects terminal signaling data from the commercial operation nodes, calculates the movement rate and generates dwell points based on the timestamps and coordinates of the terminal signaling data, and calculates the dwell stickiness index by counting the number and duration of dwell points. The expansion guide vector generation module identifies the non-commercial grid surrounding the commercial operation node, extracts the geometric center of the non-commercial grid with the residency stickiness index, and constructs the expansion guide vector based on the geometric center and the commercial operation node; The boundary displacement correction module calculates the simulated displacement parameters of the boundary of the commercial operation node based on the expansion guide vector and the residence viscosity index, and corrects the simulated displacement parameters based on the physical blocking information detected by the water flow boundary and the road median layer, thereby generating a corrected displacement vector. The commercial spillover zone reconstruction module reconstructs the boundaries of commercial operation nodes based on the corrected displacement vector, generates a commercial spillover zone covering the non-commercial grid, and renders the simulation results of commercial land morphology evolution.

[0016] Compared with the prior art, the advantages and positive effects of the present invention are as follows: In this invention, by integrating building energy consumption characteristics and mobile terminal signaling data, the actual functional attributes of urban space and the intensity of population gathering are mapped in real time. The dwell stickiness index is quantified to accurately identify potentially high-activity areas around commercial nodes. The heat of discrete populations is transformed into a geometric expansion guide vector to determine the specific direction and magnitude of commercial functions penetrating outward. The physical blockage detection mechanism is used to automatically correct the simulated displacement parameters for water flow and road median strips, realizing the dynamic reconstruction of land use morphology driven by multi-dimensional data. This ensures that the simulation results conform to the organic growth logic of urban economy and physical environment, and significantly improves the timeliness of land use evolution prediction and the rationality of spatial layout. Attached Figure Description

[0017] Figure 1 This is a flowchart of the method of the present invention; Figure 2 This is a detailed flowchart of step S1 of the present invention; Figure 3 This is a detailed flowchart of step S2 of the present invention; Figure 4This is a detailed flowchart of step S3 of the present invention; Figure 5 This is a detailed flowchart of step S4 of the present invention; Figure 6 This is a detailed flowchart of step S5 of the present invention; Figure 7 This is a system module diagram of the present invention. Detailed Implementation

[0018] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0019] Please see Figure 1 This invention provides a technical solution: a method for simulating urban land use based on digital twins, comprising the following steps: S1: Obtain the instantaneous power load sequence of the building, construct the load feature vector and calculate the difference with the commercial and residential load benchmarks, and mark the commercial operation nodes; S2: Collect terminal signaling data from commercial operation nodes, calculate the movement rate and generate dwell points based on the timestamp and coordinates of the terminal signaling data, and calculate the dwell stickiness index by counting the number and duration of dwell points. S3: Identify the non-commercial grid surrounding the commercial operation node, extract the geometric center of the non-commercial grid with the residency stickiness index, and construct an expansion guidance vector based on the geometric center and the commercial operation node; S4: Calculate the simulated displacement parameters of the commercial operation node boundary based on the expansion guide vector and the residence viscosity index, detect physical blocking information based on the water flow boundary and road median layer, correct the simulated displacement parameters, and generate the corrected displacement vector; S5: Reconstruct the boundaries of commercial operation nodes based on the corrected displacement vector, generate a commercial spillover zone covering the non-commercial grid, and render the simulation results of commercial land morphology evolution.

[0020] Commercial operation nodes include unique index identifiers for target buildings and commercial function attribute labels determined based on load characteristic differences. The dwell stickiness index is specifically calculated based on the total number of dwell points and the total dwell time within the spatial grid, resulting in regional population attractiveness. The expansion guidance vector includes the vector starting coordinates anchored at the center of the commercial operation node and the guiding orientation attribute pointing to the geometric center of the non-commercial spatial grid. The corrected displacement vector includes the displacement direction that avoids the obstruction of water flow boundaries and road medians and the displacement modulus corrected based on physical obstruction information. The simulation results of commercial land morphology evolution include the reconstructed vector polygon vertex coordinate sequence, the updated geometric layer of the commercial spillover area, and the rendered dynamic evolution visual view of the land boundary.

[0021] Please see Figure 2 Step S1 is as follows: S101: Obtain urban geographic vector data and building energy operation logs of the target area, construct a digital twin model of urban land use based on the urban geographic vector data and map energy data, extract the instantaneous power load sequence of buildings in the target area within the preset monitoring period in the model, perform extreme value difference calculation and fluctuation frequency statistics on the instantaneous power load sequence, and generate load feature vector; The system retrieves urban geographic vector data for the target area, including building outline coordinates, floor heights, and functional zoning codes. Simultaneously, it retrieves energy operation logs for all buildings within the corresponding area from the power supply management database, recording active power readings with a 15-minute time resolution. A geographic coordinate matching algorithm is used to associate the building outlines in the urban geographic vector data with the metering point numbers in the energy operation logs, thereby achieving spatial mapping of energy data within the urban land use digital twin model. A preset monitoring period of 30 natural days is selected, and the instantaneous power load sequence for each building in the target area within this period is extracted. For each instantaneous power load sequence, all values ​​in the sequence are traversed, the maximum and minimum load values ​​are identified, and the difference between the two is calculated as the extreme value difference. For example, if a building's maximum load value during the monitoring period is 150 kW and its minimum load value is 20 kW, the calculation process is as follows: The extreme value difference was found to be 130 kilowatts. Simultaneously, a fluctuation judgment threshold was set. This threshold was determined based on: statistically analyzing the load data of known residential buildings within the target area during the nighttime period from 0:00 to 4:00, calculating the standard deviation of the load data for this period, and taking three times this standard deviation as the fluctuation judgment threshold. Assuming the statistically obtained standard deviation of the nighttime load of residential buildings is 5 kilowatts, the fluctuation judgment threshold is calculated as follows: That is, 15 kilowatts. The instantaneous power load sequence is scanned sequentially by time, and the number of times the absolute value of the load change amplitude at adjacent time points exceeds the fluctuation threshold is counted. This number is defined as the fluctuation frequency. If the load in the sequence is 50 kilowatts at a certain moment and 70 kilowatts at the next moment, the change amplitude is... kilowatt, because If the extreme value difference is not recorded, it is counted as one fluctuation. Finally, the calculated extreme value difference is combined with the fluctuation frequency value to construct a load characteristic vector that characterizes the building's electricity consumption behavior.

[0022] S102: Map the load feature vector to a multi-dimensional feature measurement space, call the preset commercial operation load benchmark and residential load benchmark as spatial reference points, calculate the Euclidean distance between the load feature vector and the commercial operation load benchmark and the residential load benchmark in the measurement space, and generate a load attribute difference measurement value. A two-dimensional feature metric space is constructed, defining the extreme value difference as the first dimension and the fluctuation frequency as the second dimension. Historical load data of confirmed standard commercial complexes and standard residential communities are selected from the historical database, and the arithmetic mean of their extreme value difference and fluctuation frequency are calculated respectively. The mean coordinates of the standard commercial complex are set as the commercial operation load benchmark point, and the mean coordinates of the standard residential community are set as the residential load benchmark point. Assume the coordinates of the commercial operation load benchmark point are (200, 100) and the coordinates of the residential load benchmark point are (50, 20). The load feature vector of the building to be detected generated in step S101 is projected into this metric space, assuming the feature vector coordinates of the building to be detected are (180, 90). Using the Euclidean distance calculation logic, the straight-line distance between the feature point coordinates of the building to be detected and the coordinates of the commercial operation load benchmark point is calculated respectively. The specific calculation process is as follows: First, the square of the difference in the first dimension is calculated. Then calculate the square of the difference in the second dimension. Add the two together Finally, perform the square root operation on the sum. The load attribute difference metric is approximately 22.36. Similarly, the distance between this feature point and the residential load benchmark point is calculated. If the residential benchmark point is (50, 20), then the calculation... , The distance is These two distance values ​​will be used as the difference measure for business model and the difference measure for residential pattern, respectively.

[0023] S103: Set filtering conditions based on load attribute difference measurement value, traverse and search for building units matching commercial electricity consumption pattern characteristics from buildings in the target area, extract the unique identifier and spatial location information of the target building unit in the urban land use digital twin model, and mark the attribute status of the target building unit as commercial use to generate commercial operation nodes. A threshold for the difference metric used for screening is set. This threshold is determined by calculating the statistical distribution of the distances from all known commercial building feature points in the metric space to the commercial operating load benchmark point, specifically taking the 95th percentile value of this distance set. Assuming the statistically calculated 95th percentile value of the known commercial building distance set is 30, the difference metric threshold is set to 30. The buildings in the target area are traversed and searched to check their corresponding business model difference metric values. If a building's business model difference metric value is less than the above-mentioned difference metric threshold, and its business model difference metric value is less than its residential model difference metric value, then the building is determined to meet the characteristics of a commercial electricity consumption pattern. Taking the aforementioned calculation as an example, the building's business model difference metric value is 22.36, and its residential model difference metric value is 147.65, and... as well as The building meets the screening criteria. From the data table of the urban land use digital twin model, extract the unique identification code of the eligible building and the longitude and latitude coordinates of its geometric center. In the model's attribute database, update the current land use field of the building from its original state to "commercial use" and mark it as a commercial operation node.

[0024] Please see Figure 3 Step S2 is as follows: S201: Collect mobile terminal location signaling data covering the area where the commercial operation node is located, extract the timestamp sequence and latitude and longitude coordinate sequence from the signaling data stream, calculate the latitude and longitude coordinate displacement distance corresponding to adjacent timestamps to obtain the movement rate, filter latitude and longitude coordinates whose movement rate exceeds the preset speed threshold, and generate a terminal dwell point sequence. By interfacing with the server interface of mobile communication operators, location signaling data of mobile terminals within the administrative region covering commercial operation nodes is collected in batches. The signaling data stream is parsed to separate the anonymous user identifier, timestamp sequence, and corresponding base station triangulation latitude and longitude coordinate sequence. For consecutive records with the same user identifier, the time difference between two adjacent timestamps is calculated, and the displacement distance between the corresponding two latitude and longitude coordinates is obtained using spherical distance calculation logic. Dividing the displacement distance by the time difference yields the terminal's movement speed between the two points. A movement speed threshold is set based on the statistical distribution of average pedestrian movement speeds on urban non-motorized vehicle lanes, using the 99th percentile value of the distribution curve as the threshold. Assuming the 99th percentile of the pedestrian speed distribution is 5 kilometers per hour, the movement speed threshold is set to 5 kilometers per hour. If a terminal moves 1 kilometer in 0.1 hours, the movement speed is calculated as follows: km / h, because If a data point is determined to be in a vehicle-mounted or rapidly passing state, it is removed. If a terminal moves 0.4 kilometers in 0.1 hours, its movement rate is calculated as follows: km / h, because Retain the coordinate point. Reorganize the retained coordinate points into a sequence of terminal dwell points in chronological order.

[0025] S202: Establish a spatial grid based on the commercial operation nodes, map the terminal dwell point sequence to the corresponding spatial grid, traverse and count the cumulative distribution number of dwell points in each spatial grid, calculate the cumulative dwell time of dwell points based on the timestamp difference, associate the cumulative distribution number and cumulative dwell time with the corresponding spatial grid ID, and generate a grid dwell distribution statistics set. Using the geometric center of the commercial operation node as the origin, the surrounding area is divided into grids according to a preset side length scale, generating a spatial grid matrix containing row and column numbers. The preset side length scale is obtained by calculating the square root of the average floor area of ​​buildings within the target area and taking an integer multiple of that value (e.g., 2 times). Assuming the average floor area is 625 square meters, its square root is 25 meters, then the grid side length is set to... Meters. Read each coordinate point in the terminal dwell point sequence, determine its specific row and column position in the grid matrix, and map the coordinate point to the corresponding spatial grid. Traverse each spatial grid and count the total number of dwell points falling within that grid, as the cumulative distribution count. Simultaneously, for consecutive dwell points of the same terminal within the same grid, calculate the difference between the timestamp of its last appearance and the timestamp of its first appearance, and use this difference as the dwell time of that terminal in that grid. Sum the dwell times of all terminals within the grid to obtain the cumulative dwell time of that grid. Assuming there are 3 terminals dwelling in a grid with durations of 10 minutes, 20 minutes, and 30 minutes respectively, the cumulative dwell time is calculated as follows: Minutes. Finally, the cumulative distribution quantity and cumulative dwell time are written into the attribute table of the spatial grid, and a grid dwell distribution statistics set is generated by indexing the grid's unique ID.

[0026] S203: Call the grid dwell distribution statistics set, perform weighted normalization calculation on the cumulative distribution quantity and cumulative dwell time, quantify the degree of population gathering and dwell stickiness in each spatial grid, and assign a corresponding dwell stickiness index to each spatial grid; The residence viscosity index is calculated using a weighted normalization algorithm by calling the grid residence distribution statistics set. This process aims to eliminate the influence of data with different dimensions and introduces a smoothing parameter to prevent denominator invalidation caused by extreme values. The residence viscosity index is defined. The calculation formula is: ,in, Represents the calculated first... The dwelling viscosity index of a spatial grid; Indicates the first The cumulative distribution count obtained from statistics within each grid; It is the minimum value of the cumulative distribution in the entire statistical set. These two parameters represent the maximum cumulative distribution count across the entire statistical set and are derived from the statistical analysis of the entire grid data. Indicates the first The cumulative dwell time obtained from statistics within each grid; This is the minimum cumulative stay duration. This represents the maximum cumulative duration of stay. The numerical stability smoothing parameter is introduced to prevent the denominator from being zero, and is usually taken as a very small positive number (such as 0.001). and These are the cumulative distribution quantity weight and the cumulative stay duration weight, respectively, and satisfy the following conditions: . and The basis for setting the index is principal component analysis: A variance analysis is performed on the pedestrian flow and dwell time in historical commercial hotspots, and the normalized contribution rates of the two to the variance of the commercial vitality index are used as weights.

[0027] Assume the entire region is statistically concentrated, and the cumulative distribution has the minimum value. maximum value Minimum cumulative stay maximum value (Unit: minutes). Quantitative weights were determined using principal component analysis. Duration weight Set smoothing parameters For the currently computed grid The cumulative distribution number it monitored Cumulative duration of stay Substitute the numerical values ​​into the first part of the formula (the quantity term): Substitute the values ​​into the second part of the formula (duration term): Finally, the residence viscosity index of the mesh is calculated: This result indicates that the grid has a moderately high commercial residency potential.

[0028] Please see Figure 4 Step S3 is as follows: S301: Identify non-commercial spatial grids located adjacent to the geometric boundary of commercial operation nodes from the spatial grid, obtain the residence stickiness index corresponding to the non-commercial spatial grids, filter spatial grids whose residence stickiness index exceeds the preset commercial spillover activation threshold and calculate the geometric center coordinates of the spatial grids, and generate the center coordinates of the spillover activation grids. In the spatial topology of the digital twin model, non-commercial spatial grids adjacent to the geometric boundaries of commercial operation nodes are retrieved. The residency stickiness index of these non-commercial spatial grids is then read. A commercial spillover activation threshold is set, based on the arithmetic mean of the residency stickiness indices of all grids within the identified commercial operation nodes, with 60% of this mean used as the activation threshold. Assuming the average stickiness index of the grids within the commercial operation nodes is 0.8, the activation threshold is calculated as follows: Non-commercial grids with a residency stickiness index greater than the commercial spillover activation threshold are selected. If the stickiness index of a non-commercial grid is 0.62, then... This area is identified as an active region influenced by commercial radiation. For each selected grid, the coordinates of its four vertices are extracted, and the average of the x-coordinates and y-coordinates is calculated to obtain the geometric center coordinates of the grid. This set of coordinates is defined as the center coordinates of the spillover active grid.

[0029] S302: Spatial pairing of the geometric center coordinates of the commercial operation node with the center coordinates of the spillover activation grid, establishing an azimuth connection line from the commercial center to the center of the spillover area, calculating the azimuth angle and distance data between the two points, and determining the expansion trend guiding coordinate pair; Using the geometric center coordinates of the commercial operation node as the starting point and the center coordinates of each spillover activation grid as the ending point, a line is established between the two points in the spatial coordinate system. The angle between this line and true north is calculated using inverse trigonometric functions to obtain the azimuth data; the length of the line is calculated using the distance formula between two points to obtain the distance data. Assuming the commercial center coordinates are (100, 100) and the spillover grid center coordinates are (130, 140), the horizontal distance is... The vertical distance is Distance data is calculated as follows: The azimuth angle is obtained by calculating the arctangent value. Combining the azimuth angle with distance data determines the coordinate pairs guiding the expansion trend.

[0030] S303: Based on the spatial direction and distance defined by the expansion trend guide coordinate pair, analyze the direction of commercial land penetrating into the peripheral non-commercial area, as the expansion guide vector; Analyze the azimuth data in the coordinate pair guided by the expansion trend, and identify the continuous grid sequence traversed by the ray along the direction indicated by the azimuth in the digital twin model. Extract the residence viscosity index of each grid in this sequence, and calculate the gradient of the viscosity index change between adjacent grids. Set a gradient smoothing threshold, which is based on: statistically analyzing the average decay of the residence viscosity index of mature commercial districts along the main road direction at every grid interval (e.g., 50 meters), and taking 1.5 times this decay as the gradient smoothing threshold. Assuming the statistically obtained average decay is 0.1, the threshold is set to... If the resident viscosity exponent of the current grid A in the sequence is 0.8, and the exponent of the next grid B is 0.7, calculate the gradient difference as follows: .because This change is determined to be a normal natural decline in commercial activity, confirming that grid B is on the path of effective commercial function penetration. If the index of grid C in the sequence is 0.6 and the index of grid D is 0.3, the gradient difference is, because The data shows a sharp drop in popularity, indicating it's not a valid penetration path. All grid lines meeting the smoothing threshold are identified as specific commercial function penetration paths, eliminating isolated high-population points to ensure the continuity and rationality of the expansion direction.

[0031] Please see Figure 5 Step S4 is as follows: S401: Based on the expansion guide vector and the residence viscosity index, calculate the outward extension distance of the boundary of the commercial operation node, generate outward radiating virtual rays in combination with the initial position of the boundary, define the initial trajectory and magnitude of the boundary expansion, and generate simulated displacement parameters. The process of calculating the outward extrapolation distance from the boundary of a business operation node, based on the expansion guidance vector and the residency stickiness index, is as follows: Extract the geometric modulus of the expansion guide vector as the basic expansion step size, and at the same time obtain the preset global maximum viscosity reference value; Dividing the residence viscosity index by the global maximum viscosity reference value yields the dimensionless weighting coefficient characterizing the expansion intensity. The outward extrapolation distance is obtained by using a weighted product of the basic expansion step size with dimensionless weight coefficients. The emission orientation of the virtual ray is determined based on the direction cosine of the expansion guiding vector, and the termination point of the virtual ray is set by combining the outward extrapolation distance, thus obtaining simulated displacement parameters including the starting point coordinates, the ending point coordinates, and the emission orientation. The geometric modulus of the expansion guide vector is extracted as the basic expansion step size. Simultaneously, the resident viscosity index of all grids in the entire city digital twin model is traversed, and the global maximum value is selected as the global maximum viscosity reference value. Assuming the modulus of the expansion guide vector is 50 meters and the global maximum viscosity reference value is 1.0, the resident viscosity index of the target grid in the current expansion direction is read, assumed to be 0.8. The resident viscosity index is divided by the global maximum viscosity reference value to obtain the dimensionless weighting coefficient characterizing the expansion intensity, calculated as follows: The outward extrapolation distance is calculated by using a dimensionless weighted coefficient to perform a weighted product of the basic expansion step size. The calculation process is as follows: The virtual ray's emission direction is determined based on the direction cosine of the expansion guiding vector. The termination point of the virtual ray is set by combining the outward extrapolation of a 40-meter distance, resulting in simulated displacement parameters including the starting point coordinates, the ending point coordinates, and the emission direction.

[0032] S402: Control the virtual ray corresponding to the simulated displacement parameters to extend in the digital twin model of urban land use, detect whether the virtual ray intersects with the water flow boundary layer or the road median layer in space, if the intersection occurs, extract the coordinates of the intersection point between the virtual ray and the boundary of each type of layer, mark the physical restriction position of the blockage expansion, and generate physical blockage information; The simulated displacement parameters correspond to the extension of virtual rays in the urban land use digital twin model. A water flow boundary layer (containing coordinates of water flow and lake edges) and a road median layer (containing coordinates of highways and railway guardrails) are loaded. A spatial geometric intersection detection algorithm is executed to determine whether the generated virtual ray segment intersects with any linear features in the aforementioned layers. If the virtual ray extends from coordinates (100, 100) to (140, 140), and a water flow boundary line exists along this path, the coordinates of the intersection point between the virtual ray and the water flow boundary line are calculated. If an intersection point is detected, its coordinates are extracted, and the physical constraint location for blocking expansion is marked, generating physical blocking information.

[0033] S403: Based on the intersection coordinates of the physical blocking information, the simulated displacement parameters are truncated to restrict the virtual ray from crossing the physical isolation zone. The displacement endpoint in the blocked direction is adjusted to the intersection position, the original displacement parameters in the unblocked direction are retained, the magnitude and shape of the boundary expansion are corrected, and a corrected displacement vector is generated. The simulated displacement parameters are truncated based on the intersection coordinates of the physical obstruction information. If the intersection coordinates calculated in step S402 are (120, 120), it indicates that the expansion path is blocked at (120, 120). In this case, the endpoint of the displacement in the obstructed direction is adjusted from the original (140, 140) to the intersection position (120, 120). The original displacement parameters in the unobstructed direction are retained. The magnitude and shape of the boundary expansion are corrected to generate a corrected displacement vector. This vector represents the displacement from the starting point (100, 100) to the corrected endpoint (120, 120), ensuring that the expansion behavior does not cross the physical isolation zone.

[0034] Please see Figure 6 Step S5 is as follows: S501: Extract the vector boundary coordinates of commercial operation nodes in the digital twin model of urban land use, apply the data of the corrected displacement vector to the corresponding boundary nodes, calculate the new spatial coordinate positions of the boundary nodes after correction, connect the new coordinate points in sequence to close and form the updated polygon boundary, and generate the reconstructed boundary node sequence. The process of applying the corrected displacement vector data to the corresponding boundary nodes and calculating the new spatial coordinates of the boundary nodes after correction is as follows: Iterate through and extract the x and y coordinates of the vector boundary coordinates, and analyze and correct the horizontal and vertical components of the displacement vector in the plane coordinate system. Perform a linear superposition operation of the horizontal coordinate value and the horizontal component, and a linear superposition operation of the vertical coordinate value and the vertical component to obtain the target displacement coordinates corresponding to the boundary node. The target displacement coordinates are reorganized according to the original connection order of the boundary nodes, and polygon topology rule detection is performed to remove abnormal coordinate points that cause the edge lines to self-intersect, and to establish the new spatial coordinate position. Extract the vector boundary coordinates of commercial operation nodes in the urban land use digital twin model. Apply the corrected displacement vector data to the corresponding boundary nodes and calculate the new spatial coordinates of the boundary nodes after correction. Assume that the original coordinates of a boundary node are (100, 100), and the corresponding corrected displacement vector has a horizontal component of 20 and a vertical component of 20. Iterate through and extract the x-coordinate and y-coordinate values ​​from the vector boundary coordinates, and perform a linear superposition operation of the x-coordinate value and the horizontal component to calculate the result. Perform a linear superposition operation on the vertical axis value and the vertical component, and calculate as follows: The target displacement coordinates corresponding to the boundary nodes are obtained as (120, 120). The target displacement coordinates are reorganized according to the original connection order of the boundary nodes, and polygon topology rule detection is performed to remove abnormal coordinate points that cause the edge lines to self-intersect. The new spatial coordinate positions are established, and the new coordinate points are connected in sequence to close the updated polygon boundary, generating a reconstructed boundary node sequence.

[0035] S502: Obtain the geographic spatial range enclosed by the sequence of reconstructed boundary nodes, incorporate the non-commercial spatial grids covered by the geographic spatial range into commercial land, change the corresponding land use attribute fields, establish the geometric shape and area of ​​the newly added commercial land, and generate the commercial spillover zone. The computational geometry algorithm is used to obtain the geospatial area enclosed by the reconstructed boundary node sequence. The non-commercial spatial grids covered by this geospatial area are then incorporated into the commercial land use area. The corresponding land use attribute fields are modified, and the geometric shape and area of ​​the newly added commercial land are established. Assuming the original commercial land area is 10,000 square meters and the total area of ​​the newly covered grids is 2,000 square meters, the updated total area of ​​the commercial spillover zone is... square meters.

[0036] S503: Overlay the commercial spillover area onto the map layer of the digital twin visual interface, update the coverage of the rendered commercial land, and generate simulation results of the evolution of commercial land morphology. Commercial spillover areas are overlaid onto the map layer of the digital twin visual interface to update the coverage of commercial land and generate simulation results of commercial land morphology evolution. Using a graphics rendering engine, color gradation mapping technology is employed to visualize the newly added commercial areas, based on the residency viscosity index of the grid within that area. An RGB color channel mapping formula is defined to linearly map the residency viscosity index to a red channel value. Assume a newly added grid has a residency viscosity index of 0.8 and a maximum color depth of 255. The calculation process is as follows: If the fill color of this grid is set to RGB(204, 0, 0), it represents a high-heat commercial expansion area. If the index of another grid is 0.3, calculate... (Rounded down to 76), the color is set to RGB(76, 0, 0), representing the low-heat expansion area. The reconstructed boundary is divided into multiple rendering primitives using a polygon triangulation algorithm, and the primitive pixels are filled according to the calculated color values, intuitively presenting the expansion pattern and internal vitality differences of commercial land under the drive of pedestrian traffic and physical environment constraints.

[0037] Please see Figure 7 A digital twin-based urban land use simulation system, comprising: The load characteristic analysis module obtains the instantaneous power load sequence of the building, constructs the load characteristic vector, calculates the difference with the commercial and residential load benchmarks, and marks the commercial operation nodes; The dwell stickiness calculation module collects terminal signaling data from commercial operation nodes, calculates the movement rate and generates dwell points based on the timestamps and coordinates of the terminal signaling data, and calculates the dwell stickiness index by counting the number and duration of dwell points. The expansion guide vector generation module identifies the non-commercial grid surrounding the commercial operation node, extracts the geometric center of the non-commercial grid based on the residency stickiness index, and constructs the expansion guide vector based on the geometric center and the commercial operation node. The boundary displacement correction module calculates the simulated displacement parameters of the commercial operation node boundary based on the expansion guide vector and the residence viscosity index, and corrects the simulated displacement parameters and generates the corrected displacement vector based on the physical blocking information detected by the water flow boundary and the road median layer. The commercial spillover zone reconstruction module reconstructs the boundaries of commercial operation nodes based on the corrected displacement vector, generates a commercial spillover zone covering the non-commercial grid, and renders the simulation results of commercial land morphology evolution.

[0038] The above are merely preferred embodiments of the present invention and are not intended to limit the present invention in any other way. Any person skilled in the art may make changes or modifications to the above-disclosed technical content to create equivalent embodiments that can be applied to other fields. However, any simple modifications, equivalent changes, and modifications made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the protection scope of the present invention.

Claims

1. A method for simulating urban land use based on digital twins, characterized in that, Includes the following steps: S1: Obtain the instantaneous power load sequence of the building, construct the load feature vector and calculate the difference with the commercial and residential load benchmarks, and mark the commercial operation nodes; S2: Collect terminal signaling data of the commercial operation node, calculate the movement rate and generate dwell points based on the timestamp and coordinates of the terminal signaling data, and calculate the dwell stickiness index by counting the number and duration of dwell points. S3: Identify the non-commercial grid surrounding the commercial operation node, extract the geometric center of the non-commercial grid with the residency stickiness index, and construct an expansion guide vector based on the geometric center and the commercial operation node; S4: Calculate the simulated displacement parameters of the boundary of the commercial operation node based on the expansion guide vector and the residence viscosity index; detect physical blocking information based on the water flow boundary and road median layer; correct the simulated displacement parameters; and generate a corrected displacement vector. S5: Reconstruct the boundary of the commercial operation node based on the modified displacement vector, generate a commercial spillover area covering the non-commercial grid, and render the simulation results of the commercial land morphology evolution.

2. The urban land use simulation method based on digital twins according to claim 1, characterized in that, The commercial operation node includes a unique index identifier for the target building and a commercial function attribute label determined based on the load characteristic difference. The dwell stickiness index is specifically the regional population attraction calculated based on the total number of dwell points and the total dwell time within the spatial grid. The expansion guidance vector includes the vector starting coordinates anchored to the center of the commercial operation node and the guiding orientation attribute pointing to the geometric center of the non-commercial spatial grid. The corrected displacement vector includes the displacement direction that avoids the obstruction of water flow boundaries and road medians and the displacement modulus corrected based on physical obstruction information. The commercial land morphology evolution simulation results include the reconstructed vector polygon vertex coordinate sequence, the updated commercial spillover zone geometric layer, and the rendered dynamic evolution visual view of the land boundary.

3. The urban land use simulation method based on digital twins according to claim 1, characterized in that, Step S1 is as follows: S101: Obtain urban geographic vector data and building energy operation logs of the target area, construct a digital twin model of urban land use based on the urban geographic vector data and map energy data, extract the instantaneous power load sequence of buildings in the target area within a preset monitoring period in the model, perform extreme value difference calculation and fluctuation frequency statistics on the instantaneous power load sequence, and generate a load feature vector. S102: Map the load feature vector to a multi-dimensional feature measurement space, call the preset commercial operation load benchmark and residential load benchmark as spatial reference points, calculate the Euclidean distance between the load feature vector and the commercial operation load benchmark and the residential load benchmark in the measurement space, and generate a load attribute difference measurement value. S103: Based on the load attribute difference measurement value, set the filtering conditions, traverse and retrieve building units that match the characteristics of commercial electricity consumption patterns from the buildings in the target area, extract the unique identifier and spatial location information of the target building unit in the urban land use digital twin model, and mark the attribute status of the target building unit as commercial use to generate a commercial operation node.

4. The urban land use simulation method based on digital twins according to claim 1, characterized in that, Step S2 is as follows: S201: Collect mobile terminal location signaling data covering the area where the commercial operation node is located, extract the timestamp sequence and latitude and longitude coordinate sequence from the signaling data stream, calculate the latitude and longitude coordinate displacement distance corresponding to adjacent timestamps to obtain the movement rate, filter latitude and longitude coordinates whose movement rate exceeds a preset speed threshold, and generate a terminal dwell point sequence. S202: Establish a spatial grid based on the commercial operation nodes, map the terminal dwell point sequence to the corresponding spatial grid, traverse and count the cumulative distribution number of dwell points in each spatial grid, calculate the cumulative dwell time of dwell points based on the timestamp difference, associate the cumulative distribution number and cumulative dwell time with the corresponding spatial grid ID, and generate a grid dwell distribution statistics set. S203: Call the grid dwell distribution statistics set, perform weighted normalization calculation on the cumulative distribution quantity and cumulative dwell time, quantify the degree of crowd gathering and dwell stickiness in each spatial grid, and assign a corresponding dwell stickiness index to each spatial grid.

5. The urban land use simulation method based on digital twins according to claim 1, characterized in that, Step S3 is as follows: S301: Identify non-commercial spatial grids located adjacent to the geometric boundary of the commercial operation node from the spatial grid, obtain the residence stickiness index corresponding to the non-commercial spatial grid, filter spatial grids whose residence stickiness index exceeds the preset commercial spillover activation threshold and calculate the geometric center coordinates of the spatial grids, and generate the center coordinates of the spillover activation grids. S302: Spatially pair the geometric center coordinates of the commercial operation node with the center coordinates of the spillover activation grid, establish an azimuth connection line from the commercial center to the center of the spillover area, calculate the azimuth angle and distance data between the two points, and determine the expansion trend guiding coordinate pair; S303: Based on the spatial direction and distance defined by the expansion trend guide coordinate pair, analyze the direction of commercial land penetrating into the surrounding non-commercial areas, and use it as the expansion guide vector.

6. The urban land use simulation method based on digital twins according to claim 1, characterized in that, Step S4 is as follows: S401: Based on the expansion guide vector and the residence viscosity index, calculate the outward extension distance of the boundary of the commercial operation node, generate outward radiating virtual rays in combination with the initial position of the boundary, define the initial trajectory and amplitude of the boundary expansion, and generate simulated displacement parameters. S402: Control the virtual ray corresponding to the simulated displacement parameter to extend in the digital twin model of urban land use, detect whether the virtual ray intersects with the water flow boundary layer or the road isolation strip layer in space, if the intersection occurs, extract the coordinates of the intersection point of the virtual ray with the boundary of each type of layer, mark the physical restriction position of the blockage expansion, and generate physical blockage information; S403: Based on the intersection coordinates of the physical blocking information, the simulated displacement parameters are truncated to restrict the virtual ray from crossing the physical isolation zone. The displacement endpoint in the blocked direction is adjusted to the intersection position, the original displacement parameters in the unblocked direction are retained, the magnitude and shape of the boundary expansion are corrected, and a corrected displacement vector is generated.

7. The urban land use simulation method based on digital twins according to claim 6, characterized in that, The process of calculating the outward extrapolation distance from the boundary of the business operation node based on the expansion guidance vector and the residence stickiness index is as follows: Extract the geometric modulus of the expansion guide vector as the basic expansion step size, and at the same time obtain the preset global maximum viscosity reference value; Dividing the residence viscosity index by the global maximum viscosity reference value yields the dimensionless weighting coefficient characterizing the expansion intensity. The outward extrapolation distance is obtained by using a weighted product of the basic expansion step size with dimensionless weight coefficients. The emission orientation of the virtual ray is determined based on the direction cosine of the expansion guiding vector, and the termination point of the virtual ray is set by combining the outward extrapolation distance, thus obtaining simulated displacement parameters including the starting point coordinates, the ending point coordinates, and the emission orientation.

8. The urban land use simulation method based on digital twins according to claim 1, characterized in that, Step S5 is as follows: S501: Extract the vector boundary coordinates of the commercial operation node in the digital twin model of urban land use, apply the data of the corrected displacement vector to the corresponding boundary node, calculate the new spatial coordinate position of the boundary node after correction, connect the new coordinate points in sequence to close and form the updated polygon boundary, and generate a reconstructed boundary node sequence. S502: Obtain the geographic spatial range enclosed by the reconstructed boundary node sequence, incorporate the non-commercial spatial grid covered by the geographic spatial range into commercial land, change the corresponding land use attribute fields, and establish the geometric shape and area of ​​the newly added commercial land to generate a commercial spillover zone. S503: Overlay the commercial spillover area onto the map layer of the digital twin visual interface, update the coverage of the rendered commercial land, and generate simulation results of the evolution of commercial land morphology.

9. The urban land use simulation method based on digital twins according to claim 8, characterized in that, The process of applying the corrected displacement vector data to the corresponding boundary nodes and calculating the new spatial coordinate positions of the boundary nodes after correction is as follows: Iterate through and extract the x and y coordinates of the vector boundary coordinates, and analyze and correct the horizontal and vertical components of the displacement vector in the plane coordinate system. Perform a linear superposition operation of the horizontal coordinate value and the horizontal component, and a linear superposition operation of the vertical coordinate value and the vertical component to obtain the target displacement coordinates corresponding to the boundary node. The target displacement coordinates are reorganized according to the original connection order of the boundary nodes, and polygon topology rule detection is performed to remove abnormal coordinate points that cause the edge lines to self-intersect, thus establishing the new spatial coordinate position.

10. A digital twin-based urban land use simulation system, characterized in that, The system, which implements the urban land use simulation method based on digital twins according to any one of claims 1-9, comprises: The load characteristic analysis module obtains the instantaneous power load sequence of the building, constructs the load characteristic vector, calculates the difference with the commercial and residential load benchmarks, and marks the commercial operation nodes; The dwell stickiness calculation module collects terminal signaling data from the commercial operation nodes, calculates the movement rate and generates dwell points based on the timestamps and coordinates of the terminal signaling data, and calculates the dwell stickiness index by counting the number and duration of dwell points. The expansion guide vector generation module identifies the non-commercial grid surrounding the commercial operation node, extracts the geometric center of the non-commercial grid with the residency stickiness index, and constructs the expansion guide vector based on the geometric center and the commercial operation node; The boundary displacement correction module calculates the simulated displacement parameters of the boundary of the commercial operation node based on the expansion guide vector and the residence viscosity index, and corrects the simulated displacement parameters based on the physical blocking information detected by the water flow boundary and the road median layer, thereby generating a corrected displacement vector. The commercial spillover zone reconstruction module reconstructs the boundaries of commercial operation nodes based on the corrected displacement vector, generates a commercial spillover zone covering the non-commercial grid, and renders the simulation results of commercial land morphology evolution.

Images