Mountain park ecological restoration digital twin method and system

Abstract

The present application relates to natural system simulation technical field, specifically to a mountain park ecological restoration digital twin method and system, in the present application, by constructing multi-directional elevation profile sequence and calculating micro-landform fractal dimension, combining with ellipse fitting technology to generate roughness tensor matrix, accurately quantifying the structural difference of micro-landform in different directions, using water flow velocity vector projection in roughness tensor to analyze equivalent fractal dimension, realizing accurate calculation of surface runoff resistance coefficient dynamic response with flow direction, thereby significantly improving the solving accuracy of hydrodynamic equation under complex terrain, at the same time, introducing environmental biological stress intensity index, based on spatial Euclidean distance as distance attenuation factor and interspecific competition coefficient weighted summation, deeply coupling physical habitat hydrodynamic condition and biological community interspecific competition mechanism, using biological stress correction state transition probability, realizing dynamic prediction of vegetation community succession path.

Description

Technical Field

[0001] This invention relates to the field of natural system simulation technology, and in particular to a digital twin method and system for ecological restoration of mountain parks. Background Technology

[0002] The field of natural system simulation technology refers to a set of technologies that use mathematical models, physical mechanism models and computer simulation methods to digitally represent, extrapolate the temporal evolution and reproduce the state of natural elements such as mountains, hydrology, soil, vegetation and biological communities and their interaction processes. This field of technology usually covers natural environment data acquisition, digital modeling of mountain topography and geomorphology, description of ecological process mechanisms, simulation of environmental evolution processes, and virtual-real mapping and dynamic updating.

[0003] Among them, the digital twin method for ecological restoration of mountain parks refers to the process of acquiring elevation data, slope and aspect data, and land cover type through field surveys around the ecological restoration object of the mountain park. A three-dimensional mountain terrain model is constructed based on discrete elevation points. Combined with historical meteorological records, soil physicochemical parameters, and vegetation distribution survey results, the vegetation configuration, slope reinforcement measures, and drainage paths are entered into the simulation system as parameters according to the preset ecological restoration plan. The model state is corrected by periodically updating the monitoring data, so as to realize the corresponding expression and synchronous description of the ecological restoration process of the mountain park in the digital space.

[0004] Existing technologies primarily rely on macroscopic topographic factors and static vegetation survey data to construct models, and only use periodic monitoring data to correct the model state, ignoring the anisotropic influence of microscopic geomorphological features on surface processes. This results in the inability to distinguish the differences in surface roughness under different flow directions when simulating hydrodynamic processes, making it difficult to set the surface runoff resistance coefficient without physical mechanism support. Consequently, the accuracy of water flow velocity and water depth calculations is insufficient, and the spatial competition pressure between neighboring vegetation and its dynamic constraints on community succession are not quantified. This means that the simulation of ecological restoration processes lacks the coupling mechanism between hydrodynamics and biological stress, making it difficult to truly reflect the natural succession laws of vegetation and long-term restoration effects under complex habitats, leading to deviations between simulation results and actual ecological evolution paths. Summary of the Invention

[0005] The purpose of this invention is to address the shortcomings of existing technologies by proposing a digital twin method and system for ecological restoration of mountain parks.

[0006] To achieve the above objectives, the present invention adopts the following technical solution: a digital twin method for ecological restoration of mountain parks, comprising the following steps: S1: Construct a digital elevation model of the mountain terrain using lidar mapping equipment, divide the mountain terrain grid area in the model, extract surface elevation data from multiple directions, and generate a multi-directional elevation profile sequence. S2: Calculate the micro-topographic fractal dimension of the multi-azimuth elevation profile sequence, and fit the micro-topographic fractal dimension to construct the roughness tensor matrix; S3: Establish the fluid control equation and obtain the water flow velocity vector. Project the water flow velocity vector onto the roughness tensor matrix to calculate the equivalent fractal dimension, determine the surface runoff resistance coefficient, update the fluid control equation, and obtain the updated hydrodynamic parameters. S4: Identify vegetation species type identifiers within the mountain topographic grid area and its neighboring areas. Based on the spatial Euclidean distance between the neighboring mountain topographic grid area and the central area, use the spatial Euclidean distance as a distance decay factor to perform a weighted summation of the preset interspecific competition coefficients to obtain the environmental biological stress intensity index. S5: Based on the updated hydrodynamic parameters, match the theoretical state transition probability, correct the theoretical state transition probability according to the environmental biological stress intensity index, obtain the actual growth transition probability, determine the vegetation community succession state at the next moment, and generate digital twin simulation results for the ecological restoration of the mountain park.

[0007] As a further aspect of the present invention, the multi-azimuth elevation profile sequence includes discrete azimuth markers of grid computing units and a set of surface elevation values ​​distributed along the azimuth direction; the roughness tensor matrix includes the principal axis deflection angle of the micro-topographic fractal fitting ellipse, the fractal dimension components in the first principal axis direction and the fractal dimension components in the second principal axis direction; the updated hydrodynamic parameters include the velocity vector corrected by the surface runoff resistance coefficient and the updated water depth value; the environmental biological stress intensity index is specifically a scalar value obtained by weighted summation of the interspecific competition inhibition coefficient based on the distance attenuation factor between the neighboring mountain topographic grid area and the central area; and the digital twin simulation results of the mountain park ecological restoration include the vegetation community succession state of the grid area at the next moment and the spatial distribution mapping information of the state on the mountain topographic digital elevation model.

[0008] As a further aspect of the present invention, the steps for obtaining the multi-azimuth elevation profile sequence are specifically as follows: S111: Collect point cloud data of the target area of ​​the mountain park, convert the three-dimensional spatial coordinate information in the point cloud data into a digital elevation model of the mountain terrain, perform regional grid segmentation processing on the digital elevation model of the mountain terrain, and obtain the mountain terrain grid area. S112: Locate the geometric center of the mountain terrain grid area, establish a plane coordinate system with the geometric center as the origin, set multiple ray directions on the plane with a preset angle interval, and obtain discrete azimuth directions; S113: Using the geometric center of the mountain terrain grid area as the extraction origin, extract the surface elevation data corresponding to the surface of the digital elevation model of the mountain terrain one by one along each of the discrete azimuth directions and arrange them according to spatial distance to generate a multi-azimuth elevation profile sequence.

[0009] As a further aspect of the present invention, the step of obtaining the roughness tensor matrix specifically includes: S211: Call the box-counting dimension algorithm, set grid boxes with multiple side lengths as the measurement scale, cover the surface elevation data curves in the multi-azimuth elevation profile sequence, count the number of non-empty boxes required to cover the curve at each box side length scale, establish a linear regression relationship between the logarithm of the reciprocal of the box side length and the logarithm of the number of boxes, perform an overall correlation analysis between the natural logarithm of the reciprocal of the box side length and the natural logarithm of the corresponding number of non-empty boxes, and combine the cooperative change characteristics of the two in all scale levels and the discrete distribution characteristics of the natural logarithm of the reciprocal of the box side length to extract the regression slope and calculate the micro-topographic fractal dimension; S212: For the fractal dimension of the micro-topography corresponding to each discrete azimuth direction, the least squares fitting algorithm is used to approximate the geometric shape of the ellipse, and the length of the major axis, the length of the minor axis, and the rotation angle of the major axis relative to the preset reference coordinate system are calculated to obtain the geometric feature parameters of the ellipse fitting. S213: Call the major axis length, minor axis length, and rotation angle from the geometric feature parameters of the ellipse fitting, define the component weights of the surface roughness in each direction of the plane coordinate system, and generate the roughness tensor matrix.

[0010] As a further aspect of the present invention, the step of obtaining the updated hydrodynamic parameters specifically comprises: S311: Establish fluid control equations based on the physical spatial properties of the mountain topographic grid region, set the current simulation time step, perform numerical discretization to solve the fluid control equations, obtain the velocity component values ​​of fluid particles at the current time and perform vector synthesis to obtain the water flow velocity vector. S312: Project the water flow velocity vector onto the principal axis coordinate system defined by the roughness tensor matrix, calculate the roughness fractal dimension value matching the current water flow direction, quantify the surface friction effect of the mountain topographic grid area under the current specified flow direction based on the roughness fractal dimension value, and calculate the surface runoff resistance coefficient. S313: The fluid control equations are corrected using the surface runoff resistance coefficients. The fluid state variables in the fluid control equations are recalculated based on the current simulation time step to obtain the updated hydrodynamic parameters.

[0011] As a further aspect of the present invention, the step of obtaining the environmental biological stress intensity index specifically includes: S411: Identify the vegetation species type identifiers growing in the mountain terrain grid area within the current simulation time step, and simultaneously scan the neighboring vegetation species type identifiers growing in the neighborhood of the mountain terrain grid area, and set the interspecific competition coefficient between the neighboring identifiers and the central area identifiers. S412: Calculate the Euclidean distance between the geometric center of the neighboring mountain terrain grid region and the geometric center of the current mountain terrain grid region to obtain the spatial Euclidean distance; S413: The spatial Euclidean distance is used as a distance decay factor and the interspecific competition coefficient is weighted and summed to quantify the competitive pressure exerted by the neighboring vegetation on the central vegetation, thus obtaining the environmental biological stress intensity index.

[0012] As a further aspect of the present invention, the steps for obtaining the digital twin simulation results of the ecological restoration of the mountain park are as follows: S511: Call the preset vegetation succession probability matrix, use the updated hydrodynamic parameters and the vegetation species type identifier in the mountain topography grid area as the joint index key value, search for the matching succession rule in the vegetation succession probability matrix, extract the probability value of vegetation evolving from the current state to the next succession stage, and obtain the theoretical state transition probability. S512: Using the environmental biological stress intensity index as a probability penalty coefficient, calculate the difference between the theoretical state transition probability and the environmental biological stress intensity index to correct the succession probability, obtain the actual growth transition probability, compare it with the preset random judgment threshold, and determine the vegetation community succession state of the mountain terrain grid area in the next simulation time step. S513: Map the vegetation community succession state to the grid of the digital elevation model of the mountain terrain, perform spatial fusion of vegetation dynamic evolution information and static geographical information of terrain, and generate digital twin simulation results of ecological restoration of the mountain park.

[0013] A digital twin system for ecological restoration of a mountain park, the system comprising: The terrain construction module uses lidar mapping equipment to build a digital elevation model of the mountain terrain, divides the mountain terrain grid area in the model and extracts surface elevation data from multiple directions to generate a multi-directional elevation profile sequence. The roughness modeling module calculates the micro-topographic fractal dimension of the multi-azimuth elevation profile sequence, and performs fitting processing on the micro-topographic fractal dimension to construct a roughness tensor matrix. The hydrodynamic parameter correction module establishes the fluid control equation and obtains the water flow velocity vector. It projects the water flow velocity vector onto the roughness tensor matrix to calculate the equivalent fractal dimension, determines the surface runoff resistance coefficient, and updates the fluid control equation to obtain the updated hydrodynamic parameters. The biological stress calculation module identifies the vegetation species type identifiers in the mountain topographic grid area and its neighboring areas. Based on the spatial Euclidean distance between the neighboring mountain topographic grid area and the central area, the spatial Euclidean distance is used as a distance decay factor to perform a weighted summation on the preset interspecific competition coefficient to obtain the environmental biological stress intensity index. The community succession simulation module matches the theoretical state transition probability with the updated hydrodynamic parameters, corrects the theoretical state transition probability according to the environmental biological stress intensity index, obtains the actual growth transition probability, determines the vegetation community succession state at the next moment, and generates digital twin simulation results for the ecological restoration of the mountain park.

[0014] Compared with the prior art, the advantages and positive effects of the present invention are as follows: In this invention, by constructing a multi-directional elevation profile sequence and calculating the fractal dimension of micro-topography, and combining it with ellipse fitting technology to generate a roughness tensor matrix, the structural differences of micro-topography in different directions are accurately quantified. By projecting the water flow velocity vector into the roughness tensor to analyze the equivalent fractal dimension, the dynamic response of surface runoff resistance coefficient with flow direction is accurately calculated, thereby significantly improving the solution accuracy of hydrodynamic equations under complex terrain. At the same time, an environmental biological stress intensity index is introduced, and the spatial Euclidean distance is used as a distance attenuation factor and weighted summed with the interspecific competition coefficient to deeply couple the hydrodynamic conditions of physical habitat with the interspecific competition mechanism of biological community. By using the biological stress correction theory state transition probability, the dynamic prediction of vegetation community succession path is realized, effectively improving the realism of the entire process of ecological restoration of mountain parks in digital space and the ability to express evolutionary laws. Attached Figure Description

[0015] Figure 1 This is a schematic diagram of the workflow of the present invention; Figure 2 This is a flowchart of the steps for obtaining the multi-directional elevation profile sequence of the present invention; Figure 3 This is a flowchart of the roughness tensor matrix construction steps of the present invention; Figure 4 This is a flowchart of the updated hydrodynamic parameter acquisition steps of the present invention; Figure 5 This is a flowchart of the steps for calculating the environmental biological stress intensity index of the present invention; Figure 6 This is a flowchart illustrating the steps involved in generating digital twin simulation results for ecological restoration of a mountain park, as described in this invention. Detailed Implementation

[0016] 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.

[0017] Please see Figure 1 This invention provides a technical solution, a digital twin method for ecological restoration of mountain parks, comprising the following steps: S1: Construct a digital elevation model of the mountain terrain using lidar mapping equipment, divide the mountain terrain grid area in the model, extract surface elevation data from multiple directions, and generate a multi-directional elevation profile sequence. S2: Calculate the fractal dimension of micro-topography in the multi-azimuth elevation profile sequence, and fit the fractal dimension of micro-topography to construct the roughness tensor matrix; S3: Establish the fluid control equations and obtain the water flow velocity vector. Project the water flow velocity vector onto the roughness tensor matrix to calculate the equivalent fractal dimension, determine the surface runoff resistance coefficient, update the fluid control equations, and obtain the updated hydrodynamic parameters. S4: Identify vegetation species type identifiers within the mountain topographic grid area and its neighboring areas. Based on the spatial Euclidean distance between the neighboring mountain topographic grid area and the central area, use the spatial Euclidean distance as a distance decay factor to perform a weighted summation of the preset interspecific competition coefficients to obtain the environmental biological stress intensity index. S5: Based on the updated hydrodynamic parameters, the theoretical state transition probability is matched and corrected according to the environmental biological stress intensity index to obtain the actual growth transition probability. The vegetation community succession state at the next moment is determined, and the digital twin simulation results of the ecological restoration of the mountain park are generated.

[0018] The multi-azimuth elevation profile sequence includes discrete azimuth markers of grid computing units and a set of surface elevation values ​​distributed along the azimuth direction. The roughness tensor matrix includes the principal axis deflection angle of the micro-topographic fractal fitting ellipse, the fractal dimension components in the first principal axis direction, and the fractal dimension components in the second principal axis direction. The updated hydrodynamic parameters include the velocity vector corrected by the surface runoff resistance coefficient and the updated water depth values. The environmental biological stress intensity index is specifically a scalar value obtained by weighted summation of the interspecific competition inhibition coefficient based on the distance attenuation factor between the neighboring mountain topographic grid area and the central area. The digital twin simulation results of the mountain park ecological restoration include the vegetation community succession status of the grid area at the next moment and the spatial distribution mapping information of the status on the mountain topographic digital elevation model.

[0019] Please see Figure 2 The specific steps for obtaining the multi-azimuth elevation profile sequence are as follows: S111: Collect point cloud data of the target area of ​​the mountain park, convert the three-dimensional spatial coordinate information in the point cloud data into a digital elevation model of the mountain terrain, perform regional grid segmentation processing on the digital elevation model of the mountain terrain, and obtain the mountain terrain grid area. The digital elevation model of the mountain terrain includes a two-dimensional raster matrix covering the target area of ​​the mountain park, constructed based on point cloud data, and the absolute elevation value of the ground surface at each raster position in the two-dimensional raster matrix. The absolute elevation value of the ground surface is calculated by spatial interpolation based on the three-dimensional spatial coordinate information in the point cloud data. An airborne lidar scan is initiated to perform high-density mapping of the target mountain area, acquiring a raw point cloud dataset containing 3D coordinates and their reflection intensity attributes. Statistical filtering and denoising are performed on the raw point cloud data, setting a threshold of 50 neighboring points. This threshold is determined by statistical analysis of historical mountain topographic point cloud samples, calculating the density distribution differences between noise and valid points in different neighborhood ranges, and selecting the density boundary value with the highest discriminative power. The average distance from each point to its neighbors is calculated, and a distance histogram is generated based on the average distance of all points. The global average distance and standard deviation are calculated based on the Gaussian distribution assumption. Points with an average distance greater than the sum of the global average distance and twice the standard deviation are marked as outliers and removed. For example, when processing a specific coordinate point, the average distance between that point and its neighbors is first calculated. Then, this distance is compared with a pre-calculated removal threshold consisting of the global mean plus twice the standard deviation. If the average distance of the point is greater than this threshold, it is determined to be a noise point and removed from the dataset. Subsequently, the inverse distance weighted interpolation algorithm is invoked, with a search radius set to 5 meters. This radius is determined based on half of the range value in the terrain feature scale variogram analysis. Valid point cloud data within the neighborhood of the raster to be interpolated is searched, and the Euclidean distance between each valid point cloud data point and the raster center is calculated. The reciprocal square of the Euclidean distance is used as a weight to assign to the corresponding absolute elevation value of the ground surface, and a weighted summation operation is performed. For example, when interpolating a blank raster, multiple valid points within the search radius are found. The distances from these points to the raster center are calculated, and weights are generated accordingly. The known elevation of each point is multiplied by its corresponding weight and then summed to obtain the interpolated elevation of that raster location. This model is constructed as a two-dimensional raster matrix with a resolution of 0.5 meters, where each grid precisely stores the absolute elevation value of the ground surface at that location. Perform a regional gridding operation, traverse the rows and columns of the two-dimensional raster matrix, and cut the entire digital elevation model into multiple independent sub-regions according to the preset calculation grid size of 20 meters by 20 meters. Assign a unique index number to each sub-region, and finally obtain a regularized mountain terrain grid region.

[0020] S112: Locate the geometric center of the mountain terrain grid area, establish a plane coordinate system with the geometric center as the origin, set multiple ray directions with preset angle intervals on the plane, and obtain discrete azimuth directions; The boundary coordinate extreme values ​​of the currently processed mountain terrain grid region are read. The center x-coordinate is obtained by calculating half of the sum of the maximum and minimum x-coordinates, and the center y-coordinate is obtained by calculating half of the sum of the maximum and minimum y-coordinates, thus determining the geometric center coordinates. Using this geometric center as the origin, a local plane coordinate system is established based on the Cartesian coordinate system rules, where true north is defined as the zero-degree azimuth. Multiple ray directions with preset angular intervals are set on the plane. The preset angular interval is set based on: statistically analyzing historical terrain complexity data and calculating the standard deviation of terrain slope. There are two possibilities: the first possibility is that when the standard deviation of terrain slope exceeds 20 degrees, the angular interval is set to 10 degrees to capture subtle terrain changes; the second possibility is that when the standard deviation of terrain slope is less than or equal to 20 degrees, the angular interval is set to 15 degrees. In this embodiment, the calculated standard deviation belongs to the second possibility, therefore the preset angular interval is set to 15 degrees. Starting from zero degrees, the interval is accumulated clockwise until it reaches 360 degrees, thus obtaining a discrete azimuth direction sequence. For example, when determining the geometric center of a specific grid, the boundary coordinates of the region in the horizontal and vertical directions are extracted respectively. The center point coordinates are determined by taking the intermediate value. Then, the corresponding angle interval is selected based on the calculated standard deviation of the terrain slope, and a ray sequence containing multiple discrete angles is generated with the center point as the origin.

[0021] S113: Using the geometric center of the mountain terrain grid area as the extraction origin, extract the surface elevation data corresponding to the surface of the digital elevation model of the mountain terrain one by one along each discrete azimuth direction and arrange them according to spatial distance to generate a multi-azimuth elevation profile sequence. A sampling trajectory is generated along each discrete azimuth direction. For a specific discrete azimuth direction, the sampling step size is set to 0.5 meters, based on the minimum resolution of the digital elevation model (DEM) raster, to ensure that the sampling accuracy is not lower than that of the original data. The sampling is extended point by point along the ray direction from the geometric center until the boundary of the mountain terrain grid area is reached. At each sampling point, the surface elevation data at the corresponding coordinates in the DEM is read using bilinear interpolation. Specifically, the four nearest raster points around the sampling point are found, and distance weights are calculated based on the relative positions of the sampling point and these four raster points. The weighted average of the absolute surface elevation values ​​of the four raster points is then used to obtain the precise elevation of the sampling point. All extracted surface elevation data are stored sequentially in an array according to their spatial distance from the geometric center to the boundary, constructing the terrain profile curve data in that direction. The above steps are repeated for all discrete azimuth directions, and finally, a multi-azimuth elevation profile sequence is generated. For example, when extracting a profile along a specific azimuth angle (such as 45 degrees), starting from the center point, a sampling coordinate is determined at fixed step intervals. The weighted average of the elevations of the four grid points around that coordinate is calculated as the elevation at that point, and this series of elevation values ​​are arranged in sequence to form a complete topographic profile data chain in that direction.

[0022] Please see Figure 3 The specific steps for obtaining the roughness tensor matrix are as follows: S211: Call the box-counting dimension algorithm, set grid boxes with multiple side lengths as the measurement scale, cover the surface elevation data curves in the multi-azimuth elevation profile sequence, count the number of non-empty boxes required to cover the curve at each box side length scale, establish a linear regression relationship between the logarithm of the reciprocal of the box side length and the logarithm of the number of boxes, perform an overall correlation analysis between the natural logarithm of the reciprocal of the box side length and the natural logarithm of the corresponding number of non-empty boxes, and combine the cooperative change characteristics of the two in all scale levels and the discrete distribution characteristics of the natural logarithm of the reciprocal of the box side length to extract the regression slope and calculate the micro-topographic fractal dimension. Normalization is performed on each surface elevation data curve in the multi-azimuth elevation profile sequence, mapping spatial distance and surface elevation data to a numerical range of 0 to 1. A grid box with multiple side lengths is set as a metric scale, initially set to half the total data length, then decreasing by powers of 2 to generate a series of grid box sequences at different scales. Each scale of grid box is overlaid on the normalized surface elevation data curve, and a counting process is initiated to count the number of non-empty boxes that intersect with the surface elevation data curve at the current box side length scale. A linear regression relationship is established between the logarithm of the reciprocal of the box side length and the logarithm of the number of boxes. The slope of the regression line is calculated using the least squares method as the micro-topographic fractal dimension, calculated using the following formula: In the formula, This represents the calculated fractal dimension of the micro-topography, and its value is usually greater than 1. The larger the value, the rougher the surface. This represents the total number of scale levels used for regression analysis, determined by taking the logarithm of the ratio of data resolution to total length. Representing the The natural logarithm of the reciprocal of the side length of the box is obtained by dividing 1 by the current side length of the box and taking the natural logarithm. Representative at the The natural logarithm of the number of non-empty boxes required to cover the curve at the level of statistics is obtained by taking the natural logarithm after counting the number of non-empty boxes.

[0023] For example, select Calculations are performed at each scale level, with the first-level box having a side length of [missing value]. The reciprocal of the side length of the box is , Calculated as The count revealed 500 non-empty boxes. Calculated as The second-level box has a side length of 0.05, and the reciprocal of the side length is 20. Calculated as Statistics show that non-empty boxes are obtained. indivual, Calculated as .

[0024] Substituting into the formula, the molecule is calculated as follows: The denominator is calculated as follows: The final calculation yielded... .

[0025] S212: For the fractal dimension of the micro-topography corresponding to each discrete azimuth direction, the least squares fitting algorithm is used to approximate the geometric shape of the ellipse. The length of the major axis, the length of the minor axis, and the rotation angle of the major axis relative to the preset reference coordinate system of the fitted ellipse are calculated to obtain the geometric feature parameters of the fitted ellipse. For each discrete azimuth direction, the fractal dimension of the micro-topography is plotted in a polar coordinate system with the geometric center as the origin. The polar angle corresponds to the discrete azimuth direction, and the polar radius corresponds to the fractal dimension value. A least-squares fitting algorithm is used to approximate the ellipse's geometry, constructing the general equation of the ellipse. The objective function is defined as the sum of the squares of the algebraic distances from all data points to the ellipse boundary. By taking the partial derivatives of the ellipse parameters (center coordinates, major axis length, minor axis length, and rotation angle) and setting them to zero, the system of equations is solved iteratively to minimize the objective function. The major axis length, minor axis length, and rotation angle of the major axis relative to a preset reference coordinate system (e.g., true north) of the best-fit ellipse are calculated. These parameters are obtained as the geometric feature parameters of the ellipse fitting. For example, after completing the fitting operation for all azimuth fractal dimension points, the geometric parameters of the best-fit ellipse are output, including the specific lengths of the major and minor axes, and the specific angle of clockwise rotation of the major axis relative to true north. These parameters collectively describe the anisotropic characteristics of the terrain roughness.

[0026] S213: Call the major axis length, minor axis length, and rotation angle from the geometric feature parameters of the ellipse fitting, define the component weights of the surface roughness in each direction of the plane coordinate system, and generate the roughness tensor matrix. We define the component weights of surface roughness in each direction of a planar coordinate system and construct a second-order symmetric tensor based on the geometric properties of an ellipse. First, we calculate the rotation matrix, constructing a two-dimensional rotation matrix using the cosine and sine values ​​of the rotation angle. Next, we construct a diagonal matrix, using the squares of the major and minor axes as diagonal elements. Through matrix multiplication, we sequentially multiply the rotation matrix, the diagonal matrix, and the transpose of the rotation matrix to obtain the omnidirectional roughness tensor matrix. The main diagonal elements of this matrix represent the roughness components along the coordinate axes, while the off-diagonal elements represent the shear roughness components. For example, based on the obtained lengths of the major and minor axes of the ellipse and the rotation angle, we first construct the corresponding rotation matrix and a diagonal matrix with the squares of the axis lengths as elements. Then, through continuous matrix multiplication transformations, we calculate the specific element values ​​at each position in the roughness tensor matrix, where the values ​​on the main diagonal represent the lateral and longitudinal roughness components, respectively.

[0027] Please see Figure 4 The specific steps for obtaining the updated hydrodynamic parameters are as follows: S311: Establish fluid control equations based on the physical spatial properties of the mountain topographic grid region, set the current simulation time step, perform numerical discretization to solve the fluid control equations, obtain the velocity component values ​​of fluid particles at the current time and perform vector synthesis to obtain the water flow velocity vector. A two-dimensional shallow water equation set is used as the core model, which includes continuity and momentum equations. The current simulation time step is set to 0.1 seconds, calculated based on the Courant-Friedrichs-Lewy (CFL) stability condition to ensure that the information propagation speed during numerical solution is less than the mesh generation speed. Numerical discretization is performed on the fluid control equations, and the computational domain is divided into regular control volumes using the finite volume method. The flux through each control volume interface is calculated. Within the current time step, based on the fluid state at the previous moment, the momentum equation is iteratively solved to obtain the velocity components of the fluid particles along the horizontal and vertical axes of the planar coordinate system. Vector synthesis is performed on these two components. The arithmetic square root of the sum of the squares of the velocity components is used to obtain the velocity magnitude. The ratio of the vertical to horizontal velocity components is calculated using inverse trigonometric functions to obtain the velocity direction angle, thus yielding the complete flow velocity vector. For example, when calculating the flow velocity at a certain grid point, first solve for the velocity components of that point on the horizontal and vertical axes respectively, then use the Pythagorean theorem to synthesize these two components to obtain the total flow velocity, and use the arctangent function to calculate the deflection angle of the synthesized velocity vector relative to the horizontal axis.

[0028] S312: Project the water flow velocity vector onto the principal axis coordinate system defined by the roughness tensor matrix, calculate the roughness fractal dimension value matching the current water flow direction, quantify the surface friction effect of the mountain topographic grid area under the current specified flow direction based on the roughness fractal dimension value, and calculate the surface runoff resistance coefficient. First, the angle difference between the direction angle of the water flow velocity vector and the rotation angle of the principal axis of the roughness tensor matrix is ​​calculated. The roughness fractal dimension value matching the current water flow direction is calculated. Using the elliptic polar coordinate equation, the angle difference is used as the independent variable and substituted into the ellipse radius calculation logic with the major and minor axis lengths as parameters to obtain the polar radius value corresponding to that angle. This polar radius value is the equivalent roughness fractal dimension under the current flow direction. Based on the roughness fractal dimension value, the surface friction effect of the mountain topographic grid area under the current specified flow direction is quantified. Through a preset empirical conversion relationship, the fractal dimension value is converted into the Manning roughness coefficient. This conversion relationship is based on: measuring the resistance of surfaces with different fractal dimensions through laboratory flume experiments and fitting an exponential function relationship between the fractal dimension and the Manning coefficient. Finally, the surface runoff resistance coefficient is calculated. For example, based on the difference between the current flow direction angle and the angle of the principal axis of the ellipse, the fractal dimension in that direction is calculated by substituting it into the ellipse equation. Then, according to the preset conversion formula, the fractal dimension is mapped to a specific Manning roughness coefficient, and the drag coefficient value used for the hydrodynamic equation is further calculated.

[0029] S313: The fluid control equations are corrected using the surface runoff resistance coefficient. The fluid state variables in the fluid control equations are recalculated according to the current simulation time step to obtain the updated hydrodynamic parameters. The calculated surface runoff drag coefficient is substituted into the source term of the momentum equation as a drag term to consume fluid momentum. Based on the current simulation time step of 0.1 seconds, an explicit time integration scheme is used to recalculate the fluid state variables in the fluid control equations, including the water depth and velocity at the next moment. Specifically, the current flow rate is added to the net flux, then the drag loss is subtracted, and multiplied by the time step to update the momentum and mass within the fluid control volume, yielding the updated hydrodynamic parameters. For example, when updating the fluid state at a certain moment, the calculated drag coefficient is used to derive the drag term. The momentum loss caused by drag is subtracted from the current momentum, and the velocity and water depth at that grid point at the next moment are calculated using the time step, reflecting the deceleration effect of drag on the water flow.

[0030] Please see Figure 6 The specific steps for obtaining the environmental biological stress intensity index are as follows: S411: Identify the vegetation species type identifiers growing in the mountain terrain grid area within the current simulation time step, and simultaneously scan the neighboring vegetation species type identifiers growing in the neighborhood of the mountain terrain grid area, and set the interspecific competition coefficient between the neighboring identifiers and the central area identifiers. A vegetation classification model based on convolutional neural networks is used to process remote sensing imagery or field survey data. The model consists of an input layer, three convolutional layers, two pooling layers, and a fully connected output layer. The input layer receives normalized multispectral image data; the first to third convolutional layers use 32, 64, and 128 3x3 convolutional kernels, respectively, to extract progressively abstract texture features, with modified linear units as the activation function; the pooling layers use a 2x2 max pooling strategy to reduce dimensionality; the output layer outputs the probability of each species using the Softmax function, selecting the species number with the highest probability as the vegetation species type identifier. Simultaneously, the model scans the neighboring vegetation species type identifiers within a 3x3 grid area of ​​the mountain terrain. An interspecific competition coefficient is defined between the neighboring identifier and the central region identifier. This coefficient is derived from a pre-set ecological interaction matrix, whose values ​​are based on long-term ecological observation data statistics, calculating the biomass inhibition rate when different species coexist. For example, when processing the central grid, its vegetation type is identified, and the vegetation types of its surrounding neighboring grids are traversed. For each pair of "center-neighbor" species combinations, the corresponding competition coefficient value is found in the ecological interaction matrix.

[0031] S412: Calculate the Euclidean distance between the geometric center of the neighboring mountain terrain grid region and the geometric center of the current mountain terrain grid region to obtain the spatial Euclidean distance; Extract the x and y coordinates of the centers of two grid cells, calculate the square of the difference between the x and y coordinates, add them together, and then take the square root to obtain the spatial Euclidean distance. For example, extract the center coordinates of the central grid and the neighboring grid cells, calculate the difference between their coordinates in the horizontal and vertical directions, and then take the square root of the sum of the squares of these two differences to obtain the straight-line distance between the center points of the two grid cells.

[0032] S413: The spatial Euclidean distance is used as a distance decay factor and the interspecific competition coefficient is weighted and summed to quantify the competitive pressure exerted by the neighboring vegetation on the central vegetation, thus obtaining the environmental biological stress intensity index. The calculation formula, which considers the spatial attenuation accumulation of vegetation cover abundance, is as follows: In the formula, The index represents the calculated environmental biological stress intensity, with a value range of 0 or higher. The higher the index, the greater the survival competition pressure on the vegetation. This represents the total number of grid cells identified as having valid vegetation within the neighborhood, obtained by traversing the neighborhood grid and counting non-empty vegetation cells; Representing the The interspecific competition coefficient between vegetation in each neighborhood unit and the central vegetation is derived from the ecological interaction matrix, with a value ranging from 0 to 1. Representing the The vegetation cover correction factor for each neighborhood unit is obtained directly by calculating the normalized vegetation index (NDVI) of that unit, with a value ranging from 0 to 1. Representing the The spatial Euclidean distance between the center of each neighboring unit and the center of the current unit is calculated using the distance formula between two points and is used as a distance decay factor. This is the distance smoothing coefficient, set to 1. This coefficient is determined based on the mathematical constraints in the spatial interaction model to prevent the denominator from being zero.

[0033] For example, if there are 3 valid vegetation units in the neighborhood, a correction factor is set for each. Smoothing coefficient The first unit competition coefficient ,distance Individual calculation is The second unit competition coefficient ,distance Individual calculation is The third unit competition coefficient ,distance Individual calculation is Add these three components together and calculate. The environmental biological stress intensity index was obtained. for .

[0034] Please see Figure 6 The specific steps for obtaining the digital twin simulation results of the ecological restoration of the mountain park are as follows: S511: Call the preset vegetation succession probability matrix, use the updated hydrodynamic parameters and the vegetation species type identifier in the mountain topography grid area as the joint index key value, search for the matching succession rule in the vegetation succession probability matrix, extract the probability value of vegetation evolving from the current state to the next succession stage, and obtain the theoretical state transition probability. The updated hydrodynamic parameters (such as flow velocity levels) and vegetation species type identifiers within the mountain topographic grid area are used as joint index keys. Continuous hydrodynamic parameters are discretized into level identifiers, with multiple possibilities: the first possibility is a flow velocity less than 1 m / s, corresponding to level 1; the second possibility is a flow velocity between 1 and 3 m / s, corresponding to level 2; and the third possibility is a flow velocity greater than 3 m / s, corresponding to level 3. Matching succession rules are retrieved from the vegetation succession probability matrix. This matrix stores the probability distribution of evolution to other species or maintaining the current state under specific flow velocity and current species conditions. This probability distribution is calculated based on Markov chain analysis of historical ecological succession case data. The probability values ​​of vegetation evolving from the current state to the next succession stage are extracted to obtain the theoretical state transition probability. For example, first, it is determined which preset flow velocity range level the current flow velocity value falls into. Then, combined with the vegetation species ID of the current grid, the specific row and column are located in the probability matrix, and the theoretical probability value of that species evolving into the target species (e.g., from herbaceous to shrub) under the current flow velocity conditions is read.

[0035] S512: Using the environmental biological stress intensity index as the probability penalty coefficient, calculate the difference between the theoretical state transition probability and the environmental biological stress intensity index to correct the succession probability, obtain the actual growth transition probability, compare it with the preset random judgment threshold, and determine the vegetation community succession state of the mountain terrain grid area in the next simulation time step. The actual growth transition probability is obtained by subtracting the environmental biological stress intensity index from the theoretical state transition probability. If the difference is less than zero, the actual growth transition probability is set to zero; if the difference is greater than one, the actual growth transition probability is set to one; otherwise, the difference is taken as the actual growth transition probability. Subsequently, a uniformly distributed random number between zero and one is generated and compared with the actual growth transition probability. If the generated random number is less than or equal to the actual growth transition probability, succession is determined to have occurred; otherwise, succession is determined not to have occurred, thus determining the vegetation community succession state of the mountain terrain grid area in the next simulation time step. For example, if the theoretical state transition probability of a certain grid under current hydrodynamic parameters is 80%, while the environmental biological stress intensity index, obtained by weighted summation of neighborhood distance decay and interspecific competition coefficient, is 30%, then subtracting the stress index from the theoretical probability yields 50%, with the difference between zero and one. Therefore, the actual growth transition probability is taken as 50%. Subsequently, a random decimal between zero and one is generated. If this random decimal does not exceed 50%, the vegetation in that grid is considered to have undergone succession; otherwise, the original state remains unchanged. When the theoretical state transition probability is low and the environmental biological stress intensity index is high, resulting in a negative difference, the actual growth transition probability is set to zero according to the boundary conditions. In this case, regardless of the generated random decimal between zero and one, succession is considered not to have occurred.

[0036] S513: Map the succession status of vegetation communities to the grid of the digital elevation model of mountain terrain, perform spatial fusion of vegetation dynamic evolution information and static geographic information of terrain, and generate digital twin simulation results of ecological restoration of mountain park. The attribute data of the corresponding grid is updated based on the determined succession state, including species code, vegetation height, and cover parameters. Spatial fusion is performed between the dynamic evolution information of vegetation and the static geographic information of the terrain. Using a geographic information visualization engine, the updated vegetation layer is overlaid on the 3D terrain model, and different species are assigned corresponding texture and color attributes. State snapshots are recorded at each time step, and these snapshots are played continuously to demonstrate the ecological restoration process, generating digital twin simulation results of the ecological restoration of the mountain park. For example, when a grid is determined to have undergone succession, its attribute data in the digital model is updated to the parameters of the new species. At the same time, in the visualization interface, the rendering material at that location is replaced from the texture and color representing the original species (such as light green grass) to the texture and color representing the new species (such as dark green shrubs), and this frame is saved as part of the simulation animation.

[0037] A digital twin system for ecological restoration of mountain parks, the system includes: The terrain construction module uses lidar mapping equipment to build a digital elevation model of the mountain terrain, divides the mountain terrain grid area in the model and extracts surface elevation data from multiple directions to generate a multi-directional elevation profile sequence. The roughness modeling module calculates the fractal dimension of micro-topography in a multi-azimuth elevation profile sequence, and performs fitting processing on the fractal dimension of micro-topography to construct a roughness tensor matrix. The hydrodynamic parameter correction module establishes the fluid control equations and obtains the water flow velocity vector. It then projects the water flow velocity vector onto the roughness tensor matrix to calculate the equivalent fractal dimension, determines the surface runoff resistance coefficient, updates the fluid control equations, and obtains the updated hydrodynamic parameters. The biological stress calculation module identifies the vegetation species type identifiers in the mountain topographic grid area and its neighboring areas. Based on the spatial Euclidean distance between the neighboring mountain topographic grid area and the central area, the spatial Euclidean distance is used as a distance decay factor to perform a weighted summation on the preset interspecific competition coefficient to obtain the environmental biological stress intensity index. The community succession simulation module matches the theoretical state transition probability with the updated hydrodynamic parameters, corrects the theoretical state transition probability according to the environmental biological stress intensity index, obtains the actual growth transition probability, determines the vegetation community succession state at the next moment, and generates digital twin simulation results for the ecological restoration of the mountain park.

[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 digital twin method for ecological restoration of mountain parks, characterized in that, Includes the following steps: S1: Construct a digital elevation model of the mountain terrain using lidar mapping equipment, divide the mountain terrain grid area in the model, extract surface elevation data from multiple directions, and generate a multi-directional elevation profile sequence. S2: Calculate the micro-topographic fractal dimension of the multi-azimuth elevation profile sequence, and fit the micro-topographic fractal dimension to construct the roughness tensor matrix; S3: Establish the fluid control equation and obtain the water flow velocity vector. Project the water flow velocity vector onto the roughness tensor matrix to calculate the equivalent fractal dimension, determine the surface runoff resistance coefficient, update the fluid control equation, and obtain the updated hydrodynamic parameters. S4: Identify vegetation species type identifiers within the mountain topographic grid area and its neighboring areas. Based on the spatial Euclidean distance between the neighboring mountain topographic grid area and the central area, use the spatial Euclidean distance as a distance decay factor to perform a weighted summation of the preset interspecific competition coefficients to obtain the environmental biological stress intensity index. S5: Based on the updated hydrodynamic parameters, match the theoretical state transition probability, correct the theoretical state transition probability according to the environmental biological stress intensity index, obtain the actual growth transition probability, determine the vegetation community succession state at the next moment, and generate digital twin simulation results for the ecological restoration of the mountain park.

2. The digital twin method for ecological restoration of mountain parks according to claim 1, characterized in that, The multi-azimuth elevation profile sequence includes discrete azimuth markers of grid computing units and a set of surface elevation values ​​distributed along the azimuth direction. The roughness tensor matrix includes the principal axis deflection angle of the micro-topographic fractal fitting ellipse, the fractal dimension components in the first principal axis direction, and the fractal dimension components in the second principal axis direction. The updated hydrodynamic parameters include the velocity vector corrected by the surface runoff resistance coefficient and the updated water depth values. The environmental biological stress intensity index is specifically a scalar value obtained by weighted summation of the interspecific competition inhibition coefficient based on the distance attenuation factor between the neighboring mountain topographic grid area and the central area. The digital twin simulation results of the mountain park ecological restoration include the vegetation community succession state of the grid area at the next moment and the spatial distribution mapping information of the state on the mountain topographic digital elevation model.

3. The digital twin method for ecological restoration of mountain parks according to claim 1, characterized in that, The specific steps for obtaining the multi-azimuth elevation profile sequence are as follows: S111: Collect point cloud data of the target area of ​​the mountain park, convert the three-dimensional spatial coordinate information in the point cloud data into a digital elevation model of the mountain terrain, perform regional grid segmentation processing on the digital elevation model of the mountain terrain, and obtain the mountain terrain grid area. S112: Locate the geometric center of the mountain terrain grid area, establish a plane coordinate system with the geometric center as the origin, set multiple ray directions on the plane with a preset angle interval, and obtain discrete azimuth directions; S113: Using the geometric center of the mountain terrain grid area as the extraction origin, extract the surface elevation data corresponding to the surface of the digital elevation model of the mountain terrain one by one along each of the discrete azimuth directions and arrange them according to spatial distance to generate a multi-azimuth elevation profile sequence.

4. The digital twin method for ecological restoration of mountain parks according to claim 3, characterized in that, The specific steps for obtaining the roughness tensor matrix are as follows: S211: Call the box-counting dimension algorithm, set grid boxes with multiple side lengths as the measurement scale, cover the surface elevation data curves in the multi-azimuth elevation profile sequence, count the number of non-empty boxes required to cover the curve at each box side length scale, establish a linear regression relationship between the logarithm of the reciprocal of the box side length and the logarithm of the number of boxes, perform an overall correlation analysis between the natural logarithm of the reciprocal of the box side length and the natural logarithm of the corresponding number of non-empty boxes, and combine the cooperative change characteristics of the two in all scale levels and the discrete distribution characteristics of the natural logarithm of the reciprocal of the box side length to extract the regression slope and calculate the micro-topographic fractal dimension; S212: For the fractal dimension of the micro-topography corresponding to each discrete azimuth direction, the least squares fitting algorithm is used to approximate the geometric shape of the ellipse, and the length of the major axis, the length of the minor axis, and the rotation angle of the major axis relative to the preset reference coordinate system are calculated to obtain the geometric feature parameters of the ellipse fitting. S213: Call the major axis length, minor axis length, and rotation angle from the geometric feature parameters of the ellipse fitting, define the component weights of the surface roughness in each direction of the plane coordinate system, and generate the roughness tensor matrix.

5. The digital twin method for ecological restoration of mountain parks according to claim 4, characterized in that, The specific steps for obtaining the updated hydrodynamic parameters are as follows: S311: Establish fluid control equations based on the physical spatial properties of the mountain topographic grid region, set the current simulation time step, perform numerical discretization to solve the fluid control equations, obtain the velocity component values ​​of fluid particles at the current time and perform vector synthesis to obtain the water flow velocity vector. S312: Project the water flow velocity vector onto the principal axis coordinate system defined by the roughness tensor matrix, calculate the roughness fractal dimension value matching the current water flow direction, quantify the surface friction effect of the mountain topographic grid area under the current specified flow direction based on the roughness fractal dimension value, and calculate the surface runoff resistance coefficient. S313: The fluid control equations are corrected using the surface runoff resistance coefficients. The fluid state variables in the fluid control equations are recalculated based on the current simulation time step to obtain the updated hydrodynamic parameters.

6. The digital twin method for ecological restoration of mountain parks according to claim 5, characterized in that, The specific steps for obtaining the environmental biological stress intensity index are as follows: S411: Identify the vegetation species type identifiers growing in the mountain terrain grid area within the current simulation time step, and simultaneously scan the neighboring vegetation species type identifiers growing in the neighborhood of the mountain terrain grid area, and set the interspecific competition coefficient between the neighboring identifiers and the central area identifiers. S412: Calculate the Euclidean distance between the geometric center of the neighboring mountain terrain grid region and the geometric center of the current mountain terrain grid region to obtain the spatial Euclidean distance; S413: The spatial Euclidean distance is used as a distance decay factor and the interspecific competition coefficient is weighted and summed to quantify the competitive pressure exerted by the neighboring vegetation on the central vegetation, thus obtaining the environmental biological stress intensity index.

7. The digital twin method for ecological restoration of mountain parks according to claim 6, characterized in that, The specific steps for obtaining the digital twin simulation results of the ecological restoration of the mountain park are as follows: S511: Call the preset vegetation succession probability matrix, use the updated hydrodynamic parameters and the vegetation species type identifier in the mountain topography grid area as the joint index key value, search for the matching succession rule in the vegetation succession probability matrix, extract the probability value of vegetation evolving from the current state to the next succession stage, and obtain the theoretical state transition probability. S512: Using the environmental biological stress intensity index as a probability penalty coefficient, calculate the difference between the theoretical state transition probability and the environmental biological stress intensity index to correct the succession probability, obtain the actual growth transition probability, compare it with the preset random judgment threshold, and determine the vegetation community succession state of the mountain terrain grid area in the next simulation time step. S513: Map the vegetation community succession state to the grid of the digital elevation model of the mountain terrain, perform spatial fusion of vegetation dynamic evolution information and static geographical information of terrain, and generate digital twin simulation results of ecological restoration of the mountain park.

8. The digital twin method for ecological restoration of mountain parks according to claim 3, characterized in that, The digital elevation model of the mountain terrain includes a two-dimensional raster matrix covering the target area of ​​the mountain park, constructed based on point cloud data, and the absolute elevation value of the ground surface at each raster position in the two-dimensional raster matrix. The absolute elevation value of the ground surface is calculated by spatial interpolation based on the three-dimensional spatial coordinate information in the point cloud data.

9. A digital twin system for ecological restoration of mountain parks, characterized in that, The digital twin method for ecological restoration of mountain parks according to any one of claims 1-8, wherein the system comprises: The terrain construction module uses lidar mapping equipment to build a digital elevation model of the mountain terrain, divides the mountain terrain grid area in the model and extracts surface elevation data from multiple directions to generate a multi-directional elevation profile sequence. The roughness modeling module calculates the micro-topographic fractal dimension of the multi-azimuth elevation profile sequence, and performs fitting processing on the micro-topographic fractal dimension to construct a roughness tensor matrix. The hydrodynamic parameter correction module establishes the fluid control equation and obtains the water flow velocity vector. It projects the water flow velocity vector onto the roughness tensor matrix to calculate the equivalent fractal dimension, determines the surface runoff resistance coefficient, and updates the fluid control equation to obtain the updated hydrodynamic parameters. The biological stress calculation module identifies the vegetation species type identifiers in the mountain topographic grid area and its neighboring areas. Based on the spatial Euclidean distance between the neighboring mountain topographic grid area and the central area, the spatial Euclidean distance is used as a distance decay factor to perform a weighted summation on the preset interspecific competition coefficient to obtain the environmental biological stress intensity index. The community succession simulation module matches the theoretical state transition probability with the updated hydrodynamic parameters, corrects the theoretical state transition probability according to the environmental biological stress intensity index, obtains the actual growth transition probability, determines the vegetation community succession state at the next moment, and generates digital twin simulation results for the ecological restoration of the mountain park.

Images