Garden landscape intelligent planning method based on digital twinning

Abstract

The application relates to the technical field of data processing, and discloses a garden landscape intelligent planning method based on digital twinning. The method acquires three-dimensional geological profile data of a planning area, establishes an underground soil porosity and saturation tensor field, and divides the underground space into a multi-scale hexahedral microenvironment grid; root system morphological parameters of a candidate plant community are acquired, and a plant physiological feature library containing a root system water absorption rate curve and a nutrient absorption spectrum is constructed; a spatial mapping topological relationship between an aboveground plant three-dimensional canopy model and the underground grid is established, and water potential gradient changes of the microenvironment grid under different aboveground layout schemes are calculated. The application quantitatively constrains the surface layout and the dynamic distribution of underground soil water and fertilizer, solves the problem that the surface layout is disconnected from the underground ecological carrying capacity in the prior art, leading to plant decline, and improves the ecological stability of the digital twinning planning scheme.

Description

Technical Field

[0001] This invention relates to the field of data processing technology and discloses a smart planning method for garden landscape based on digital twins. Background Technology

[0002] Existing digital twin-based landscape planning systems typically acquire surface elevation data and remote sensing imagery to construct a three-dimensional geometric model of the land surface. During the planning process, the system places different types of plants as three-dimensional geometric shapes within the surface model. By calculating the geometric spacing between plants, visual occlusion, and the distribution of sunlight and shadow on the surface, the system adjusts the spatial positions of the plants. This planning process involves geometric layout calculations entirely within the three-dimensional space above the ground. By continuously iterating the spatial coordinates of the plants, the system achieves the predetermined green coverage rate and spatial hierarchy in the visual presentation of the landscape.

[0003] Based on the specific implementation methods of the aforementioned existing technologies, the core technical problem lies in the fact that existing digital twin planning methods only arrange plants geometrically above the ground, completely severing the physical interaction between the surface plant layout and the underground soil microenvironment. This results in a disconnect between the planned plant community layout and the actual ecological carrying capacity of the underground soil. Because the root systems of different plants exhibit complex spatial interweaving underground and vary in their ability to absorb soil moisture and nutrients, relying solely on the geometric spacing at the surface cannot reflect the true physical state of underground water transport and nutrient consumption. Consequently, after actual planting, the generated digital twin planning scheme is prone to localized plant decline or even death due to imbalances in water and fertilizer resource competition among underground roots, failing to guarantee the long-term physical stability of the entire garden landscape ecosystem. Summary of the Invention

[0004] The purpose of this invention is to provide a digital twin-based intelligent planning method for garden landscapes, which can effectively solve the problems mentioned in the background art.

[0005] To achieve the above objectives, the technical solution adopted by the present invention is as follows: A digital twin-based intelligent planning method for garden landscape includes: acquiring three-dimensional geological profile data of the planning area, establishing a tensor field of underground soil porosity and saturation based on the three-dimensional geological profile data, and dividing the underground space into a multi-scale hexahedral microenvironment grid in the digital twin system. Root morphology parameters of candidate plant communities were obtained, and a plant physiological feature library containing root water absorption rate curves and nutrient uptake spectra was constructed. A spatial mapping topological relationship between the three-dimensional canopy model of aboveground plants and the multi-scale hexahedral microenvironment mesh is established in the digital twin system; Based on the spatial mapping topology and the tensor fields of underground soil porosity and saturation, the water potential gradient changes of the multi-scale hexahedral microenvironment grid under different aboveground plant spatial layout schemes are calculated. A graph neural network model is constructed with the optimization objectives of minimizing the water potential gradient range and minimizing the area of ​​nutrient-deficient regions in the multi-scale hexahedral microenvironment grid. The graph neural network model is used to aggregate node features and update edge weights. The optimal spatial combination configuration of the aboveground plant community is solved by backpropagation gradient descent. The optimal spatial combination configuration is then mapped to a digital twin surface layer to complete the landscape planning.

[0006] Preferably, the three-dimensional geological profile data of the planning area is obtained, and a tensor field of underground soil porosity and saturation is established based on the three-dimensional geological profile data. The underground space is divided into a multi-scale hexahedral micro-environment grid in the digital twin system, including: obtaining electromagnetic wave reflection signals of the planning area through ground penetrating radar, performing time-frequency domain transformation on the electromagnetic wave reflection signals to extract the dielectric constants of strata at different depths, and calculating the initial soil porosity distribution based on the dielectric constants. Multi-point saturation data collected by a soil moisture sensor network are acquired, and the underground soil porosity and saturation tensor field is constructed by combining the initial soil porosity distribution. The underground space is initially divided using an octree structure. The porosity variance within each initially divided grid is calculated. When the porosity variance is greater than a preset geological bedding threshold, the initially divided grid is subdivided into eight equal parts until the porosity variance of all grids is less than the preset geological bedding threshold, thereby generating the multi-scale hexahedral microenvironment grid.

[0007] Preferably, the root morphology parameters of the candidate plant community are obtained, and a plant physiological feature library containing root water absorption rate curves and nutrient absorption spectra is constructed. This includes: performing a three-dimensional tomographic scan of the root system of the candidate plant community, extracting the root spatial orientation coordinate set and the radius dimensions of each root level, substituting the coordinate set and the radius dimensions into the modified van der Hoff soil hydrodynamic equation, inputting different soil water potential boundary conditions to solve, and generating the root water absorption rate curve with soil water potential as the independent variable and water flux as the dependent variable. Ion adsorption kinetic parameters of nitrogen, phosphorus and potassium elements in the candidate plant communities are extracted, and the nutrient absorption spectrum is generated by fitting it along the time dimension. The root water absorption rate curve and the nutrient absorption spectrum are associated with the plant species number and stored in the plant physiological feature database.

[0008] Preferably, establishing a spatial mapping topological relationship between the three-dimensional canopy model of aboveground plants and the multi-scale hexahedral microenvironment mesh in the digital twin system includes: obtaining the vertical projection polygon of the three-dimensional canopy model of aboveground plants on the ground surface plane, and determining the multi-scale hexahedral microenvironment mesh within the coverage area of ​​the vertical projection polygon as candidate associated meshes; Construct the root envelope equation of the candidate plant community, perform spatial Boolean intersection operation between the root envelope equation and the candidate associated grid, and extract the set of intersecting grids as the actual associated grids; Calculate the Euclidean distance from the center point of the three-dimensional canopy model of the aboveground plants to the center point of each grid in the substantially associated grid, and construct a spatial mapping topological adjacency matrix describing the hydraulic connectivity strength between the aboveground plants and the underground grid using the reciprocal of the Euclidean distance as the weight.

[0009] Preferably, based on the spatial mapping topology and the tensor field of underground soil porosity and saturation, the water potential gradient change of the multi-scale hexahedral microenvironment grid under different aboveground plant spatial layout schemes is calculated, including: for any aboveground plant spatial layout scheme, retrieving the root water absorption rate curve of the corresponding plant from the plant physiological feature database as the source and sink term; Using the underground soil porosity and saturation tensor fields as initial boundary conditions, and combining them with the spatial mapping topological adjacency matrix, Darcy's law and Richards' equations are substituted to construct the underground three-dimensional transient seepage partial differential equation. The partial differential equation of the underground three-dimensional transient seepage is discretized using the finite volume method, and iteratively solved within a set time step. The water potential value matrix of the multi-scale hexahedral microenvironment grid at each time node is output, and the water potential gradient change is calculated based on the water potential values ​​of the adjacent grids.

[0010] Preferably, a graph neural network model is constructed with the optimization objectives of minimizing the water potential gradient range and minimizing the area of ​​nutrient deficit regions in the multi-scale hexahedral microenvironment grid. The model includes: using the multi-scale hexahedral microenvironment grid as graph nodes, and using the water potential matrix and soil nutrient concentration as the initial feature vectors of the graph nodes. The spatial mapping topological adjacency matrix is ​​used as the edge connection relationship and initial edge weight of the graph neural network model; Construct a graph neural network model containing three layers of graph convolutional network, and update the node state by weighted aggregation of the feature vectors of adjacent graph nodes based on the initial edge weights; A joint loss function is constructed, in which the difference between the maximum and minimum water potential values ​​in the graph node state vector is taken as the water potential gradient range term, and the number of graph nodes in the state vector with nutrient concentrations below the minimum threshold for plant growth is taken as the nutrient deficit area term. The two terms are weighted and summed to generate the joint loss function.

[0011] Preferably, after generating the multi-scale hexahedral microenvironment mesh, the method further includes: extracting the porosity difference between adjacent meshes in the multi-scale hexahedral microenvironment mesh, and calculating the distribution gradient of the porosity difference in three-dimensional space; When the porosity difference between two adjacent grids is greater than a preset abrupt change threshold, it is determined that there is a geological fault or rock stratum boundary between the two grids, and a transition grid is inserted between the two grids. The node coordinates of the transition mesh are smoothed using Laplacian to eliminate mesh distortion, and the hydraulic parameters of the transition mesh are set to the interpolation calculation results of the parameters of two adjacent meshes. In subsequent calculations of the water potential gradient change, the transition grid is used as an independent calculation node to participate in the Darcy's law solution.

[0012] Preferably, after generating the root water absorption rate curve, the method further includes: acquiring historical meteorological data and real-time microclimate monitoring data of the planning area, and extracting saturated water vapor pressure difference and solar radiation flux from the historical meteorological data and the real-time microclimate monitoring data. A stomatal conductance model of the canopy is established, and the saturated water vapor pressure difference and the solar radiation flux are input into the stomatal conductance model of the canopy to calculate the canopy transpiration rate; The canopy transpiration rate is superimposed as a boundary constraint on the water flux, the dependent variable of the root water uptake rate curve, to generate a canopy-root coupled water uptake rate curve. When calculating the water potential gradient change, the canopy-root coupled water absorption rate curve is used to replace the root water absorption rate curve as the source and sink term.

[0013] Preferably, after constructing a spatial mapping topological adjacency matrix describing the hydraulic connectivity strength between aboveground plants and underground grids, the method further includes: obtaining data on the types and concentrations of secondary metabolites secreted by the candidate plant communities, and constructing a plant allelopathy feature vector. Calculate the cosine similarity of the plant allelopathic feature vectors between different plant species. When the cosine similarity is lower than a preset antagonism threshold, it is determined that the two corresponding plant species have root allelopathic antagonism. In the spatial mapping topological adjacency matrix, the edge weights between the substantially associated grids corresponding to the two plants with root allelopathic antagonism are multiplied by a preset penalty factor for attenuation processing to generate an allelopathic constraint topological adjacency matrix. In the graph neural network model, the chemical constraint topological adjacency matrix is ​​used to replace the spatial mapping topological adjacency matrix for feature aggregation.

[0014] Preferably, after constructing the joint loss function, the method further includes: extracting the predicted water potential and predicted nutrient concentration of each graph node output by the graph neural network model during the process of solving the optimal spatial combination configuration of the aboveground plant community using backpropagation gradient descent; Based on the law of conservation of mass, calculate the residual between the sum of the flux flowing into each graph node and the sum of the flux flowing out of the graph node; When the residual is greater than the preset physical tolerance, the physical constraint penalty term is activated, and the square of the residual is multiplied by the preset penalty coefficient and added to the joint loss function; When calculating the gradient during backpropagation, the gradient propagation path outside the physical constraint penalty term is truncated, and the backpropagation gradient descent is forced to update the network weight parameters of the graph neural network model along the direction of reducing the residual, until the residual converges to zero.

[0015] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. This invention establishes a tensor field of underground soil porosity and saturation by acquiring three-dimensional geological profile data and dividing it into multi-scale hexahedral microenvironment grids. A spatial mapping topological relationship between the three-dimensional canopy model of aboveground plants and the underground grid is established within a digital twin system, extending landscape planning from surface geometric layout to underground physical field calculation. The scheme calculates the water potential gradient changes under different layouts based on Darcy's law and Richards' equations, constructing a graph neural network model to solve the problem with the objectives of minimizing the water potential gradient range and minimizing the area of ​​nutrient-deficient regions. This technique directly establishes quantitative constraints between the layout of surface plants and the dynamic distribution of water and fertilizer in underground soil, ensuring that the output spatial configuration conforms to the natural laws of groundwater and soil resource balance. It solves the technical problem of plant decline caused by the disconnect between surface layout and underground ecological carrying capacity in existing technologies, improving the ecological stability of digital twin planning schemes in real physical environments.

[0016] 2. By introducing porosity variance constraints into the underground grid division, performing octree subdivision, and inserting Laplace smoothed transition grids, the numerical error caused by grid distortion at geological faults in water potential calculations was reduced. By constructing a canopy stomatal conductance model, meteorological factors were converted into transpiration rates and superimposed on the root water absorption rate curve, making the source-sink boundary conditions of underground seepage calculations consistent with actual atmospheric physical exchange processes. By calculating the cosine similarity of plant secondary metabolite feature vectors to determine allelopathic antagonism and penalizing the edge weights of the topological adjacency matrix, the adjacent layout of plants with antagonistic effects in the underground grid was excluded. By introducing a flux residual penalty term based on mass conservation into the graph neural network loss function and truncating the gradient propagation path, the model parameter updates were forced to follow physical conservation laws, avoiding the output of abnormal layout coordinates that could not be implemented due to deviations from real physical laws during the neural network optimization process. Attached Figure Description

[0017] Figure 1 This is a flowchart illustrating the overall process of the intelligent landscape planning method of the present invention. Figure 2 This is a flowchart illustrating the construction process of the underground multi-scale hexahedral microenvironment mesh according to the present invention. Figure 3 This is a flowchart illustrating the construction process of the plant physiological characteristic library of this invention. Figure 4 This is a flowchart of the spatial mapping topology construction and allelopathic constraint processing of the present invention; Figure 5 This is a flowchart illustrating the calculation of water potential gradient changes and meteorological coupling in this invention. Figure 6 This is a flowchart of the graph neural network model construction and physical constraint solution of the present invention. Detailed Implementation

[0018] Please refer to the attached document. Figure 1 This embodiment provides a method for acquiring three-dimensional geological profile data of a planning area. Based on the three-dimensional geological profile data, a tensor field for underground soil porosity and saturation is established. In a digital twin system, the underground space is divided into a multi-scale hexahedral micro-environment grid. Electromagnetic wave reflection signals of the planning area are acquired using ground-penetrating radar. The electromagnetic wave reflection signals are then subjected to time-frequency domain transformation to extract the dielectric constants of strata at different depths. Based on the dielectric constants, the initial soil porosity distribution is calculated. The time-frequency domain transformation employs short-time Fourier transform processing, segmenting the electromagnetic wave reflection signals into windowed segments with a fixed time window. Fourier transform is performed on each segment to obtain the corresponding spectral information. Based on the correspondence between the propagation speed of electromagnetic waves in different media and the dielectric constant, the two-way travel time in the spectral information is converted into stratum depth information to obtain the dielectric constant values ​​corresponding to different depth locations.

[0019] Specifically, the relationship between dielectric constant and soil volumetric water content is calculated using the following formula:

[0020] in, Soil volumetric water content, This is the measured value of the formation dielectric constant.

[0021] The soil porosity value is equal to the soil saturation water content value. The initial soil porosity distribution at different depths is calculated using the dielectric constant saturation values ​​of strata at different depths. Multi-point saturation data collected by a soil moisture sensor network are acquired and combined with the initial soil porosity distribution to construct a tensor field for underground soil porosity and saturation. Specifically, the saturation data is the ratio of the actual volumetric water content of the soil to the corresponding soil porosity. Using the collected multi-point saturation data as measured control points, Kriging interpolation is used to interpolate the three-dimensional space corresponding to the initial soil porosity distribution, obtaining continuous porosity and saturation distribution data covering the entire planned area's underground space. Due to the anisotropy of the soil medium in three-dimensional space, the conduction characteristics of porosity and saturation differ in different spatial directions. Therefore, porosity and saturation are constructed as three-dimensional second-order symmetric tensors, forming a tensor field for underground soil porosity and saturation. The tensor elements at each spatial location in the tensor field correspond to the porosity and saturation conduction parameters at that location in different spatial directions.

[0022] Specifically, the formula for calculating the second-order saturation tensor is:

[0023] in, Spatial coordinates The second-order tensor of saturation at that point, , , The location is respectively at Saturation values ​​along the three principal axes. , , For the cross-coupling terms between different principal axes, the tensor has a symmetric structure and satisfies = , = , = The structure of the porosity tensor is the same as that of the saturation tensor, and the porosity tensor element at a corresponding position is the porosity value at that position in the corresponding direction.

[0024] Please refer to the attached document. Figure 2In the digital twin system, the underground space is divided into a multi-scale hexahedral microenvironment grid. Using the surface boundary of the planned area as the top surface and the preset maximum root distribution depth as the bottom surface, a three-dimensional cuboid computational domain of the underground space is constructed, and this domain is then meshed. An octree structure is used for initial subdivision of the underground space. The porosity variance within each initial subdivision grid is calculated. When the porosity variance exceeds a preset geological bedding threshold, the initial subdivision grid is further subdivided into eight equal parts until the porosity variance of all grids is less than the preset geological bedding threshold, generating a multi-scale hexahedral microenvironment grid. Each hexahedral microenvironment grid corresponds to a unique spatial coordinate range. The porosity tensor and saturation tensor within the grid are taken as the average tensor value within that grid's spatial range, serving as the hydraulic calculation parameters for that grid.

[0025] Please refer to the attached document. Figure 3 The root morphology parameters of candidate plant communities were obtained, and a plant physiological feature library including root water absorption rate curves and nutrient absorption spectra was constructed. Three-dimensional tomographic scanning was performed on the root systems of individual plants in the candidate plant communities to extract the spatial orientation coordinate set and root radius dimensions at each level. The spatial orientation coordinate set contains the three-dimensional spatial coordinates of all root nodes, and the root radius dimensions at each level correspond to the cross-sectional radii of root segments at different levels. Root segments are classified according to the hierarchy of taproot, primary lateral roots, secondary lateral roots, and tertiary lateral roots, with each root segment corresponding to a unique coordinate interval and radius dimension.

[0026] Substituting the coordinate set and radius dimensions into the modified van der Hoff soil hydrodynamic equation, and inputting different soil water potential boundary conditions, the solution is obtained, generating root water uptake rate curves with soil water potential as the independent variable and water flux as the dependent variable. The modified van der Hoff soil hydrodynamic equation is as follows:

[0027] in, Water flux absorbed by the root system per unit root length. To correspond to the root hydraulic conductivity under soil water potential, For rhizosphere soil water potential, For the water potential of the root xylem, This refers to the osmotic pressure difference between the rhizosphere soil solution and the root cell sap. This represents the length of the corresponding root segment.

[0028] Root hydraulic conductivity is positively correlated with root segment radius. Different root segments correspond to different initial values ​​of root hydraulic conductivity. The range of rhizosphere soil water potential covers the critical water potential range for plant survival. Different soil water potential boundary conditions are set with a fixed step size. The water flux value corresponding to each soil water potential value is obtained by substituting into the equation. The discrete data points obtained by the solution are fitted by cubic spline interpolation to generate a continuous root water absorption rate curve. The independent variable of the curve is soil water potential, and the dependent variable is the root water absorption flux per unit root length.

[0029] Ion adsorption kinetic parameters for nitrogen, phosphorus, and potassium in candidate plant communities were extracted and fitted over time to generate nutrient uptake spectra. The ion adsorption kinetic parameters include the maximum uptake rate and the Michaelis constant. The Michaelis equation was used to fit the relationship between the uptake rate of different nutrient ions and environmental concentration. The Michaelis equation is as follows:

[0030] in, To correspond to the absorption rate of nutrient ions, This represents the maximum absorption rate of the nutrient ion. The concentration of this nutrient ion in the rhizosphere environment. Let Michaelis-Menten constant be the constant of the nutrient ion. These are the identifiers for the types of nutrient molecules, corresponding to the three elements: nitrogen, phosphorus, and potassium.

[0031] For each candidate plant, the maximum absorption rates and Michaelis constants for nitrogen, phosphorus, and potassium were measured. Different nutrient concentration variation scenarios were set along the time dimension, and nutrient absorption spectra for the corresponding plants were generated. The horizontal axis of the nutrient absorption spectrum represents time, and the vertical axis represents the absorption rate of the corresponding nutrient ions, reflecting the absorption characteristics of different nutrients by plants at different time dimensions. The root water absorption rate curve and the nutrient absorption spectrum were associated with the plant species ID and stored in a plant physiological feature database. Each data point in the database corresponds to a unique plant species ID and contains complete data on the plant's root morphology parameters, root water absorption rate curve, and nutrient absorption spectrum.

[0032] Table 1. Basic physiological characteristics of candidate plant species

[0033] As shown in Table 1, the plant species numbers correspond one-to-one with the stored numbers in the plant physiological characteristic database. The average radius of the taproot and first-order lateral roots is the statistical average value of the root morphology parameters of the corresponding plants. The initial value of the root hydraulic conductivity is the measured value of the root hydraulic conductivity when the soil water potential is 0 MPa. The maximum absorption rates of nitrogen, phosphorus and potassium elements are the measured values ​​of the corresponding plant nutrient absorption kinetic parameters, providing basic data support for the construction of root water absorption rate curves and nutrient absorption spectra.

[0034] Please refer to the attached document. Figure 4 In a digital twin system, a spatial mapping topological relationship is established between a three-dimensional canopy model of aboveground plants and a multi-scale hexahedral microenvironment mesh. The three-dimensional canopy model of aboveground plants is a standardized three-dimensional model corresponding to the plant species. The model contains complete parameters of the three-dimensional geometry of the plant canopy, canopy height, and canopy width. The bottom center point of the model corresponds to the coordinates of the plant's planting position on the ground surface. The vertical projection polygon of the three-dimensional canopy model of aboveground plants on the ground surface plane is obtained. The vertical projection polygon is a closed convex polygon formed by the vertical projection points of all vertices of the canopy model on the xy plane of the ground surface. The multi-scale hexahedral microenvironment mesh within the coverage area of ​​the vertical projection polygon is determined as a candidate associated mesh. The top center point of the candidate associated mesh is located inside or on the boundary of the vertical projection polygon.

[0035] The root envelope equation for candidate plant communities is constructed. A spatial Boolean intersection operation is performed between the root envelope equation and the candidate associated grids, and the intersecting grid set is extracted as the substantially associated grid. The root envelope equation adopts an ellipsoidal equation form, specifically as follows:

[0036] in,( , , () represents the coordinates of the center point of the root envelope, where the center point is... , The coordinates are consistent with the coordinates of the bottom center point of the corresponding three-dimensional canopy model of aboveground plants. The coordinates are the coordinates of the midpoint of the root system's depth. , For the root envelope surface in , The length of the horizontal semi-axis in the direction corresponds to the radius of the plant's crown. For the root envelope surface in The length of the vertical semi-axis in the direction corresponds to the maximum root distribution depth of the plant.

[0037] Please refer to the attached document. Figure 2 The root envelope equation defines the three-dimensional spatial distribution range of the corresponding plant roots. A spatial Boolean intersection operation is performed between the root envelope equation and the hexahedral boundary equation of each candidate associated grid. If any part of the candidate associated grid lies within the root envelope, the grid is considered to have an intersection relationship with the root envelope and is included in the substantially associated grid set. The substantially associated grid set corresponds to the soil grid range actually affected by the plant roots in the underground space, providing spatial constraints for the calculation of hydraulic connectivity between aboveground plants and the underground soil microenvironment.

[0038] Calculate the Euclidean distance from the center point of the 3D canopy model of aboveground plants to the center point of each grid in the substantially associated grid. Using the reciprocal of the Euclidean distance as weights, construct a spatial mapping topological adjacency matrix describing the hydraulic connectivity between aboveground plants and the subsurface grid. The weights are calculated as follows:

[0039] in, For the first The above-ground plants and the first Hydraulic connectivity weights between underground hexahedral microenvironment grids For the first The center point of each above-ground plant and the first The Euclidean distance between the center points of each underground grid This is a local constant used to avoid cases where the denominator is 0.

[0040] The spatial mapping topological adjacency matrix has rows corresponding to the number of aboveground plants and columns corresponding to the number of subsurface hexahedral microenvironment grids. The elements in the matrix represent the hydraulic connectivity weights between the aboveground plants in the corresponding row and the subsurface grids in the corresponding column. When a subsurface grid does not belong to the substantially associated grid set of its corresponding aboveground plant, the element at that position in the matrix has a value of 0. The spatial mapping topological adjacency matrix comprehensively describes the spatial mapping topological relationship between aboveground plant communities and subsurface multi-scale hexahedral microenvironment grids, providing topological constraints for subsequent subsurface seepage calculations and graph neural network feature aggregation.

[0041] Please refer to the attached document. Figure 5 Based on spatial mapping topological relationships and the tensor fields of subsurface soil porosity and saturation, this study calculates the water potential gradient changes in a multi-scale hexahedral microenvironment grid under different aboveground plant spatial layout schemes. Specifically, for any aboveground plant spatial layout scheme, the scheme includes the planting location coordinates, planting quantity, and species configuration for each plant. The root water absorption rate curves of the corresponding plants are retrieved from a plant physiological characteristic database as source and sink terms, which describe the water consumption caused by root water absorption in the subsurface seepage field. Using the subsurface soil porosity and saturation tensor fields as initial boundary conditions, and combining them with the spatial mapping topological adjacency matrix, Darcy's law and Richards equations are substituted to construct the subsurface three-dimensional transient seepage partial differential equations. The Richards equations for three-dimensional transient seepage are:

[0042] in, Soil water content is given by , and soil volumetric water content is given by , which represents the effect of soil volumetric water content on matrix potential. The derivative, For soil matrix potential, For time, This represents the unsaturated soil hydraulic conductivity tensor, which is related to the soil porosity tensor, saturation tensor, and matrix potential. The coordinates are in the vertical direction, with upward being positive. This is the root water source and sink term, which is related to the root water absorption rate curve of the corresponding plant and the weight of the spatial mapping topological adjacency matrix.

[0043] The unsaturated soil hydraulic conductivity tensor was calculated using the van Genuchten model based on the soil porosity tensor, saturation tensor, and matrix potential. The root water uptake source-sink term was taken as the sum of the products of the root water flux of all plants associated with the corresponding grid and the corresponding hydraulic connectivity weight. The partial differential equations of subsurface three-dimensional transient seepage were discretized using the finite volume method. Each hexahedral microenvironment grid was treated as an independent control volume, and flux integration was performed over the six faces of the control volume to transform the partial differential equations into a system of linear algebraic equations.

[0044] The solution is iteratively solved within a set time step. The convergence criterion is that the maximum difference between the calculated water potential values ​​of two adjacent iterations is less than the preset convergence threshold. The output is a water potential value matrix of a multi-scale hexahedral microenvironment mesh at each time node. Each element in the water potential value matrix corresponds to the calculated water potential value of a hexahedral microenvironment mesh at the corresponding time node.

[0045] The change in water potential gradient is calculated based on the water potential values ​​of adjacent grids. The formula for calculating the water potential gradient is:

[0046] in, For the first The water potential gradient vector of a hexahedral microenvironment grid For the first The volume of each grid, For the first The set of adjacent grids of a grid. Refers to the first The water potential value of each target grid. For adjacent grids water potential value, For the first The grid and the first The distance between the center points of each grid For the first The grid and the first Contact area between grids Let be the unit normal vector of the contact surface, derived from the _th The grid points to the first Each grid.

[0047] The magnitude of the water potential gradient vector is the magnitude of the water potential gradient change in the corresponding grid, and the direction is the direction of the maximum water potential change. By calculating the water potential gradient vectors of all grids, the distribution of water potential gradient changes in the entire underground space under different aboveground vegetation spatial layout schemes can be obtained.

[0048] A graph neural network model was constructed with the optimization objectives of minimizing the range of water potential gradients and minimizing the area of ​​nutrient-deficient regions within a multi-scale hexahedral microenvironment grid. The multi-scale hexahedral microenvironment grid was used as the graph nodes, with the number of nodes matching the number of grid cells. Each graph node corresponded to a unique spatial location within the subsurface grid. The water potential matrix and soil nutrient concentration were used as the initial feature vectors for the graph nodes, with four dimensions: water potential, nitrogen concentration, phosphorus concentration, and potassium concentration. A spatial mapping topological adjacency matrix was used as the edge connections and initial edge weights for the graph neural network model. Non-zero elements in the adjacency matrix corresponded to edge connections between graph nodes, and the element values ​​represented the initial edge weights.

[0049] A graph neural network model with three graph convolutional layers is constructed. The node state is updated by weighted aggregation of feature vectors from adjacent graph nodes based on initial edge weights. The feature aggregation calculation formula for the graph convolutional layer is as follows:

[0050] in, For the first The graph node at the th ... The feature vector output by the layered graph convolutional layer It is a non-linear activation function. For the first The set of neighboring nodes of a graph node For the first The node and the first Edge weights between nodes , The first The node and the first The degree of each node, For the first The trainable weight matrix of a layered graph convolutional layer For the first The node at the th The feature vector output by the layered graph convolutional layer For the first Bias vector of a layered graph convolutional layer.

[0051] The three-layer graph convolutional network sequentially aggregates and updates the feature vectors of graph nodes. The first graph convolutional layer takes the initial feature vector of the graph node as input and outputs the first-layer node features. The second graph convolutional layer takes the first-layer node features as input and outputs the second-layer node features. The third graph convolutional layer takes the second-layer node features as input and outputs the final graph node state vector. The graph node state vector contains the predicted water potential value and predicted nutrient concentration value of the corresponding grid.

[0052] A joint loss function is constructed, using the difference between the maximum and minimum water potential values ​​in the graph node state vector as the water potential gradient range term, and the number of graph nodes in the state vector with nutrient concentrations below the minimum threshold for plant growth as the nutrient deficit area term. These two terms are then weighted and summed to generate the joint loss function. The formula for the joint loss function is:

[0053] in, For the joint loss function value, The maximum water potential value among all graph node state vectors. The minimum water potential value among all graph node state vectors. This represents the number of graph nodes where the nutrient concentration is below the minimum threshold for plant growth. , These are the weighting coefficients for the water potential gradient range term and the nutrient deficit area term, respectively.

[0054] The minimum threshold for plant growth is the minimum concentration of nitrogen, phosphorus, and potassium required for the normal growth of the corresponding plant species. When the concentration of any nutrient element in the grid corresponding to a node is lower than the corresponding minimum threshold, that node is counted as a nutrient-deficient node. The magnitude of the joint loss function corresponds to the degree of deviation between the current plant spatial layout scheme and the optimization objective; the smaller the loss function value, the closer the scheme is to the optimization objective.

[0055] A graph neural network model is used to aggregate node features and update edge weights. Backpropagation gradient descent is then employed to find the optimal spatial configuration of aboveground plant communities. This optimal configuration is mapped onto a digital twin land surface layer to complete landscape planning. The Adam optimizer is used for backpropagation gradient descent, with the goal of minimizing the joint loss function. It iterative updates are performed on the trainable weight matrix and bias vector of the graph neural network model, while simultaneously optimizing the planting location coordinates, species configuration, and planting quantity of aboveground plants.

[0056] During the iteration process, after each weight update, the partial differential equation of the underground three-dimensional transient seepage is recalculated, the initial feature vectors of the graph nodes are updated, graph convolution feature aggregation is performed again, and the value of the joint loss function is calculated until the value of the joint loss function converges to the preset convergence threshold. The iteration then stops, and the current spatial configuration of the aboveground plant community is output as the optimal spatial configuration. The optimal spatial configuration contains complete information on the species number, planting location coordinates, and planting quantity of each plant. The planting location coordinates of the plants in the optimal spatial configuration are mapped to the corresponding three-dimensional canopy model of the aboveground plants onto the surface layer of the digital twin system. The arrangement of the three-dimensional plant models is completed in the surface layer, generating the final landscape planning scheme.

[0057] In this embodiment, a digital computing platform for the underground soil microenvironment was constructed by establishing a tensor field of underground soil porosity and saturation and a multi-scale hexahedral microenvironment grid. Through the spatial mapping topological relationship between the three-dimensional canopy model of aboveground plants and the underground grid, a quantitative correlation between the layout of aboveground plants and the transport of water and fertilizer in underground soil was established. By solving the three-dimensional transient seepage equation, the changes in groundwater potential gradient under different plant layouts were accurately simulated. A graph neural network model was constructed with the optimization objectives of minimizing the maximum difference of water potential gradient and minimizing the area of ​​nutrient deficit region. The solution obtained a spatial combination configuration of plants that conforms to the law of groundwater and soil resource balance. This extends the landscape planning from the geometric arrangement on the ground surface to the quantitative calculation of the underground physical field, avoiding the problem of the disconnect between plant layout and the ecological carrying capacity of underground soil, and ensuring the long-term stability of the landscape ecosystem.

[0058] In a preferred embodiment, an octree structure is used to initially divide the underground space, and the porosity variance within each initially divided grid is calculated. When the porosity variance is greater than a preset geological bedding threshold, the initially divided grid is subdivided into eight equal parts until the porosity variance of all grids is less than the preset geological bedding threshold, thereby generating a multi-scale hexahedral microenvironment grid.

[0059] The preset geological bedding threshold is determined based on the geological survey report of the planning area. The value of the geological bedding threshold is negatively correlated with the homogeneity of the strata; the lower the homogeneity of the strata, the smaller the value of the geological bedding threshold. The initial subdivision level of the octree structure is level 1. The size of the initial grid is consistent with the horizontal dimension of the planning area, and the depth dimension is consistent with the preset maximum root distribution depth. Each level 1 grid is divided into eight level 2 grids, and each level 2 grid can be further divided into eight level 3 grids, and so on, until the minimum size of the grid reaches the preset minimum grid size threshold, at which point the subdivision stops.

[0060] After generating a multi-scale hexahedral microenvironment mesh, the porosity difference between adjacent meshes is extracted, and the distribution gradient of the porosity difference in three-dimensional space is calculated. When the porosity difference between two adjacent meshes exceeds a preset abrupt change threshold, a geological fault or rock stratum boundary is determined to exist between the two meshes, and a transition mesh is inserted between them. The preset abrupt change threshold is three times the preset geological bedding threshold. When the porosity difference between adjacent meshes exceeds the abrupt change threshold, it indicates a sudden change in soil medium properties between the two meshes, which can cause numerical oscillations during water potential calculations. Therefore, a transition mesh needs to be inserted for smoothing.

[0061] The number of transition grids is determined by the ratio of the porosity difference between two adjacent grids to a preset abrupt change threshold. A larger ratio results in more transition grids being inserted. The porosity values ​​of the transition grids are distributed between the porosity values ​​of two adjacent grids using linear interpolation. Laplacian smoothing is applied to the node coordinates of the transition grids to eliminate grid distortion. The Laplacian smoothing process involves adjusting the coordinates of each node in the transition grid to the average of the coordinates of its adjacent nodes, iteratively adjusting until all interior angles of the grid are within a preset reasonable angle range, thus avoiding calculation errors caused by grid distortion.

[0062] The hydraulic parameters of the transition grid are set as interpolated results of the parameters of two adjacent grids. These hydraulic parameters include the porosity tensor, saturation tensor, and hydraulic conductivity tensor. Three-dimensional linear interpolation is used to ensure that the parameters of the transition grid change continuously in space. In subsequent calculations of water potential gradient changes, the transition grid is used as an independent computational node in the Darcy's law solution, employing the same discretization and iterative solution methods as other hexahedral microenvironment grids.

[0063] Table 2. Correspondence between multi-scale hexahedral microenvironment mesh subdivision levels and parameter thresholds.

[0064] As shown in Table 2, the mesh subdivision levels correspond one-to-one with the subdivision levels of the octree structure. The porosity variance threshold is the threshold for determining when to stop subdividing the mesh at the corresponding level. The porosity mutation threshold is the threshold for determining when to insert a transition mesh into an adjacent mesh at the corresponding level. The maximum number of subdivision iterations is the maximum number of subdivisions allowed for the corresponding level. This provides complete parameter constraints for the generation and optimization of multi-scale hexahedral microenvironment meshes, ensuring a balance between mesh generation accuracy and computational efficiency.

[0065] A three-dimensional tomographic scan of the root system of candidate plant communities was performed to extract the spatial orientation coordinate set and root radius dimensions at each level. These coordinate sets and radius dimensions were then substituted into a modified van der Hoff soil hydrodynamic equation, and different soil water potential boundary conditions were input to solve the equation, generating root water uptake rate curves with soil water potential as the independent variable and water flux as the dependent variable. After generating the root water uptake rate curves, historical meteorological data and real-time microclimate monitoring data for the planning area were acquired, and saturated vapor pressure difference and solar radiation flux were extracted from these data.

[0066] Historical meteorological data comprises hourly meteorological observations of the planning area over the past 10 years, including temperature, relative humidity, solar radiation, wind speed, and precipitation. Real-time microclimate monitoring data consists of minute-by-minute monitoring data collected by microclimate sensors deployed within the planning area, including temperature, relative humidity, solar radiation, and wind speed above the canopy. Saturated vapor pressure difference is calculated based on temperature and relative humidity, representing the difference between saturated vapor pressure and actual vapor pressure. Solar radiation flux is the solar radiation energy received per unit area per unit time.

[0067] A stomatal conductance model for the canopy was established. The saturated vapor pressure difference and solar radiation flux were input into the model to calculate the canopy transpiration rate. The Jarvis model was adopted, which expresses stomatal conductance as a product function of solar radiation, saturated vapor pressure difference, air temperature, and soil water potential. Specifically, the calculation process involves first calculating the stomatal conductance constraint function for each environmental factor, multiplying all constraint functions together, and then multiplying by the maximum stomatal conductance to obtain the canopy stomatal conductance. Based on the canopy stomatal conductance, the Penman-Monteith equation was used to calculate the canopy transpiration rate. The Penman-Monteith equation comprehensively considers the influence of radiation and aerodynamic terms on the transpiration rate, accurately simulating the transpiration water consumption process of the plant canopy.

[0068] Please refer to the attached document. Figure 5The canopy transpiration rate is superimposed as a boundary constraint on the water flux, the dependent variable, in the root water uptake rate curve, generating a canopy-root coupled water uptake rate curve. Water consumption from canopy transpiration needs to be replenished through root uptake. Under steady-state conditions, the total water flux from root uptake equals the total water consumption from canopy transpiration. Therefore, the time-series data of canopy transpiration rate is used as a boundary constraint to correct the water flux in the root water uptake rate curve. The corrected water flux is the original water flux from the root water uptake rate curve multiplied by the correction coefficient corresponding to the canopy transpiration rate. The correction coefficient is the ratio of the actual canopy transpiration rate to the potential transpiration rate. When calculating changes in water potential gradient, the canopy-root coupled water uptake rate curve is used instead of the root water uptake rate curve as the source and sink term. This ensures that the source and sink terms in the groundwater seepage calculation reflect the impact of atmospheric environmental changes on plant root water uptake, improving the accuracy of seepage calculations.

[0069] After constructing a spatial mapping topological adjacency matrix describing the hydraulic connectivity between aboveground plants and underground grids, data on the types and concentrations of secondary metabolites secreted by candidate plant communities are obtained to construct a plant allelopathic feature vector. Specifically, secondary metabolites include substances such as phenolic acids, terpenes, and alkaloids that promote or inhibit the growth of surrounding plants. Each dimension of the allelopathic feature vector corresponds to the relative concentration of a type of secondary metabolite, and the dimension of the vector is consistent with the number of detected secondary metabolite types.

[0070] The cosine similarity of allelopathic feature vectors between different plant species is calculated. When the cosine similarity is lower than a preset antagonism threshold, it is determined that the two corresponding plant species have root allelopathic antagonism. The cosine similarity ranges from -1 to 1. The lower the value, the greater the difference in secondary metabolite characteristics between the two plants, and the higher the probability of antagonism. The preset antagonism threshold is determined based on measured data of plant allelopathic effects. When the cosine similarity is lower than the antagonism threshold, it indicates that there is mutual inhibition between the root exudates of the two plants, and they are not suitable for adjacent planting.

[0071] In the spatially mapped topological adjacency matrix, the edge weights between the substantially associated grids corresponding to two plants with root allelopathic antagonism are multiplied by a preset penalty factor for attenuation, generating an allelopathic constraint topological adjacency matrix. The preset penalty factor ranges from (0,1), and its value is negatively correlated with the strength of the antagonism; the stronger the antagonism, the smaller the penalty factor. When two plants have allelopathic antagonism, the edge weights between the substantially associated grids corresponding to the two plants are multiplied by the penalty factor, reducing the feature aggregation strength between the two grids. This avoids planting the two plants adjacently during the optimization process, reducing the impact of allelopathic antagonism on plant growth. In the graph neural network model, the allelopathic constraint topological adjacency matrix replaces the spatially mapped topological adjacency matrix for feature aggregation, allowing the graph neural network's feature aggregation process to consider the allelopathic effects between plants. The optimized plant layout scheme avoids the adjacent planting of antagonistic plants.

[0072] In this embodiment, a multi-scale hexahedral microenvironment mesh is generated through octree subdivision constrained by porosity variance. Transitional meshes are inserted at locations of geological abrupt changes and smoothed using Laplace's algorithm, reducing the numerical error in water potential calculations caused by abrupt changes in formation properties and ensuring the stability and accuracy of seepage calculations. The canopy transpiration rate is calculated using a canopy stomatal conductance model, generating a canopy-root coupled water uptake rate curve. This ensures that the source and sink terms of root water uptake closely reflect the actual atmosphere-plant-soil water exchange process, improving the accuracy of groundwater and nutrient transport simulation. Antagonistic effects are determined using the cosine similarity of plant allelopathy feature vectors, and the edge weights of the topological adjacency matrix are attenuated. This process eliminates adjacent layouts of antagonistic plants during optimization, further improving the ecological adaptability of the planning scheme.

[0073] In another preferred embodiment, a multi-scale hexahedral microenvironment grid is used as graph nodes, the water potential matrix and soil nutrient concentration are used as initial feature vectors of the graph nodes, and the spatial mapping topological adjacency matrix is ​​used as the edge connection relationship and initial edge weights of the graph neural network model. A graph neural network model containing three layers of graph convolutional networks is constructed. The feature vectors of adjacent graph nodes are weighted and aggregated according to the initial edge weights to update the node states, and a joint loss function is constructed. After constructing the joint loss function, in the process of solving the optimal spatial combination configuration of the aboveground plant community using backpropagation gradient descent, the predicted water potential value and predicted nutrient concentration of each graph node output by the graph neural network model are extracted. Specifically, the graph node state vector output by the third layer of the graph neural network model contains the predicted water potential value, predicted nitrogen concentration, predicted phosphorus concentration, and predicted potassium concentration corresponding to each graph node, which correspond to the predicted water potential and nutrient concentration of the underground grid, respectively.

[0074] Based on the law of conservation of mass, the residual between the sum of fluxes flowing into each graph node and the sum of fluxes flowing out of each graph node is calculated. The law of conservation of mass requires that, under steady-state conditions, the difference between the total flux flowing into and out of the control volume equals the sum of the source and sink terms within the control volume. For soil water transport, the difference between the sum of water fluxes flowing into and out of each grid equals the total water consumed by the roots within that grid. For nutrient transport, the difference between the sum of nutrient fluxes flowing into and out of each grid equals the total nutrients absorbed by the plant roots within that grid. The formula for calculating the mass conservation residual is:

[0075] in, For the first The mass conservation residuals corresponding to each graph node For the inflow of the first Total flux of each grid For the outflow of the first Total flux of each grid For the first The sum of source and sink terms within each grid. For water transport, the source and sink terms are the total water flux absorbed by the roots. For nutrient transport, the source and sink terms are the total nutrient flux absorbed by the roots.

[0076] The total inflow and outflow fluxes are calculated based on the predicted water potential, predicted nutrient concentration, and hydraulic conductivity tensors of adjacent grids. The sum of source and sink terms is calculated based on the root water absorption rate curves, nutrient uptake spectra, and weights of the spatially mapped topological adjacency matrix of the plants associated with the corresponding grids. The magnitude of the residuals reflects the degree of conformity between the prediction results and physical conservation laws; the larger the residuals, the greater the deviation of the prediction results from physical laws.

[0077] When the residual exceeds the preset physical tolerance, the physical constraint penalty term is activated, and the square of the residual multiplied by the preset penalty coefficient is added to the joint loss function. The preset physical tolerance is determined based on the required computational accuracy; the smaller the physical tolerance, the higher the requirement for physical consistency of the prediction results. The preset penalty coefficient is a constant greater than 0; the larger the penalty coefficient, the stronger the physical constraint. The formula for calculating the physical constraint penalty term is:

[0078] in, This represents the value of the physical constraint penalty term. To preset the penalty coefficient, This represents the total number of nodes in the graph. For the first The quality-conserving residuals of each graph node.

[0079] Please refer to the attached document. Figure 6 When the residual of any graph node exceeds the preset physical tolerance, a physical constraint penalty term is activated. The value of the physical constraint penalty term is added to the joint loss function, and the updated joint loss function is the sum of the original joint loss function and the physical constraint penalty term. During backpropagation gradient calculation, gradient propagation paths other than the physical constraint penalty term are truncated, forcing backpropagation gradient descent to update the network weight parameters of the graph neural network model along the direction that reduces the residual, until the residual converges to zero.

[0080] The gradient truncation process involves retaining only the gradient of the physical constraint penalty term with respect to the network weight parameters during backpropagation, while setting the gradient of the original joint loss function with respect to the network weight parameters to 0. This ensures that the update of the network weight parameters is determined solely by the gradient of the physical constraint penalty term, forcing the network's prediction results to satisfy the law of conservation of mass and preventing abnormal results that deviate from the true physical laws during network optimization.

[0081] The training process of the graph neural network model is divided into two stages. The first stage is the pre-training stage without physical constraints. The original joint loss function is used as the optimization objective to pre-train the network weight parameters until it converges to a preset pre-training threshold, resulting in the pre-trained network weight parameters. The second stage is the fine-tuning stage with physical constraints. A physical constraint penalty term is activated, and the joint loss function with the added penalty term is used as the optimization objective. Gradient truncation is employed to fine-tune the pre-trained network weight parameters until the mass conservation residuals of all graph nodes are less than a preset physical tolerance, and the original joint loss function converges to a preset convergence threshold. Training then stops, and the final network weight parameters and corresponding plant space combination configuration are output.

[0082] Table 3 shows the layer structure and parameter configuration of the neural network model.

[0083] As shown in Table 3, the input feature dimensions of the input layer correspond to the four dimensions of the initial feature vector of the graph node, namely water potential value, nitrogen concentration, phosphorus concentration and potassium concentration. The three graph convolutional layers sequentially transform and aggregate the feature vectors. The output feature dimensions of the output layer are consistent with those of the input layer, corresponding to the predicted water potential value and predicted nutrient concentration of the graph node. The activation function, weight initialization method and regularization coefficient are the configuration parameters of the corresponding network layer, providing complete parameter constraints for the construction and training of the graph neural network model, and ensuring the convergence and generalization ability of the model.

[0084] In the process of solving the graph neural network model, the spatial combination configuration of the aboveground plant community is constrained. The constraints include the upper limit of plant planting density, the minimum planting spacing, and the proportion requirement of native plants. The constraints are superimposed on the joint loss function in the form of penalty terms. When the plant layout scheme violates the constraints, the corresponding penalty term is activated, increasing the value of the joint loss function. In the optimization process, layout schemes that violate the constraints are avoided.

[0085] The optimal spatial configuration obtained from the solution is mapped onto a digital twin surface layer to complete the landscape planning. The surface layer of the digital twin system contains complete data on the surface elevation model, distribution of surface structures, road system, and water system of the planning area. During the mapping process, the areas of surface structures, roads, and water systems are avoided. The three-dimensional canopy model of the plants is arranged in the available green area of ​​the surface layer according to the planting position coordinates of the optimal spatial configuration. At the same time, the planting elevation of the plants is adjusted to match the surface elevation model, ensuring that the bottom of the plant model is in close contact with the ground. The final planning scheme is visualized in three dimensions in the digital twin system, and a plant list of the planning scheme is output. The list contains complete information on each plant species number, plant name, planting quantity, planting position coordinates, and maintenance requirements, providing data support for actual garden construction and maintenance.

[0086] In this embodiment, by introducing a physical constraint penalty term based on the law of mass conservation into the joint loss function, and truncating the gradient propagation path of non-physical constraints during backpropagation, the prediction results and update process of the graph neural network model are forced to follow the law of physical conservation. This avoids abnormal layout schemes that deviate from the true physical laws during the neural network optimization process, ensuring the physical feasibility of the planning scheme. Through a phased training method, both the optimization efficiency and physical consistency of the model are taken into account. The optimal spatial combination configuration obtained by the solution satisfies both the optimization objectives of the water potential gradient range and nutrient-deficient areas, and conforms to the physical laws of groundwater and soil transport, further improving the stability and reliability of the planning scheme in practical applications.

Claims

1. A smart landscape planning method based on digital twins, characterized in that, include: Obtain three-dimensional geological profile data of the planning area, establish a tensor field of underground soil porosity and saturation based on the three-dimensional geological profile data, and divide the underground space into a multi-scale hexahedral microenvironment grid in the digital twin system. Root morphology parameters of candidate plant communities were obtained, and a plant physiological feature library containing root water absorption rate curves and nutrient uptake spectra was constructed. A spatial mapping topological relationship between the three-dimensional canopy model of aboveground plants and the multi-scale hexahedral microenvironment mesh is established in the digital twin system; Based on the spatial mapping topology and the tensor fields of underground soil porosity and saturation, the water potential gradient changes of the multi-scale hexahedral microenvironment grid under different aboveground plant spatial layout schemes are calculated. A graph neural network model is constructed with the optimization objectives of minimizing the water potential gradient range and minimizing the area of ​​nutrient-deficient regions in the multi-scale hexahedral microenvironment grid. The graph neural network model is used to aggregate node features and update edge weights. The optimal spatial combination configuration of the aboveground plant community is solved by backpropagation gradient descent. The optimal spatial combination configuration is then mapped to a digital twin surface layer to complete the landscape planning.

2. The intelligent landscape planning method based on digital twins according to claim 1, characterized in that, The process involves acquiring three-dimensional geological profile data of the planning area, establishing a tensor field for underground soil porosity and saturation based on the three-dimensional geological profile data, and dividing the underground space into a multi-scale hexahedral micro-environment grid in the digital twin system. This includes: acquiring electromagnetic wave reflection signals of the planning area through ground-penetrating radar, performing time-frequency domain transformation on the electromagnetic wave reflection signals to extract the dielectric constants of strata at different depths, and calculating the initial soil porosity distribution based on the dielectric constants. Multi-point saturation data collected by a soil moisture sensor network are acquired, and the underground soil porosity and saturation tensor field is constructed by combining the initial soil porosity distribution. The underground space is initially divided using an octree structure. The porosity variance within each initially divided grid is calculated. When the porosity variance is greater than a preset geological bedding threshold, the initially divided grid is subdivided into eight equal parts until the porosity variance of all grids is less than the preset geological bedding threshold, thereby generating the multi-scale hexahedral microenvironment grid.

3. The intelligent landscape planning method based on digital twins according to claim 1, characterized in that, Obtain root morphological parameters of candidate plant communities and construct a plant physiological feature library including root water absorption rate curves and nutrient absorption spectra. This includes: performing three-dimensional tomographic scanning on the root system of the candidate plant communities, extracting the root spatial orientation coordinate set and the radius dimensions of each root level, substituting the coordinate set and the radius dimensions into the modified van der Hoff soil hydrodynamic equation, inputting different soil water potential boundary conditions to solve, and generating the root water absorption rate curve with soil water potential as the independent variable and water flux as the dependent variable. Ion adsorption kinetic parameters of nitrogen, phosphorus and potassium elements in the candidate plant communities are extracted, and the nutrient absorption spectrum is generated by fitting it along the time dimension. The root water absorption rate curve and the nutrient absorption spectrum are associated with the plant species number and stored in the plant physiological feature database.

4. The intelligent landscape planning method based on digital twins according to claim 1, characterized in that, Establishing a spatial mapping topological relationship between the three-dimensional canopy model of aboveground plants and the multi-scale hexahedral microenvironment mesh in the digital twin system includes: obtaining the vertical projection polygon of the three-dimensional canopy model of aboveground plants on the ground surface plane, and determining the multi-scale hexahedral microenvironment mesh within the coverage area of ​​the vertical projection polygon as candidate associated meshes; Construct the root envelope equation of the candidate plant community, perform spatial Boolean intersection operation between the root envelope equation and the candidate associated grid, and extract the set of intersecting grids as the actual associated grids; Calculate the Euclidean distance from the center point of the three-dimensional canopy model of the aboveground plants to the center point of each grid in the substantially associated grid, and construct a spatial mapping topological adjacency matrix describing the hydraulic connectivity strength between the aboveground plants and the underground grid using the reciprocal of the Euclidean distance as the weight.

5. The intelligent landscape planning method based on digital twins according to claim 1, characterized in that, Based on the spatial mapping topology and the tensor field of underground soil porosity and saturation, the water potential gradient change of the multi-scale hexahedral microenvironment grid under different aboveground plant spatial layout schemes is calculated, including: for any aboveground plant spatial layout scheme, the root water absorption rate curve of the corresponding plant is retrieved from the plant physiological feature database as the source and sink term; Using the underground soil porosity and saturation tensor fields as initial boundary conditions, and combining them with the spatial mapping topological adjacency matrix, Darcy's law and Richards' equations are substituted to construct the underground three-dimensional transient seepage partial differential equation. The partial differential equation of the underground three-dimensional transient seepage is discretized using the finite volume method, and iteratively solved within a set time step. The water potential value matrix of the multi-scale hexahedral microenvironment grid at each time node is output, and the water potential gradient change is calculated based on the water potential values ​​of the adjacent grids.

6. The intelligent landscape planning method based on digital twins according to claim 1, characterized in that, With the optimization objectives of minimizing the water potential gradient range and minimizing the area of ​​nutrient deficiency region in the multi-scale hexahedral microenvironment grid, a graph neural network model is constructed, including: using the multi-scale hexahedral microenvironment grid as graph nodes, and using the water potential matrix and soil nutrient concentration as the initial feature vectors of the graph nodes; The spatial mapping topological adjacency matrix is ​​used as the edge connection relationship and initial edge weight of the graph neural network model; Construct a graph neural network model containing three layers of graph convolutional network, and update the node state by weighted aggregation of the feature vectors of adjacent graph nodes based on the initial edge weights; A joint loss function is constructed, in which the difference between the maximum and minimum water potential values ​​in the graph node state vector is taken as the water potential gradient range term, and the number of graph nodes in the state vector with nutrient concentrations below the minimum threshold for plant growth is taken as the nutrient deficit area term. The two terms are weighted and summed to generate the joint loss function.

7. The intelligent landscape planning method based on digital twins according to claim 2, characterized in that, After generating the multi-scale hexahedral microenvironment mesh, the method further includes: extracting the porosity difference between adjacent meshes in the multi-scale hexahedral microenvironment mesh, and calculating the distribution gradient of the porosity difference in three-dimensional space; When the porosity difference between two adjacent grids is greater than a preset abrupt change threshold, it is determined that there is a geological fault or rock stratum boundary between the two grids, and a transition grid is inserted between the two grids. The node coordinates of the transition mesh are smoothed using Laplacian to eliminate mesh distortion, and the hydraulic parameters of the transition mesh are set to the interpolation calculation results of the parameters of two adjacent meshes. In subsequent calculations of the water potential gradient change, the transition grid is used as an independent calculation node to participate in the Darcy's law solution.

8. The intelligent landscape planning method based on digital twins according to claim 3, characterized in that, After generating the root water absorption rate curve, the method further includes: acquiring historical meteorological data and real-time microclimate monitoring data of the planning area, and extracting saturated water vapor pressure difference and solar radiation flux from the historical meteorological data and the real-time microclimate monitoring data. A stomatal conductance model of the canopy is established, and the saturated water vapor pressure difference and the solar radiation flux are input into the stomatal conductance model of the canopy to calculate the canopy transpiration rate; The canopy transpiration rate is superimposed as a boundary constraint on the water flux, the dependent variable of the root water uptake rate curve, to generate a canopy-root coupled water uptake rate curve. When calculating the water potential gradient change, the canopy-root coupled water absorption rate curve is used to replace the root water absorption rate curve as the source and sink term.

9. The intelligent landscape planning method based on digital twins according to claim 4, characterized in that, After constructing a spatial mapping topological adjacency matrix describing the hydraulic connectivity strength between aboveground plants and underground grids, the method further includes: obtaining data on the types and concentrations of secondary metabolites secreted by the candidate plant communities, and constructing a plant allelopathy feature vector. Calculate the cosine similarity of the plant allelopathic feature vectors between different plant species. When the cosine similarity is lower than a preset antagonism threshold, it is determined that the two corresponding plant species have root allelopathic antagonism. In the spatial mapping topological adjacency matrix, the edge weights between the substantially associated grids corresponding to the two plants with root allelopathic antagonism are multiplied by a preset penalty factor for attenuation processing to generate an allelopathic constraint topological adjacency matrix. In the graph neural network model, the chemical constraint topological adjacency matrix is ​​used to replace the spatial mapping topological adjacency matrix for feature aggregation.

10. The intelligent landscape planning method based on digital twins according to claim 6, characterized in that, After constructing the joint loss function, the method further includes: extracting the predicted water potential and predicted nutrient concentration of each graph node output by the graph neural network model during the process of solving the optimal spatial combination configuration of the aboveground plant community using backpropagation gradient descent. Based on the law of conservation of mass, calculate the residual between the sum of the flux flowing into each graph node and the sum of the flux flowing out of the graph node; When the residual is greater than the preset physical tolerance, the physical constraint penalty term is activated, and the square of the residual is multiplied by the preset penalty coefficient and added to the joint loss function; When calculating the gradient during backpropagation, the gradient propagation path other than the physical constraint penalty term is truncated, and the backpropagation gradient descent is forced to update the network weight parameters of the graph neural network model along the direction of reducing the residual, until the residual converges to 0.

Images