A method for optimizing land space ecological restoration scheme based on GIS and multi-objective decision analysis

Abstract

The invention provides a land space ecological restoration scheme optimization method based on GIS and multi-objective decision analysis, which belongs to the field of ecological restoration and solves the problem of low ecological restoration efficiency. The method is specifically as follows: step S1: geographic space model and ecological environment data; step S2: determining the vegetation change rate, water flow change rate and sand content change rate in the area to be optimized, and determining the plant restoration scheme in combination with the geographic space model and ecological environment data; step S3: dividing the protection range for organisms based on the geographic space model and ecological environment data, and determining the animal restoration scheme; step S4: summarizing the plant restoration scheme and the animal restoration scheme and feeding them back to a user; the invention obtains, analyzes and processes the geographic information and ecological information of the area to be optimized, determines the planting area and the biological protection area for the area to be optimized, and improves the ecological restoration efficiency.

Description

Technical Field

[0001] The present invention discloses a land space ecological restoration program optimization method based on GIS and multi-objective decision analysis, and relates to the field of ecological restoration. Background Art

[0002] The existing methods of ecological restoration have the following shortcomings:

[0003] Strong data dependence: The realization of GIS functions is highly dependent on the accuracy and completeness of data. Data quality issues may lead to deviations in analysis results.

[0004] Regional differences: Ecological problems in different regions are different and need to be adjusted appropriately according to specific circumstances.

[0005] High time cost: Ecological restoration is a long process that requires continuous monitoring and evaluation, and has a relatively high time cost.

[0006] There is subjectivity: In multi-objective decision analysis, the weight allocation of each objective may be affected by the subjective judgment of the decision maker. Summary of the invention

[0007] In view of the shortcomings of the existing technology, the purpose of the present invention is to provide a land space ecological restoration plan optimization method based on GIS and multi-objective decision analysis, aiming to solve the problem of slow ecological restoration rate.

[0008] In order to achieve the above object, the present invention is implemented through the following technical solutions: A method for optimizing land space ecological restoration scheme based on GIS and multi-objective decision analysis includes:

[0009] Step S1: Using the GIS system, obtain the terrain information of the area to be optimized and obtain a geographic space model;

[0010] Obtain soil quality information, climate information, river pollution information and biodiversity information of the area to be optimized, and obtain ecological environment data;

[0011] Step S2: Obtain the distribution area information and coverage area information of vegetation and rivers in the area to be optimized, and determine the vegetation change rate, water flow change rate and sediment content change rate;

[0012] Using the GIS model analysis component, combined with ecological environment data, based on the geographic spatial model, simulate and analyze the soil and water loss in the area to be optimized, and determine the plant restoration plan;

[0013] Step S3: Based on the geospatial model, the diffusion trend of vegetation and the accumulation trend of biological harmful substances are simulated according to the soil quality information, river pollution information and biodiversity information;

[0014] Based on the diffusion and accumulation trends, the protection range is divided for organisms and animal restoration plans are determined;

[0015] Step S4: Summarize the plant restoration plan and the animal restoration plan and feed back to the user;

[0016] Furthermore, the specific steps of step S1 are as follows:

[0017] Step S11: Obtain elevation points and contour lines of the area to be optimized to obtain vector data;

[0018] Step S12: load the scipy library and use the Delaunay triangulation function to convert the vector data into TIN data;

[0019] Step S13: Load the ArcGIS software package and use the TIN to Raster tool to convert the TIN data into a DEM to obtain a geospatial model;

[0020] Step S14: Acquire soil quality information, climate information, river pollution information and biodiversity information in the area to be optimized to obtain ecological environment data;

[0021] Soil quality information includes: average microbial content, inorganic content and soil moisture content of the soil in the area to be optimized;

[0022] River pollution information includes: pollutant content corresponding to each river in the area to be optimized;

[0023] Climate information includes: the start time and duration of the dry season in the area to be optimized, as well as the average, maximum and minimum values ​​of rainfall;

[0024] The starting time and duration of the rainy season, as well as the average, maximum and minimum values ​​of rainfall;

[0025] Biodiversity information includes: the population size and distribution range of all protected plants and animals in the area to be optimized;

[0026] Step S15: Summarize the geographic space model and ecological environment data and proceed to step S2.

[0027] Furthermore, the specific steps of step S2 are as follows:

[0028] Step S21: using the same steps of obtaining the geographic space model, obtain the distribution area and coverage area of ​​the vegetation in the area to be optimized; the vegetation coverage area is recorded as s;

[0029] Obtain the river basin, water flow and sediment content in the area to be optimized; record the water flow as f; record the sediment content as sc;

[0030] Step S22: according to the distribution area of ​​vegetation, the historical vegetation coverage area of ​​the area to be optimized in the past five years is obtained, which is recorded as hs1, hs2, hs3, hs4 and hs5;

[0031] According to the river basin, obtain the historical water flow and sediment content of the area to be optimized in the past five years;

[0032] Historical water flows, denoted as hf1, hf2, hf3, hf4 and hf5;

[0033] Historical sediment content, recorded as hsc1, hsc2, hsc3, hsc4 and hsc5;

[0034] Step S23: Define calculation formula 21 to calculate the vegetation coverage area growth rate corresponding to the area to be optimized: us 1～2 、us 2～3 、us 3～4 and us 4～5 ;

[0035] Water flow rate growth rate: uf 1～2 uf 2～3 uf 3～4 and uf 4～5 ;

[0036] Sediment content growth rate: usc 1～2 、usc 2～3 、usc 3～4 and usc 4～5 ;

[0037] Calculation 21: Where X represents the historical vegetation coverage area (hs1~hs5), historical water flow (hf1~hf5) and historical sediment content (hsc1~hsc5); u(X) represents the growth rate corresponding to X;

[0038] Step S24: determining the vegetation change rate qs;

[0039] Step S25: Repeat the same steps of determining the vegetation change rate qs to determine the change rates of water flow and sediment content, which are recorded as qf and qsc respectively;

[0040] Step S26: Simulate and analyze the soil erosion in the area to be optimized and determine a phytoremediation plan.

[0041] Furthermore, the specific steps of step S24 are as follows:

[0042] Step S241: extract and calculate us 1～2 ~us 4～5The maximum, minimum and average values ​​are recorded as mus, lus and aus respectively; compare (mus+aus) with 2us 4～5 size;

[0043] Step S242: If (mus+aus)<2us 4～5 , it indicates that the vegetation coverage growth rate of the area to be optimized is increasing, and the subsequent steps S243 to S244 are skipped. The value of qs is us 4～5 ;

[0044] If (mus+aus)>2us 4～5 , it means that the vegetation coverage area growth rate of the area to be optimized is on a downward trend, and the subsequent steps S243 to S244 are skipped, and the value of qs is lus;

[0045] If (mus+aus)=2us 4～5 , it means that the vegetation coverage area growth rate of the area to be optimized is unstable, and the process goes to step S243;

[0046] Step S243: Calculate us 1～2 ~us 4～5 In the calculation, the sum of positive growth rates is recorded as zus; the sum of negative growth rates is recorded as fus; comparison (zus) 2 and (fus) 2 The size of qs is determined, and the process goes to step S244;

[0047] Step S244: If (zus) 2 ≥(fus) 2 , then the calculation formula of qs is as follows:

[0048]

[0049] If (zus) 2 ≥(fus) 2 , then the calculation formula of qs is as follows:

[0050]

[0051] Furthermore, the specific steps of step S26 are as follows:

[0052] Step S261: Set the vegetation area as color 1 and the river area as color 2; read the RGB values ​​corresponding to color 1 and color 2, and record them as R1 and R2 respectively;

[0053] Step S262: Determine the influence coefficient aa1 of the dry season on the water flow;

[0054] Step S2621: the average value, the maximum value and the minimum value corresponding to the dry season rainfall are recorded as ad, md and ld respectively;

[0055] Step S2622: define relation 22: lim(ad)→md;

[0056] Define relation 23: lim(ad)→ld;

[0057] Step S2623: If definition formula 22 is satisfied, then the influence coefficient aa1, aa1 = (ad + md) / 2;

[0058] If definition 23 is satisfied, then the influence coefficient aa1, aa1 = (ad + ld) / 2;

[0059] If neither definition 22 nor 23 is satisfied, then the influence coefficient is aa1, aa1 = (md / ad+ad / ld) / 2;

[0060] Step S263: Repeat the same steps of determining the influence coefficient aa1 to determine the influence coefficient aa2 of the rainy season on the water flow;

[0061] Step S264: According to the historical vegetation coverage area (hs1-hs5) and the historical sand content (hsc1-hsc5), construct the relationship equation 1: y=a1x1+b1; wherein y represents the sand content; x1 represents the vegetation coverage area; a1 and b1 represent the slope and intercept of the relationship equation 1 respectively;

[0062] Step S265: Repeat the same steps of constructing relational equation 1, and construct relational equation 2 according to historical water flow (hf1~hf5) and historical sediment content (hsc1~hsc5): y=a2x2+b2; wherein x2 represents water flow, a2 and b2 represent the slope and intercept of relational equation 2 respectively.

[0063] Furthermore, the subsequent steps of step S265 are as follows:

[0064] Step S266: taking the geographic space model as the model benchmark, taking s, f and sc as the basic parameters, taking the vegetation change rate qs and the influence coefficient aa1 as the change parameters, and taking the relationship equation 1 and the relationship equation 2 as the constraint conditions;

[0065] Taking the start time and duration of the dry season as simulation period 1, the CSLE generation component in GIS is called to generate the dry season model;

[0066] Taking the start time and duration of the rainy season as simulation period 2, the CSLE generation component in GIS is called to generate the rainy season model;

[0067] Step S267: Determine the vegetation diffusion direction and R2 change relationship in the rainy season model or the dry season model as a supplementary constraint condition;

[0068] Step S268: Substitute the supplementary constraint condition into the constraint condition (step S266); execute step S266 again to construct a new dry season model and a new rainy season model;

[0069] Step S269: Determine a phytoremediation plan;

[0070] Step S2691: using the model overlap component of DIS to overlay the new dry season model and the new rainy season model;

[0071] Step S2692: marking the area where the new dry season model and the new rainy season model overlap and the color is R1 as target area 1;

[0072] The area where the new dry season model and the new rainy season model overlap and the color is ΔR2<R2 is the planting area 2;

[0073] Step S2693: in the target area 1, extracting an area not including the original area as the planting area 1;

[0074] Planting Area 1 and Planting Area 2 are used as phytoremediation solutions.

[0075] Furthermore, the specific steps of step S264 are as follows:

[0076] Step S2641: Import the LinearRegression class from the sklearn.linear_model module, and then import the NumPy library;

[0077] Step S2642: create a receiving array 1 for storing hs1 to hs5; create a receiving array 2 for storing hsc1 to hsc5;

[0078] Step S2643: Create a receiving object model instantiated by the LinearRegression class;

[0079] Step S2644: calling the fit method through model, and substituting the receiving array 1 and the receiving array 2 as parameters into the fit method, obtaining the relationship equation 1: y=a0x1+b0; wherein y represents the sand content; x1 represents the vegetation coverage area; a0 and b0 represent the slope and intercept respectively, and the initial values ​​of a0 and b0 are both 0;

[0080] Step S2645: Call the coef method and the intercept method through the model to calculate a0 and b0 respectively, and obtain a1 and b1; where a1 is a rational number and is not equal to 0; b1 is an arbitrary rational number;

[0081] Step S2646: Summarize the processing results of steps S2641 to S2645 to obtain relationship equation 1: y=a1x1+b1.

[0082] Furthermore, the specific steps of step S267 are as follows:

[0083] Step S2671: Mark the vegetation-covered area in the geographic space model as the original area; read the water flow of the river that intersects with the original area and has the shortest spatial distance as the original flow;

[0084] Step S2672: repeat the step of reading the original area and the original flow in step S2671, and in the rainy season model or the dry season model, read the adjacent vegetation coverage area that does not intersect with the original area as the adjacent area, and read the river water flow corresponding to the adjacent area as the adjacent flow;

[0085] In the rainy season model or dry season model, the direction of vegetation diffusion is from the original area to the adjacent area where the adjacent flow rate is ≥ the original flow rate, and the vegetation area after diffusion is marked by R1;

[0086] Step S2673: Mark the river area in the geographic space model as the target river; read the sediment content of the target river as the original sediment content, recorded as osc;

[0087] Read the sediment content corresponding to the target river in the rainy season model or the dry season model as the new sediment content, recorded as gsc;

[0088] Step S2674: Compare the magnitudes of osc and gsc to determine the R2 variation relationship ΔR2;

[0089] If gsc≥osc, it means that the sediment content increases, ΔR2=[1+log255(gsc-osc)]*R2;

[0090] If gsc<osc, it means that the sand content is reduced, ΔR2=[1-log255(osc-gsc)]*R2.

[0091] Furthermore, the specific steps of step S3 are as follows:

[0092] Step S31: Obtain the slope grid, use the Spatial Analyst tool in GIS to analyze the geospatial model, and divide the terrain of the geospatial model into flat areas, transition areas, and steep areas;

[0093] Step S32: taking the calibration range as a scanning unit, randomly selecting an area covered by vegetation in the geographic space model as a monitoring area; scanning the area around the monitoring area as a neighboring area;

[0094] Step S33: reading and calculating the soil quality information corresponding to the monitoring area, and recording the average microbial content, inorganic matter content and soil moisture content as ws, nw and hh respectively;

[0095] Read and calculate the soil quality information corresponding to the adjacent area, and record the average microbial content, inorganic matter content and soil water content as bws, bnw and bhh respectively;

[0096] Step S34: using the LCG algorithm to randomly generate weight coefficients, recorded as wa1, wa2 and wa3; wa1, wa2 and wa3, satisfy: wa1, wa2 and wa3 are all rational numbers, wa1+wa2+wa3=1, and wa1>wa2>wa3;

[0097] Step S35: defining the constraint conditions of vegetation diffusion trend;

[0098] Constraint 31:

[0099] If wa1*bws+wa2*bnw+wa3*bhh>wa1*ws+wa2*nw+wa3*hh, then the direction of vegetation diffusion is: from the monitoring area to the adjacent area;

[0100] Constraint 32: If wa1*bws+wa2*bnw+wa3*bhh≤wa1*ws+wa2*nw+wa3*hh, then determine whether the terrain of the monitoring area satisfies: steep area → transition area → flat area;

[0101] If satisfied, the direction of vegetation diffusion is: from the monitoring area to the adjacent area;

[0102] If not, compare other neighboring areas of the monitoring area until constraint condition 31 or constraint condition 32 is satisfied;

[0103] Step S36: Based on step S35, supplement the constraint condition of the cumulative trend;

[0104] Sum up the populations of protected animals and plants, recorded as pp and pa;

[0105] Summarize the number of rivers in the area to be optimized, denoted as rn; summarize the pollutant content corresponding to each river, denoted as hn1~hn rn ;

[0106] Constraint 33: The migration direction of the distribution range of protected animals and plants is the same as the diffusion direction of vegetation;

[0107] Constraint 34: The cumulative trend of protected plants increases, and the added value is (phn) 1 / 2 / pp 2; Among them, phn represents the pollutant content of the river corresponding to the “distribution range of protected plants”;

[0108] The cumulative trend of protected animals is increasing, and the added value is (ahn*pp) 1 / 2 / pa; where ahn represents the pollutant content in the river corresponding to the “distribution range of protected plants”;

[0109] Step S37: Call the same steps of analyzing the plant restoration plan, and build a cumulative model of harmful substances to animals and plants based on the geographic space model and constraints 31 to 34; mark the area with the least accumulated harmful substances to animals and plants in the cumulative model of harmful substances to animals and plants as the protection range, and obtain an animal restoration plan.

[0110] Compared with the prior art, the present invention has the following beneficial effects:

[0111] Strong data integration capability: The present invention can integrate geographic data from multiple sources to provide comprehensive data support for ecological restoration programs. It uses GIS technology to display complex spatial data in a graphical form to more intuitively evaluate ecological conditions.

[0112] Multi-objective decision analysis: The present invention can comprehensively consider soil erosion, animal and plant protection, and environmental pollution to achieve the optimization of ecological restoration plans for decision makers to choose from.

[0113] High flexibility: The present invention makes dynamic adjustments according to the ecological conditions in different regions and time periods, predicts future environmental changes based on the existing ecological environment, selects the optimal adjustment area, and improves the efficiency and effectiveness of ecological restoration. BRIEF DESCRIPTION OF THE DRAWINGS

[0114] Other features, objects and advantages of the present invention will become more apparent from the detailed description of non-limiting embodiments made with reference to the following drawings:

[0115] Figure 1 It is a schematic diagram of the method of the present invention;

[0116] Figure 2 This is a schematic diagram of the diffusion of the original area to the adjacent area of ​​the present invention;

[0117] Figure 3 This is a logical diagram of the constraints of the present invention. DETAILED DESCRIPTION

[0118] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0119] See also Figure 1, a land space ecological restoration program optimization method based on GIS and multi-objective decision analysis includes:

[0120] Step S1: Using the GIS system, obtain the terrain information of the area to be optimized and obtain a geographic space model;

[0121] Obtain soil quality information, climate information, river pollution information and biodiversity information of the area to be optimized, and obtain ecological environment data;

[0122] It should be noted that “areas to be optimized” refer to areas “that need ecological restoration”.

[0123] The specific steps of step S1 are as follows:

[0124] Step S11: Associating with a geographic information collection device; obtaining elevation points and contour lines of the area to be optimized through the geographic information collection device to obtain vector data;

[0125] Geographic information collection equipment includes: total stations, GPS receivers, remote sensing satellites, etc., external equipment used for surveying terrain;

[0126] Step S12: Load the scipy library and use the Delaunay triangulation function to convert the vector data into TIN (triangulated irregular network) data;

[0127] Step S13: Load the ArcGIS software package and use the TIN to Raster tool to convert the TIN data into a DEM (digital elevation model) to obtain a geospatial model;

[0128] Step S14: obtaining soil quality information, climate information, river pollution information and biodiversity information in the area to be optimized through the environmental supervision department to obtain ecological environment data;

[0129] Soil quality information includes: average microbial content, inorganic content and soil moisture content corresponding to the soil in the area to be optimized (per 10 cubic meters);

[0130] River pollution information includes: pollutant content corresponding to each river in the area to be optimized;

[0131] Climate information includes: the start time and duration of the dry season in the area to be optimized, as well as the average, maximum and minimum values ​​of rainfall;

[0132] The starting time and duration of the rainy season, as well as the average, maximum and minimum values ​​of rainfall;

[0133] Biodiversity information includes: the population size and distribution range of all protected plants and animals in the area to be optimized;

[0134] Step S15: Summarize the geographic space model and ecological environment data and proceed to step S2.

[0135] Step S2: Obtain the distribution area information and coverage area information of vegetation and rivers in the area to be optimized, and determine the vegetation change rate, water flow change rate and sediment content change rate;

[0136] Using the GIS model analysis component, combined with ecological environment data, based on the geographic space model, simulate and analyze the soil erosion in the area to be optimized, and determine the plant restoration plan;

[0137] The specific steps of step S2 are as follows:

[0138] Step S21: using the same steps of obtaining the geographic space model, obtain the distribution area and coverage area of ​​the vegetation in the area to be optimized; the vegetation coverage area is recorded as s;

[0139] Obtain the river basin, water flow and sediment content (per cubic meter) in the area to be optimized; record the water flow as f; record the sediment content as sc;

[0140] Step S22: according to the distribution area of ​​vegetation, the historical vegetation coverage area of ​​the area to be optimized in the past five years is obtained, which is recorded as hs1, hs2, hs3, hs4 and hs5;

[0141] According to the river basin, obtain the historical water flow and sediment content of the area to be optimized in the past five years;

[0142] Historical water flows, denoted as hf1, hf2, hf3, hf4 and hf5;

[0143] Historical sediment content, recorded as hsc1, hsc2, hsc3, hsc4 and hsc5;

[0144] Step S23: Define calculation formula 21 to calculate the vegetation coverage area growth rate corresponding to the area to be optimized: us 1～2 、us 2～3 、us 3～4 and us 4～5 ;

[0145] Water flow rate growth rate: uf 1～2 uf 2～3 uf 3～4 and uf 4～5 ;

[0146] Sediment content growth rate: usc 1～2 、usc 2～3 、usc 3～4 and usc 4～5 ;

[0147] Calculation 21: Where X represents the historical vegetation coverage area (hs1~hs5), historical water flow (hf1~hf5) and historical sediment content (hsc1~hsc5); u(X) represents the growth rate corresponding to X;

[0148] Step S24: determining the vegetation change rate qs;

[0149] Step S241: extract and calculate us 1～2 ~us 4～5 The maximum, minimum and average values ​​are recorded as mus, lus and aus respectively; compare (mus+aus) with 2us 4～5 size;

[0150] Step S242: If (mus+aus)<2us 4～5 , it indicates that the vegetation coverage growth rate of the area to be optimized is increasing, and the subsequent steps S243 to S244 are skipped. The value of qs is us 4～5 ;

[0151] If (mus+aus)>2us 4～5 , it means that the vegetation coverage area growth rate of the area to be optimized is on a downward trend, and the subsequent steps S243 to S244 are skipped, and the value of qs is lus;

[0152] If (mus+aus)=2us 4～5 , it means that the vegetation coverage area growth rate of the area to be optimized is unstable, and the process goes to step;

[0153] Step S243: Calculate us 1～2 ~us 4～5 In the calculation, the sum of positive growth rates is recorded as zus; the sum of negative growth rates is recorded as fus; comparison (zus) 2 and (fus) 2 The size of qs is determined, and the process goes to step S244;

[0154] It should be noted that the “positive growth rate” in “step S243” refers to us 1～2 ~us 4～5 The vegetation coverage area growth rate is greater than or equal to 0; "negative growth rate" refers to us 1～2 ~us 4～5 The vegetation coverage area growth rate is less than 0.

[0155] Step S244: If (zus) 2 ≥(fus) 2 , then the calculation formula of qs is as follows:

[0156]

[0157] If (zus) 2 ≥(fus) 2 , then the calculation formula of qs is as follows:

[0158]

[0159] Step S25: Repeat the same steps of determining the vegetation change rate qs to determine the change rates of water flow and sediment content, which are recorded as qf and qsc respectively;

[0160] Step S26: Simulate and analyze the soil erosion in the area to be optimized and determine the plant restoration plan;

[0161] Step S261: Set the vegetation area as color 1 and the river area as color 2; read the RGB values ​​corresponding to color 1 and color 2, and record them as R1 and R2 respectively;

[0162] It should be noted that "Color 1" and "Color 2" in the present invention are only used to distinguish vegetation and rivers in the geographic space model; generally, the GGB values ​​of "Color 1" and "Color 2" are 150 and 28 respectively; the specific "Color 1" and "Color 2" can be changed by the user according to actual needs.

[0163] Step S262: Determine the influence coefficient aa1 of the dry season on the water flow;

[0164] Step S2621: the average value, the maximum value and the minimum value corresponding to the dry season rainfall are recorded as ad, md and ld respectively;

[0165] Step S2622: define relation 22: lim(ad)→md;

[0166] Define relation 23: lim(ad)→ld;

[0167] Step S2623: If definition formula 22 is satisfied, then the influence coefficient aa1, aa1 = (ad + md) / 2;

[0168] If definition 23 is satisfied, then the influence coefficient aa1, aa1 = (ad + ld) / 2;

[0169] If neither definition 22 nor 23 is satisfied, then the influence coefficient is aa1, aa1 = (md / ad+ad / ld) / 2;

[0170] Step S263: Repeat the same steps of determining the influence coefficient aa1 to determine the influence coefficient aa2 of the rainy season on the water flow;

[0171] Step S264: constructing relation equation 1 according to the historical vegetation coverage area (hs1-hs5) and the historical sand content (hsc1-hsc5);

[0172] Step S2641: Import the LinearRegression class from the sklearn.linear_model module, and then import the NumPy library;

[0173] Step S2642: create a receiving array 1 for storing hs1 to hs5; create a receiving array 2 for storing hsc1 to hsc5;

[0174] Step S2643: Create a receiving object model instantiated by the LinearRegression class;

[0175] Step S2644: calling the fit method through model, and substituting the receiving array 1 and the receiving array 2 as parameters into the fit method (wherein the receiving array 1 is the independent variable, on the left of the fit method parameters; the receiving array 2 is the dependent variable, on the right of the fit method parameters) to obtain the relationship equation 1: y=a0x1+b0; where y represents the sand content; x1 represents the vegetation coverage area; a0 and b0 represent the slope and intercept respectively, and the initial values ​​of a0 and b0 are both 0;

[0176] Step S2645: Call the coef method and the intercept method through the model to calculate a0 and b0 respectively, and obtain a1 and b1; where a1 is a rational number and is not equal to 0; b1 is an arbitrary rational number;

[0177] Step S2646: Summarize the processing results of steps S2641 to S2645 to obtain relational equation 1: y=a1x1+b1;

[0178] Step S265: Repeat the same steps of constructing relational equation 1, and construct relational equation 2 according to the historical water flow (hf1-hf5) and the historical sediment content (hsc1-hsc5): y=a2x2+b2; wherein x2 represents the water flow;

[0179] Step S266: taking the geographic space model as the model benchmark, taking s, f and sc as the basic parameters, taking the vegetation change rate qs and the influence coefficient aa1 as the change parameters, and taking the relationship equation 1 and the relationship equation 2 as the constraint conditions;

[0180] Taking the start time and duration of the dry season as simulation period 1, the CSLE generation component in GIS is called to generate the dry season model;

[0181] Taking the start time and duration of the rainy season as simulation period 2, the CSLE generation component in GIS is called to generate the rainy season model;

[0182] Step S267: Determine the vegetation diffusion direction and R2 change relationship in the rainy season model or the dry season model as a supplementary constraint condition;

[0183] Step S2671: Mark the vegetation-covered area in the geographic space model as the original area; read the water flow of the river that intersects with the original area and has the shortest spatial distance as the original flow;

[0184] Step S2672: Please refer to Figure 2 , repeat the step S2671 of reading the original area and the original flow, in the rainy season model or the dry season model, read the adjacent vegetation coverage area that does not intersect with the original area as the adjacent area, and read the river water flow corresponding to the adjacent area as the adjacent flow;

[0185] In the rainy season model or dry season model, the direction of vegetation diffusion is from the original area to the adjacent area where the adjacent flow rate is ≥ the original flow rate, and the vegetation area after diffusion is marked by R1;

[0186] Step S2673: Mark the river area in the geographic space model as the target river; read the sediment content of the target river as the original sediment content, recorded as osc;

[0187] Read the sediment content corresponding to the target river in the rainy season model or the dry season model as the new sediment content, recorded as gsc;

[0188] Step S2674: Compare the magnitudes of osc and gsc to determine the R2 variation relationship ΔR2;

[0189] If gsc≥osc, it means that the sediment content increases, ΔR2=[1+log255(gsc-osc)]*R2;

[0190] If gsc<osc, it means that the sediment content is reduced, ΔR2=[1-log255(osc-gsc)]*R2;

[0191] Step S268: Substitute the supplementary constraint condition into the constraint condition (step S266); execute step S266 again to construct a new dry season model and a new rainy season model;

[0192] Step S269: Determine a phytoremediation plan;

[0193] Step S2691: using the model overlap component of GIS to overlay the new dry season model and the new rainy season model;

[0194] Step S2692: marking the area where the new dry season model and the new rainy season model overlap and the color is R1 as target area 1;

[0195] The area where the new dry season model and the new rainy season model overlap and the color is ΔR2<R2 is the planting area 2;

[0196] Step S2693: in the target area 1, extracting an area not including the original area as the planting area 1;

[0197] Planting Area 1 and Planting Area 2 are used as phytoremediation solutions.

[0198] Step S3: Based on the geospatial model, the diffusion trend of vegetation and the accumulation trend of biological harmful substances are simulated according to the soil quality information, river pollution information and biodiversity information;

[0199] Based on the diffusion and accumulation trends, the protection range is divided for organisms and animal restoration plans are determined;

[0200] The specific steps of step S3 are as follows:

[0201] Step S31: Obtain the slope grid, use the Spatial Analyst tool in GIS to analyze the geospatial model, and divide the terrain of the geospatial model into flat areas, transition areas, and steep areas;

[0202] It should be noted that the "slope grid" of the present invention represents the slope range of flat areas, transition areas and steep areas; generally flat areas: 0°~5°, transition areas: 5°~25°, steep areas: <25o; the specific slope grid values ​​can be adjusted according to the actual needs of the user.

[0203] Step S32: taking the calibration range as a scanning unit, randomly selecting an area covered by vegetation in the geographic space model as a monitoring area; scanning the area around the monitoring area as a neighboring area;

[0204] It should be noted that the "calibration range" in the present invention refers to the minimum area of ​​the monitoring area and the adjacent area; generally, the calibration range is selected to be 660 cubic meters; the specific calibration range value can be adjusted according to the actual needs of the user.

[0205] Step S33: reading and calculating the soil quality information corresponding to the monitoring area, and recording the average microbial content, inorganic matter content and soil moisture content as ws, nw and hh respectively;

[0206] Read and calculate the soil quality information corresponding to the adjacent area, and record the average microbial content, inorganic matter content and soil water content as bws, bnw and bhh respectively;

[0207] Step S34: using the LCG algorithm to randomly generate weight coefficients, recorded as wa1, wa2 and wa3; wa1, wa2 and wa3, satisfy: wa1, wa2 and wa3 are all rational numbers, wa1+wa2+wa3=1, and wa1>wa2>wa3;

[0208] Step S35: Please refer to Figure 3 , defining the constraints on vegetation diffusion trends;

[0209] Constraint 31:

[0210] If wa1*bws+wa2*bnw+wa3*bhh>wa1*ws+wa2*nw+wa3*hh, then the direction of vegetation diffusion is: from the monitoring area to the adjacent area;

[0211] Constraint 32: If wa1*bws+wa2*bnw+wa3*bhh≤wa1*ws+wa2*nw+wa3*hh, then determine whether the terrain of the monitoring area satisfies: steep area → transition area → flat area;

[0212] If satisfied, the direction of vegetation diffusion is: from the monitoring area to the adjacent area;

[0213] If not, compare other neighboring areas of the monitoring area until constraint condition 31 or constraint condition 32 is satisfied;

[0214] Step S36: Based on step S35, supplement the constraint condition of the cumulative trend;

[0215] Sum up the populations of protected animals and plants, recorded as pp and pa;

[0216] Summarize the number of rivers in the area to be optimized, denoted as rn; summarize the pollutant content corresponding to each river, denoted as hn1~hn rn ;

[0217] Constraint 33: The migration direction of the distribution range of protected animals and plants is the same as the diffusion direction of vegetation;

[0218] Constraint 34: The cumulative trend of protected plants increases, and the added value is (phn) 1 / 2 / pp 2 ; Among them, phn represents the pollutant content of the river corresponding to the "distribution range of protected plants" (i.e. phn represents hn1~hn rn Any value in );

[0219] The cumulative trend of protected animals is increasing, and the added value is (ahn*pp) 1 / 2 / pa; where ahn represents the pollutant content of the river corresponding to the “distribution range of protected plants” (i.e., ahn represents hn1 to hn rn Any value in );

[0220] Step S37: Call the same steps of analyzing the plant restoration plan (steps S266 to S369), and build a cumulative model of harmful substances to animals and plants based on the geographic spatial model and constraints 31 to 34; mark the area with the least accumulated harmful substances to animals and plants in the cumulative model of harmful substances to animals and plants as the protection range, and obtain an animal restoration plan.

[0221] Step S4: Summarize the plant restoration plan and the animal restoration plan and feed back to the user.

[0222] The above formulas are all dimensionless and numerical calculations. The formula is a formula for the most recent real situation obtained by collecting a large amount of data and performing software simulation. The preset parameters in the formula are set by technical personnel in this field according to actual conditions. For example, if there are weight coefficients and proportional coefficients, their set sizes are to quantify each parameter to obtain a specific value, which is convenient for subsequent comparison. Regarding the size of the weight coefficient and the proportional coefficient, it is sufficient as long as it does not affect the proportional relationship between the parameter and the quantized value.

[0223] Finally, it should be noted that the above-described embodiments are only specific implementations of the present invention, which are used to illustrate the technical solutions of the present invention, rather than to limit them. The protection scope of the present invention is not limited thereto. Although the present invention is described in detail with reference to the above-described embodiments, ordinary technicians in the field should understand that any technician familiar with the technical field can still modify the technical solutions recorded in the above-described embodiments within the technical scope disclosed by the present invention, or can easily think of changes, or make equivalent replacements for some of the technical features therein; and these modifications, changes or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present invention, and should be included in the protection scope of the present invention. Therefore, the protection scope of the present invention shall be based on the protection scope of the claims.

Claims

1. A method for optimizing land space ecological restoration scheme based on GIS and multi-objective decision analysis, characterized in that: Optimization methods include: Step S1: Using the GIS system, obtain the terrain information of the area to be optimized and obtain a geographic space model; Obtain soil quality information, climate information, river pollution information and biodiversity information of the area to be optimized, and obtain ecological environment data; Step S2: Obtain the distribution area information and coverage area information of vegetation and rivers in the area to be optimized, and determine the vegetation change rate, water flow change rate and sediment content change rate; Using the GIS model analysis component, combined with ecological environment data, based on the geographic spatial model, simulate and analyze the soil and water loss in the area to be optimized, and determine the plant restoration plan; Step S3: Based on the geospatial model, the diffusion trend of vegetation and the accumulation trend of biological harmful substances are simulated according to the soil quality information, river pollution information and biodiversity information; Based on the diffusion and accumulation trends, the protection range is divided for organisms and animal restoration plans are determined; Step S4: Summarize the plant restoration plan and the animal restoration plan and feed back to the user; The specific steps of step S1 are as follows: Step S11: Obtain elevation points and contour lines of the area to be optimized to obtain vector data; Step S12: load the scipy library and use the Delaunay triangulation function to convert the vector data into TIN data; Step S13: Load the ArcGIS software package and use the TIN to Raster tool to convert the TIN data into a DEM to obtain a geospatial model; Step S14: Acquire soil quality information, climate information, river pollution information and biodiversity information in the area to be optimized to obtain ecological environment data; Soil quality information includes: average microbial content, inorganic content and soil moisture content of the soil in the area to be optimized; River pollution information includes: pollutant content corresponding to each river in the area to be optimized; Climate information includes: the start time and duration of the dry season in the area to be optimized, as well as the average, maximum and minimum values ​​of rainfall; The starting time and duration of the rainy season, as well as the average, maximum and minimum values ​​of rainfall; Biodiversity information includes: the population size and distribution range of all protected plants and animals in the area to be optimized; Step S15: Summarize the geographic space model and ecological environment data and proceed to step S2; The specific steps of step S2 are as follows: Step S21: using the same steps of obtaining the geographic space model, obtain the distribution area and coverage area of ​​the vegetation in the area to be optimized; the vegetation coverage area is recorded as s; Obtain the river basin, water flow and sediment content in the area to be optimized; record the water flow as f; record the sediment content as sc; Step S22: according to the distribution area of ​​vegetation, the historical vegetation coverage area of ​​the area to be optimized in the past five years is obtained, which is recorded as hs1, hs2, hs3, hs4 and hs5; According to the river basin, obtain the historical water flow and sediment content of the area to be optimized in the past five years; Historical water flows, denoted as hf1, hf2, hf3, hf4 and hf5; Historical sediment content, recorded as hsc1, hsc2, hsc3, hsc4 and hsc5; Step S23: Define calculation formula 21 to calculate the vegetation coverage area growth rate corresponding to the area to be optimized: us 1～2 、us 2～3 、us 3～4 and us 4～5 ; Water flow rate growth rate: uf 1～2 uf 2～3 uf 3～4 and uf 4～5 ; Sediment content growth rate: usc 1～2 、usc 2～3 、usc 3～4 and usc 4～5 ; Calculation 21: ; Where X represents the historical vegetation coverage area (hs1~hs5), historical water flow (hf1~hf5) and historical sediment content (hsc1~hsc5); u(X) represents the growth rate corresponding to X; Step S24: determining the vegetation change rate qs; Step S25: Repeat the same steps of determining the vegetation change rate qs to determine the change rates of water flow and sediment content, which are recorded as qf and qsc respectively; Step S26: Simulate and analyze the soil erosion in the area to be optimized and determine the plant restoration plan; The specific steps of step S26 are as follows: Step S261: Set the vegetation area as color 1 and the river area as color 2; read the RGB values ​​corresponding to color 1 and color 2, and record them as R1 and R2 respectively; Step S262: Determine the influence coefficient aa1 of the dry season on the water flow; Step S2621: the average value, the maximum value and the minimum value corresponding to the dry season rainfall are recorded as ad, md and ld respectively; Step S2622: define relation 22: lim(ad)→md; Define relation 23: lim(ad)→ld; Step S2623: If definition 22 is satisfied, then the influence coefficient aa1, aa1 = (ad + md) / 2; If definition 23 is satisfied, then the influence coefficient aa1, aa1 = (ad + ld) / 2; If neither definition 22 nor 23 is satisfied, then the influence coefficient is aa1, aa1 = (md / ad + ad / ld) / 2; Step S263: Repeat the same steps of determining the influence coefficient aa1 to determine the influence coefficient aa2 of the rainy season on the water flow; Step S264: constructing a relational equation 1 according to the historical vegetation coverage area (hs1-hs5) and the historical sand content (hsc1-hsc5): y=a1x1+b1; wherein y represents the sand content; x1 represents the vegetation coverage area; a1 and b1 represent the slope and intercept of the relational equation 1 respectively; Step S265: Repeat the same steps of constructing relational equation 1, and construct relational equation 2 according to the historical water flow (hf1-hf5) and the historical sediment content (hsc1-hsc5): y=a2x2+b2; wherein x2 represents the water flow, a2 and b2 represent the slope and intercept of relational equation 2 respectively; The subsequent steps of step S265 are as follows: Step S266: taking the geographic space model as the model benchmark, taking s, f and sc as the basic parameters, taking the vegetation change rate qs and the influence coefficient aa1 as the change parameters, and taking the relationship equation 1 and the relationship equation 2 as the constraint conditions; Taking the start time and duration of the dry season as simulation period 1, the CSLE generation component in GIS is called to generate the dry season model; Taking the start time and duration of the rainy season as simulation period 2, the CSLE generation component in GIS is called to generate the rainy season model; Step S267: Determine the vegetation diffusion direction and R2 change relationship in the rainy season model or the dry season model as a supplementary constraint condition; Step S268: Substitute the supplementary constraint condition into the constraint condition of step S266; execute step S266 again to construct a new dry season model and a new rainy season model; Step S269: Determine a phytoremediation plan; Step S2691: using the model overlap component of DIS to overlay the new dry season model and the new rainy season model; Step S2692: marking the area where the new dry season model and the new rainy season model overlap and the color is R1 as target area 1; The area where the new dry season model and the new rainy season model overlap and the color is ΔR2<R2 is the planting area 2; Step S2693: in the target area 1, extracting an area not including the original area as the planting area 1; Planting Area 1 and Planting Area 2 are used as phytoremediation solutions; The specific steps of step S267 are as follows: Step S2671: Mark the vegetation-covered area in the geographic space model as the original area; read the water flow of the river that intersects with the original area and has the shortest spatial distance as the original flow; Step S2672: repeat the step of reading the original area and the original flow in step S2671, and in the rainy season model or the dry season model, read the adjacent vegetation coverage area that does not intersect with the original area as the adjacent area, and read the river water flow corresponding to the adjacent area as the adjacent flow; In the rainy season model or dry season model, the direction of vegetation diffusion is: from the original area to the adjacent area where the adjacent flow rate is ≥ the original flow rate, and the vegetation area after diffusion is marked by R1; Step S2673: Mark the river area in the geographic space model as the target river; read the sediment content of the target river as the original sediment content, recorded as osc; Read the sediment content corresponding to the target river in the rainy season model or the dry season model as the new sediment content, recorded as gsc; Step S2674: Compare the magnitudes of osc and gsc to determine the R2 variation relationship ΔR2; If gsc≥osc, it means that the sediment content increases, ΔR2=[1+log255(gsc-osc)]*R2; If gsc<osc, it means that the sediment content is reduced, ΔR2=[1-log255(osc-gsc)]*R2; The specific steps of step S3 are as follows: Step S31: Obtain the slope grid, use the Spatial Analyst tool in GIS to analyze the geospatial model, and divide the terrain of the geospatial model into flat areas, transition areas, and steep areas; Step S32: taking the calibration range as a scanning unit, randomly selecting an area covered by vegetation in the geographic space model as a monitoring area; scanning the area around the monitoring area as a neighboring area; Step S33: reading and calculating the soil quality information corresponding to the monitoring area, and recording the average microbial content, inorganic matter content and soil moisture content as ws, nw and hh respectively; Read and calculate the soil quality information corresponding to the adjacent area, and record the average microbial content, inorganic matter content and soil water content as bws, bnw and bhh respectively; Step S34: using the LCG algorithm to randomly generate weight coefficients, recorded as wa1, wa2 and wa3; wa1, wa2 and wa3, satisfy: wa1, wa2 and wa3 are all rational numbers, wa1 + wa2 + wa3 = 1, and wa1> wa2> wa3; Step S35: defining the constraint conditions of vegetation diffusion trend; Constraint 31: If wa1*bws+wa2*bnw+wa3*bhh>wa1*ws+wa2*nw+wa3*hh, then the direction of vegetation diffusion is: from the monitoring area to the adjacent area; Constraint 32: If wa1*bws+wa2*bnw+wa3*bhh≤wa1*ws+wa2*nw+wa3*hh, then determine whether the terrain of the monitoring area satisfies: steep area → transition area → flat area; If satisfied, the direction of vegetation diffusion is: from the monitoring area to the adjacent area; If not, compare other neighboring areas of the monitoring area until constraint condition 31 or constraint condition 32 is satisfied; Step S36: Based on step S35, supplement the constraint condition of the cumulative trend; Sum up the populations of protected animals and plants, recorded as pp and pa; Summarize the number of rivers in the area to be optimized, recorded as rn; summarize the pollutant content corresponding to each river, recorded as hn1~hn rn ; Constraint 33: The migration direction of the distribution range of protected animals and plants is the same as the diffusion direction of vegetation; Constraint 34: The cumulative trend of protected plants increases, and the added value is (phn) 1 / 2 / pp 2 ; Among them, phn represents the pollutant content in the river corresponding to the distribution range of protected plants; The cumulative trend of protected animals is increasing, and the added value is (ahn*pp) 1 / 2 / pa; where ahn represents the pollutant content in the river corresponding to the distribution range of protected animals; Step S37: Call the same steps of analyzing the plant restoration plan, and build a cumulative model of harmful substances to animals and plants based on the geographic space model and constraints 31 to 34; mark the area with the least accumulated harmful substances to animals and plants in the cumulative model of harmful substances to animals and plants as the protection range, and obtain an animal restoration plan.

2. According to claim 1, a method for optimizing land space ecological restoration scheme based on GIS and multi-objective decision analysis is characterized in that: The specific steps of step S264 are as follows: Step S2641: Import the LinearRegression class from the sklearn.linear_model module, and then import the NumPy library; Step S2642: create a receiving array 1 for storing hs1 to hs5; create a receiving array 2 for storing hsc1 to hsc5; Step S2643: Create a receiving object model instantiated by the LinearRegression class; Step S2644: calling the fit method through model, and substituting the receiving array 1 and the receiving array 2 as parameters into the fit method, obtaining the relation equation 1: y=a0x1+b0; wherein y represents the sand content; x1 represents the vegetation coverage area; a0 and b0 represent the slope and intercept respectively, and the initial values ​​of a0 and b0 are both 0; Step S2645: Call the coef method and the intercept method through the model to calculate a0 and b0 respectively, and obtain a1 and b1; where a1 is a rational number and is not equal to 0; b1 is an arbitrary rational number; Step S2646: Summarize the processing results of steps S2641 to S2645 to obtain relationship equation 1: y=a1x1+b1.

Images