Simulation method for non-point source pollution along the way in coastal saline-alkali land irrigation area based on the network topology of drainage ditch

Abstract

The application discloses a simulation method for non-point source pollution along the way in a coastal saline-alkali land irrigation area based on a drainage ditch network topology, and belongs to the technical field of agricultural non-point source pollution control and water environment management. The method comprises the following steps: collecting basic data of the irrigation area, correcting a DEM, constructing a drainage ditch network topology model, and calculating attenuation coefficients and dynamic flow thresholds of each directed edge; obtaining nitrogen and phosphorus attenuation coefficients and salt transmission coefficients; traversing the migration path of a farmland unit to a water outlet to obtain the length of each directed edge; calculating the initial load of nitrogen, phosphorus and salt in the farmland unit, combining the path length to obtain the river load, and finally summarizing the river load of the whole irrigation area. The method can accurately identify the confluence path in a plain area with a slope of less than 0.1 degree, dynamically adjust the attenuation coefficient according to the path difference, and improve the simulation accuracy.

Description

Technical Field

[0001] This invention belongs to the field of agricultural non-point source pollution control and water environment management technology, specifically involving a method for simulating the process-by-process absorption of non-point source pollution in coastal saline-alkali irrigation areas based on drainage ditch network topology. Background Technology

[0002] Coastal saline-alkali land serves as a reserve of arable land and is of great significance for ensuring food security. When developing coastal saline-alkali land, if large-scale irrigation is used to leach salt, the loss of salt will be accompanied by an increase in nitrogen and phosphorus loss; on the other hand, water-saving and pollution-reduction measures can easily lead to salt accumulation.

[0003] Irrigation district drainage ditch systems serve a dual function: they are both primary channels for nitrogen, phosphorus, and salt migration and important sites for their natural absorption. Calculating the absorption effect of drainage ditches on nitrogen, phosphorus, and salt can effectively assess the merits of different development methods for coastal saline-alkali land. However, current terrain extraction algorithms still have the following shortcomings: 1. It is difficult to identify the true confluence path in plains areas with a slope of less than 0.1 degrees; 2. The impact of different migration paths on the absorption capacity was ignored.

[0004] Therefore, it is necessary to propose a simulation method for the along-path absorption of non-point source pollution in coastal saline-alkali irrigation areas based on drainage ditch network topology in order to solve the above problems. Summary of the Invention

[0005] This invention proposes a method for simulating the non-point source pollution mitigation along the route in coastal saline-alkali irrigation areas based on drainage ditch network topology, in order to solve the problems of low identification accuracy of traditional models in flat irrigation areas and the impact of migration paths on mitigation capacity.

[0006] To achieve the above objectives, the present invention proposes the following technical content: A simulation method for non-point source pollution mitigation along the irrigation route in coastal saline-alkali land irrigation areas based on drainage ditch network topology includes the following steps: S1: Collect basic data of the target irrigation area and perform local DEM correction preprocessing to obtain the corrected DEM data of the target irrigation area; S2: Based on the corrected target irrigation area DEM data in S1, construct a drainage ditch network topology model and calculate the attenuation coefficient and dynamic flow threshold of each directed edge in the drainage ditch network topology model. S3: Combine the attenuation coefficient and dynamic flow threshold of each directed edge to calculate the nitrogen absorption coefficient, phosphorus absorption coefficient and salt transport coefficient of each directed edge. S4: Based on the drainage ditch network topology model established in S2, traverse the unique migration path from each farmland unit to its outlet to obtain the length of each directed edge on the migration path of each farmland unit. S5: Calculate the initial nitrogen, phosphorus and salinity loads for each farmland unit; S6: Based on the initial nitrogen, phosphorus and salinity loads of each farmland unit and the length of each directed edge on the migration path of each farmland unit, calculate the nitrogen, phosphorus and salinity loads entering the river for each farmland unit. S7: Calculate the nitrogen, phosphorus, and salt inflow loads into the river for all farmland units in the target irrigation district.

[0007] Furthermore, in step S3, the first k The formula for calculating the dynamic flow threshold of a directed edge is:

[0008] In the formula, S field，k Indicates the first k The catchment area of ​​a directed edge; 0.0015 represents the corrected drainage coefficient after extreme rainfall; Q edge,k Indicates the first k A dynamic flow threshold for each directed edge is used to constrain the flow of each directed edge, as described below. No. k The formula for calculating the attenuation coefficient of a directed edge is:

[0009] In the formula, L edge，k Indicates the first k The length of each directed edge; W k Indicates the first k The average width of the water surface on each directed edge; c k Indicates the first k Measured correction factor for a directed edge; K edge,k Indicates the first k The attenuation coefficient of a directed edge.

[0010] Furthermore, in step S3, the first k The nitrogen absorption coefficient of a directed edge is given by the formula:

[0011] In the formula, K N，k Indicates the first k The nitrogen absorption coefficient of a directed edge; K edge,k Indicates the first k The attenuation coefficient of a directed edge; V f,k Indicates the first kThe denitrification rate per unit area of ​​a directed edge; W k Indicates the first k The average width of the water surface on each directed edge; Q k Indicates the first k The flow of a directed edge is less than or equal to Q edge,k ; Indicates the first k Localized comprehensive correction parameters for each directed edge; K N,k Indicates the first k The nitrogen absorption coefficient of a directed edge; No. k The phosphorus absorption coefficient for each directed edge is given by the following formula:

[0012] In the formula, K p,k Indicates the first k The phosphorus absorption coefficient of a directed edge; K edge,k Indicates the first k The attenuation coefficient of a directed edge; α p,k and β p,k Indicates the first k The localization parameters of a directed edge; T r,k Indicates the first k Hydraulic residence time of a directed edge; No. k The salt transport coefficient of a directed edge is given by the formula:

[0013] In the formula, K s,k Indicates the first k Salt transport coefficient of a directed edge; K edge,k Indicates the first k The attenuation coefficient of a directed edge; E k Indicates the first k The evaporation intensity of a directed edge of water surface; W k Indicates the first k The average width of the water surface on each directed edge; r Indicates the density of water; Q k Indicates the first k The flow of a directed edge; This represents the soil salt release coefficient on the slope.

[0014] Furthermore, step S5 includes the following steps: S5.1: Calculate the first x The initial nitrogen load for each farmland unit is calculated using the following formula:

[0015] In the formula, L N,x Indicates the first x Initial nitrogen load of each farmland unit; A x Indicates the first x The area of ​​each farmland unit; F N,x Indicates the first x Nitrogen fertilizer application rate per unit area of ​​each farmland unit; m N Represents the empirical coefficient; C season Indicates crop growth period or hydrological scenario correction factor; S5.2: Calculate the first x The initial phosphorus load for each farmland unit is calculated using the following formula:

[0016] In the formula, A x Indicates the first x The area of ​​each farmland unit; F P,x Indicates the first x Phosphate fertilizer application rate per unit area of ​​each farmland unit; m P Indicates the phosphorus loss rate; C clay This represents the correction factor for soil clay content; L P,x Indicates the first x Initial phosphorus load of each farmland unit; S5.3: Calculate the first x The initial salinity load of each farmland unit is calculated using the following formula:

[0017] In the formula ，Q drain,x Indicates the first x Drainage of each farmland unit; C salt,x Indicates the first x Initial salinity concentration of each farmland unit; Indicates the salt concentration coefficient; S 0,x Indicates the first xInitial salinity load of each farmland unit; A x Indicates the first x The area of ​​each farmland unit; P Indicates rainfall amount; I rr,x Indicates the first x Irrigation volume for each farmland unit; Indicates the first x The runoff and drainage coefficients of each farmland unit.

[0018] Furthermore, step S6 includes the following steps: S6.1: Calculate the first... x The nitrogen inflow load of each farmland unit into the river is calculated using the following formula:

[0019] In the formula, L N,out,x Indicates the first x Nitrogen load into the river in each farmland unit; L N,x Indicates the first x Initial nitrogen load of each farmland unit; D k This represents the length of the k-th directed edge; K N，k Indicates the first k The nitrogen absorption coefficient of a directed edge; Z x Indicates the first x The number of directed edges covered by the path of each farmland unit; S6.2: Calculate the first... x The phosphorus inflow load of each farmland unit into the river is calculated using the following formula:

[0020] In the formula, L P,out,x Indicates the first x Phosphorus load into the river in each farmland unit; L P,x Indicates the first x Initial phosphorus load of each farmland unit; D k This represents the length of the k-th directed edge; K p,k Indicates the first k The phosphorus absorption coefficient of a directed edge; Z x Indicates the first x The number of directed edges covered by the path of each farmland unit; S6.3: Calculate the first x The formula for the salinity load into the river for each farmland unit is:

[0021] In the formula, S out,x Indicates the first x The salt load of each farmland unit entering the river; S 0,x Indicates the first x Initial salinity load of each farmland unit; D k This represents the length of the k-th directed edge; K s,k Indicates the first k Salt transport coefficient of a directed edge; Z x Indicates the first x The number of directed edges covered by the path of a farmland unit.

[0022] Furthermore, step S7 includes the following steps: S7.1: Nitrogen inflow load into the river for all farmland units in the target irrigation district, calculated using the following formula:

[0023] In the formula, L N,out,x Indicates the first x Nitrogen load into the river in each farmland unit; Ng Indicates the number of farmland units in the target irrigation district; TL N,out,x This represents the nitrogen inflow load into the river for all farmland units in the target irrigation district; S7.2: Phosphorus inflow load into the river for all farmland units in the target irrigation area, calculated using the following formula:

[0024] In the formula, L P,out,x Indicates the first x Phosphorus load into the river in each farmland unit; Ng Indicates the number of farmland units in the target irrigation district; TL P,out,x This represents the phosphorus inflow load into the river for all farmland units in the target irrigation district; S7.3: Salt inflow load into the river for all farmland units in the target irrigation district, calculated using the following formula:

[0025] In the formula, S out,x Indicates the first x The salt load of each farmland unit entering the river; TS out This represents the salt load into the river for all farmland units in the target irrigation area; Ng This indicates the number of farmland units in the target irrigation district.

[0026] Furthermore, step S1 includes the following steps: S1.1: Obtain basic data for the target irrigation area; S1.2: The DEM local correction preprocessing algorithm is used to rasterize the drainage ditch vector data to generate drainage ditch raster data with the same spatial resolution as the original DEM data. S1.3: Based on the rasterized drainage ditch data, obtain the "downcut depth" value of each ditch pixel, and perform pixel-by-pixel correction on the DEM data of the target irrigation area. The formula is as follows:

[0027] In the formula, H corrected ( i,j ) indicates position ( i,j The corrected elevation value at point ) H original ( i,j ) indicates position ( i, j The elevation value before correction at point ) D ( i,j ) indicates position ( i,j The incision depth of the ditch pixel at point ).

[0028] The beneficial effects that can be achieved by adopting the above technologies are: 1. The problem of difficult identification of confluence paths under flat terrain was solved by local DEM correction, ensuring the accuracy of the spatial basis of the simulation.

[0029] 2. Based on directed graph networks, the entire path from farmland to outlet is tracked, enabling the simulation results to accurately reflect the spatiotemporal dynamic characteristics of water, salt, and pollutants in complex drainage systems. 3. Construct the absorption coefficient for each directed edge to scientifically characterize the absorption capacity of ditches of different sizes. Attached Figure Description

[0030] Figure 1 This is the flowchart of this method. Detailed Implementation

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

[0032] like Figure 1 As shown, the simulation method for non-point source pollution mitigation along the irrigation route in coastal saline-alkali land irrigation areas based on drainage ditch network topology includes the following steps: S1: Collect basic data of the target irrigation area and perform local DEM correction preprocessing. This includes the following steps: S1.1: Obtain basic data for the target irrigation area.

[0033] Specifically, the basic data of the target irrigation area (such as the coastal irrigation area in Yancheng, Jiangsu) is obtained, including: original digital elevation model (DEM) data, vector data of drainage ditches (including spatial location, ditch grade, cross-sectional dimensions and other attribute information), and spatial information of farmland units.

[0034] S1.2: The DEM local correction preprocessing algorithm is used to rasterize the drainage ditch vector data to generate drainage ditch raster data with the same spatial resolution as the original DEM data.

[0035] S1.3: Based on the rasterized drainage ditch data, obtain the "downcut depth" value of each ditch pixel, and perform pixel-by-pixel correction on the DEM data of the target irrigation area. The formula is as follows:

[0036] In equation (1), H corrected ( i,j ) indicates position ( i,j The corrected elevation value at point ) H original ( i,j ) indicates position ( i,j The elevation value before correction at point ) D ( i,j ) indicates position ( i,j The incision depth of the ditch pixel at point ).

[0037] S2: Based on the corrected target irrigation district DEM data in S1, construct a drainage ditch network topology model. This includes the following steps: S2.1: Define nodes and attributes.

[0038] Define the centroid of each farmland unit as the source node; Key locations such as the intersections and inflection points of ditches are defined as transmission nodes; The ecological wastewater treatment pond, the forced drainage pumping station, and the final outlet are defined as the treatment node, the control node, and the collection node, respectively.

[0039] The attributes of each node are as follows: Source nodes: area, amount of fertilizer applied, soil salinity.

[0040] Transmission nodes: slope, width, length.

[0041] Processing nodes: nitrogen absorption efficiency and phosphorus absorption efficiency.

[0042] S2.2: Define directed edges and attributes.

[0043] Based on the water flow direction and engineering connections, directed edges are constructed using a depth-first search algorithm. The construction logic for the directed edges is: source node → transmission node → processing node → control node → sink node.

[0044] Each directed edge has basic attributes (length, width, type) and dynamic attributes (dynamic flow threshold, attenuation coefficient).

[0045] In this scheme, the types of directed edges are: trunk, branch, and ditch; each directed edge belongs to one of the trunk, branch, and ditch types.

[0046] When nodes of two types are close to each other, determine whether a directed edge exists using the following method: First, calculate the Euclidean distance between the two nodes (A and B), and then calculate the distance matching degree based on the preset maximum search distance between the two nodes, using the following formula:

[0047] In the formula, R dist,AB This represents the distance matching degree between nodes A and B; distance AB This represents the Euclidean distance between nodes A and B. distance max This represents the maximum search distance between two preset nodes; Then, based on the slope and width of node B, its flow fit is calculated using the following formula:

[0048] In the formula, slope B This represents the normalized slope value of node B; width B This represents the normalized width value of node B; R flow,AB Indicates the water flow fit between nodes A and B; Next, query whether the two nodes (A and B) are within the service area of ​​the same pumping station. If they are, then the service matching degree of the two nodes is determined. R service,AB The service matching score is 1 if it is not present. R service,AB =0; Finally, the spatial correlation degree between two nodes is calculated using the following formula:

[0049] like R AB Greater than the set threshold R 阈 If the value is 0.5, then node A has a directed edge to node B; otherwise, node A does not have a directed edge to node B.

[0050] S2.3: Calculate the dynamic properties of each directed edge.

[0051] Among them, the k The formula for calculating the dynamic flow threshold of a directed edge is:

[0052] In equation (2), S field，k Indicates the first k The catchment area of ​​a directed edge; 0.0015 represents the corrected drainage coefficient after extreme rainfall; Q edge,k Indicates the first k A dynamic flow threshold for each directed edge is used to constrain the flow of each directed edge, as described later.

[0053] No. k The formula for calculating the attenuation coefficient of a directed edge is:

[0054] In equation (3), L edge，k Indicates the first k The length of each directed edge; W k Indicates the first k The average width of the water surface on each directed edge; c k Indicates the first k Measured correction factor for a directed edge; K edge,k Indicates the first k The attenuation coefficient of a directed edge.

[0055] S3: Calculate the nitrogen absorption coefficient, phosphorus absorption coefficient, and salt transport coefficient for each directed edge.

[0056] Among them, the k The nitrogen absorption coefficient of a directed edge, based on denitrification kinetics, is given by the following formula:

[0057] In the formula, K N，kIndicates the first k The nitrogen absorption coefficient of a directed edge; K edge,k Indicates the first k The attenuation coefficient of each directed edge (calculated from S2.3); V f,k Indicates the first k The denitrification rate per unit area of ​​a directed edge; W k Indicates the first k The average width of the water surface on each directed edge; Q k Indicates the first k The flow of a directed edge is less than or equal to Q edge,k ; Indicates the first k Localized comprehensive correction parameters for each directed edge; K N,k Indicates the first k The nitrogen absorption coefficient of a directed edge.

[0058] No. k The phosphorus absorption coefficient of a directed edge, based on deposition adsorption kinetics, is given by the following formula:

[0059] In equation (5), K p,k Indicates the first k The phosphorus absorption coefficient of a directed edge; K edge,k Indicates the first k The attenuation coefficient of each directed edge (calculated from S2.3); α p,k and β p,k Indicates the first k The localization parameters of a directed edge; T r,k Indicates the first k Hydraulic residence time of a directed edge.

[0060] No. k The salt transport coefficient of a directed edge is given by the formula:

[0061] In equation (6), K s,k Indicates the first k Salt transport coefficient of a directed edge; K edge,k Indicates the first k The attenuation coefficient of each directed edge (calculated from S2.3); Ek Indicates the first k The evaporation intensity of a directed edge of water surface; W k Indicates the first k The average width of the water surface on each directed edge; r Indicates the density of water; Q k Indicates the first k The flow of a directed edge; This represents the soil salt release coefficient on the slope.

[0062] S4: Based on the drainage ditch network topology model established in S2, traverse each farmland unit (i.e., the source node) and obtain the unique migration path to its outlet, thus obtaining the length of each directed edge on the migration path of each farmland unit.

[0063] S5: Calculate the initial pollution load and initial salinity load for each farmland unit. This includes the following steps: S5.1: Calculate the first x The initial nitrogen load for each farmland unit is calculated using the following formula:

[0064] In equation (7), L N,x Indicates the first x Initial nitrogen load of each farmland unit; A x Indicates the first x The area of ​​each farmland unit; F N,x Indicates the first x Nitrogen fertilizer application rate per unit area of ​​each farmland unit; m N Represents the empirical coefficient; C season This indicates a crop growth period or hydrological scenario correction factor.

[0065] S5.2: Calculate the first x The initial phosphorus load for each farmland unit is calculated using the following formula:

[0066] In equation (8), A x Indicates the first x The area of ​​each farmland unit; F P,x Indicates the first x Phosphate fertilizer application rate per unit area of ​​each farmland unit; m P Indicates the phosphorus loss rate; C clay This represents the correction factor for soil clay content;L P,x Indicates the first x Initial phosphorus load of each farmland unit.

[0067] S5.3: Calculate the first x The initial salinity load of each farmland unit is calculated using the following formula:

[0068] In equation (9), Q drain,x Indicates the first x Drainage of each farmland unit; C salt,x Indicates the first x Initial salinity concentration of each farmland unit; Indicates the salt concentration coefficient; S 0,x Indicates the first x Initial salinity load of each farmland unit.

[0069]

[0070] In equation (10), A x Indicates the first x The area of ​​each farmland unit; P Indicates rainfall amount; I rr,x Indicates the first x Irrigation volume for each farmland unit; Indicates the first x The runoff and drainage coefficients of each farmland unit.

[0071] S6: Calculate the inflow load for each farmland unit. This includes the following steps: S6.1: Calculate the first... x The nitrogen load into the river for each farmland unit follows first-order kinetics. The formula is:

[0072] In equation (11), L N,out,x Indicates the first x Nitrogen load into the river in each farmland unit; L N,x Indicates the first x Initial nitrogen load of each farmland unit; D k This represents the length of the k-th directed edge; K N，k Indicates the first k The nitrogen absorption coefficient of a directed edge; Z x Indicates the firstx The number of directed edges covered by the path of a farmland unit.

[0073] S6.2: Calculate the first... x The phosphorus inflow load of each farmland unit follows first-order kinetics. The formula is:

[0074] In equation (12), L P,out,x Indicates the first x Phosphorus load into the river in each farmland unit; L P,x Indicates the first x Initial phosphorus load of each farmland unit; D k This represents the length of the k-th directed edge; K p,k Indicates the first k The phosphorus absorption coefficient of a directed edge; Z x Indicates the first x The number of directed edges covered by the path of a farmland unit.

[0075] S6.3: Calculate the first x The salt load entering the river for each farmland unit is calculated using a headway transport model. The formula is:

[0076] In equation (13), S out,x Indicates the first x The salt load of each farmland unit entering the river; S 0,x Indicates the first x Initial salinity load of each farmland unit; D k This represents the length of the k-th directed edge; K s,k Indicates the first k Salt transport coefficient of a directed edge; Z x Indicates the first x The number of directed edges covered by the path of a farmland unit.

[0077] S7: Calculate the nitrogen, phosphorus, and salinity loads entering the river for all farmland units in the target irrigation district. This includes the following steps: S7.1: Nitrogen inflow load into the river for all farmland units in the target irrigation district, calculated using the following formula:

[0078] In equation (14), L N,out,xIndicates the first x Nitrogen load into the river in each farmland unit; Ng Indicates the number of farmland units in the target irrigation district; TL N,out,x This represents the nitrogen inflow load into the river for all farmland units in the target irrigation district.

[0079] S7.2: Phosphorus inflow load into the river for all farmland units in the target irrigation area, calculated using the following formula:

[0080] In equation (15), L P,out,x Indicates the first x Phosphorus load into the river in each farmland unit; Ng Indicates the number of farmland units in the target irrigation district; TL P,out,x This represents the phosphorus inflow load into the river for all farmland units in the target irrigation district.

[0081] S7.3: Salt inflow load into the river for all farmland units in the target irrigation district, calculated using the following formula:

[0082] In equation (16), S out,x Indicates the first x The salt load of each farmland unit entering the river; TS out This represents the salt load into the river for all farmland units in the target irrigation area; Ng This indicates the number of farmland units in the target irrigation district.

[0083] Calculation example: Taking a coastal irrigation area in Yancheng, Jiangsu Province as an example, the irrigation area has a flat terrain with a slope of less than 0.08 degrees. A typical irrigation ditch system serving 6 farmland units is selected, with each farmland unit having an area of ​​0.5 hectares. When constructing the drainage ditch network topology model, the centroids of the 6 farmland units are defined as source nodes S1 to S6, the intersection of the upstream and midstream sections of the irrigation ditch are defined as transmission nodes T1 and T2, respectively, and the outlet at the end of the irrigation ditch is defined as the collection node O1.

[0084] Based on the direction of water flow, three directed edges are defined. The first edge, Edge1, runs from the source node to the transmission node T1, with a length of 150 meters and a width of 1.5 meters. The second edge, Edge2, runs from the transmission node T1 to the transmission node T2, with a length of 100 meters and a width of 1.2 meters. The third edge, Edge3, runs from the transmission node T2 to the sink node O1, with a length of 50 meters and a width of 1.0 meter.

[0085] When calculating the dynamic flow threshold, according to the formula Q edge,k =S field,k × α The corrected drainage coefficient is taken as 0.0015 cubic meters per second per square meter. Substituting this into the farmland area of ​​5000 square meters, the flow threshold for all three directed edges is calculated to be 7.5 cubic meters per second. When calculating the ditch attenuation coefficient, the formula is used... K edge,k = c k ×( L edge,k / W k ), where the correction factor c k Based on the vegetation cover of the gully, 85% vegetation cover in the upstream section corresponds to c1 equal to 0.025, 78% vegetation cover in the midstream section corresponds to c2 equal to 0.022, and 72% vegetation cover in the downstream section corresponds to c3 equal to 0.018. Substituting the length and width of each side, the following values ​​are calculated: K edge,1 =2.500、 K edge,2 =1.833、 K edge,3 =0.900. This embodiment demonstrates that the proposed method can accurately construct pollutant migration paths based on real ditch network topology, avoiding the difficulty in identifying confluence paths caused by insignificant elevation differences in flat areas using traditional DEM methods.

[0086] The nitrogen absorption coefficients for each section were calculated based on field measurement data. The actual flow rates at the three cross-sections were measured using a portable flow meter. Q 1 = 0.26 liters per second Q 2 = 0.19 liters per second Q 3 = 0.12 liters per second. The vegetation nitrogen uptake rate was obtained through quadrat survey and biomass measurement. V f,1 =1.5 grams per square meter per day V f,2 =1.8 grams per square meter per day V f,3 =2.0 grams per square meter per day. Based on water quality analysis and isotope tracer experiments, the nitrogen conversion efficiency coefficients were determined to be... f N,1 =0.04、 f N,2 =0.05、 f N,3 =0.10. According to the formula K N,k = K edge,k × f N,k ×(V f,k ×W k ) / Q k Perform the calculation; for the first edge, K edge,1 =2.500、 f N,1 =0.04、 V f,1 =1.5、 W 1 = 1.5 Q Substituting 1=0.26, we get... K N,1 =2.500×0.04×(1.5×1.5) / 0.26=0.865, and similarly calculated to obtain K N,2 =1.042、 K N,3 =1.500.

[0087] For the phosphorus absorption coefficient, according to the formula K P,k = K edge,k × α p,k ×(1-e^(- β p,k × T r,k The adsorption rate constant was obtained by actual measurement. α p,1 =0.035、 α p,2 =0.042、 α p,3 =0.038, adsorption kinetic parameters β p,1 =0.18、 β p,2 =0.22、 β p,3 =0.25, hydraulic residence time T r,1 =2.5 hours T r,2 =3.2 hours T r,3 =4.1 hours, substituting into the calculation yields K P,1 =0.286、 K P,2 =0.348、 K P,3 =0.228.

[0088] For the salt migration coefficient, according to the formula K s,k = K edge,k ×[1+( E k × W k × p / Q k )×(1+ c soil [The measured evapotranspiration rate in June] E =5.2 mm per day, soil salinity release coefficient c soil =0.08, and the corresponding salt migration coefficient is obtained by substituting it into the calculation.

[0089] This embodiment demonstrates that the absorption coefficient calculation system established by this method can comprehensively consider multiple absorption mechanisms such as vegetation absorption, microbial action, and physical adsorption, and ensure the reliability of the calculation results through calibration with measured parameters.

[0090] The initial nitrogen load and the load entering the river were calculated and compared with the measured data for verification.

[0091] According to the formula L N,x = A x × F N,x × m N × C season Calculate the initial nitrogen load, including the farmland area. A x =5000 square meters, nitrogen application rate F N,x =800 kg / ha, nitrogen leaching coefficient m N =0.35 is a seasonal correction factor for June, determined based on soil leaching tests. C June =1.20 is calculated based on a rainfall of 286.6 mm and an irrigation amount of 250 mm. Substituting these values ​​into the formula yields... L N,x,June =5000×(800 / 10000)×0.35×1.20=168 kg. When calculating the pollution along the treatment path, the pollutants travel from the farmland unit through three directed edges to the outlet. The treatment index is = K N,1 × D 1+ K N,2 × D 2+ K N,3 ×D 3 = 0.00865 × 150 + 0.01042 × 100 + 0.01500 × 50 = 3.090, attenuation factor = e^(-3.090) = 0.0456, absorption rate = 1 - 0.0456 = 95.44%, and the absorption rate after correction based on measured data is 97.1%. According to the formula... L N,out,x = L N,x The nitrogen load entering the river for a single farmland unit is calculated as 168 × 0.029 = 4.872 kg, and the total nitrogen load entering the river for 6 farmland units is 4.872 × 6 = 29.23 kg.

[0092] According to the formula Q drain,x = A x ×( P + I rr,x )×Ψ x Calculate the drainage volume, where the drainage coefficient Ψ x =0.65, the calculated drainage volume of a single unit in June is 1743.95 cubic meters, the total drainage volume of 6 units is 10463.7 cubic meters, and the TN concentration entering the river is 29230 grams / 10463.7 cubic meters = 2.79 milligrams per liter.

[0093] Water quality monitoring was conducted at the outlet for four consecutive months. The measured TN concentrations were 2.65 mg / L in June, 5.38 mg / L in July, 2.48 mg / L in August, and 3.62 mg / L in September. Compared with the simulated values ​​of 2.79, 5.15, 2.60, and 3.48 mg / L, the relative errors were +5.3%, -4.3%, +4.8%, and -3.9%, respectively, with an average relative error of ±4.6%. The coefficient of determination R² reached 0.94, and the Nash efficiency coefficient NSE was 0.91, meeting the accuracy requirements for engineering applications.

[0094] Simultaneously, a comparative simulation using the traditional SWAT model was conducted, with a TN concentration error of ±23.5%, indicating that the accuracy of this method is 80% higher than that of the traditional DEM method. Four months of continuous monitoring quantified the absorption capacity of the agricultural ditch system. Of the total 690.6 kg of nitrogen entering the field, 570.51 kg was absorbed by the agricultural ditch system, achieving an absorption rate of 82.6%. Vegetation absorption contributed 45%, microbial denitrification 32%, and sediment adsorption 23%, demonstrating the important role of the agricultural ditch ecosystem in reducing non-point source pollution and providing a scientific basis for water environment protection in coastal saline-alkali irrigation areas.

[0095] Based on the above-described preferred embodiments of the present invention, and through the foregoing description, those skilled in the art can make various changes and modifications without departing from the inventive concept. The technical scope of this invention is not limited to the contents of the specification, but must be determined according to the scope of the claims.

Claims

1. A method for simulating the diffusion of non-point source pollution along the irrigation route in coastal saline-alkali land irrigation areas based on drainage ditch network topology, characterized in that, Includes the following steps: S1: Collect basic data of the target irrigation area and perform local DEM correction preprocessing to obtain the corrected DEM data of the target irrigation area; S2: Based on the corrected target irrigation area DEM data in S1, construct a drainage ditch network topology model and calculate the attenuation coefficient and dynamic flow threshold of each directed edge in the drainage ditch network topology model. S3: Combine the attenuation coefficient and dynamic flow threshold of each directed edge to calculate the nitrogen absorption coefficient, phosphorus absorption coefficient and salt transport coefficient of each directed edge. S4: Based on the drainage ditch network topology model established in S2, traverse the unique migration path from each farmland unit to its outlet to obtain the length of each directed edge on the migration path of each farmland unit. S5: Calculate the initial nitrogen, phosphorus and salinity loads for each farmland unit; S6: Based on the initial nitrogen, phosphorus and salinity loads of each farmland unit and the length of each directed edge on the migration path of each farmland unit, calculate the nitrogen, phosphorus and salinity loads entering the river for each farmland unit. S7: Calculate the nitrogen, phosphorus, and salt inflow loads into the river for all farmland units in the target irrigation district.

2. The method for simulating the diffusion of non-point source pollution along the irrigation route in coastal saline-alkali land irrigation areas based on drainage ditch network topology according to claim 1, characterized in that, In step S3, the first k The formula for calculating the dynamic flow threshold of a directed edge is: ； In the formula, S field，k Indicates the first k The catchment area of ​​a directed edge; 0.0015 represents the corrected drainage coefficient after extreme rainfall; Q edge,k Indicates the first k A dynamic flow threshold for each directed edge is used to constrain the flow of each directed edge, as described below. No. k The formula for calculating the attenuation coefficient of a directed edge is: ； In the formula, L edge，k Indicates the first k The length of each directed edge; W k Indicates the first k The average width of the water surface on each directed edge; c k Indicates the first k Measured correction factor for a directed edge; K edge,k Indicates the first k The attenuation coefficient of a directed edge.

3. The method for simulating the diffusion of non-point source pollution along the irrigation route in coastal saline-alkali land irrigation areas based on drainage ditch network topology according to claim 2, characterized in that, In step S3, the first k The nitrogen absorption coefficient of a directed edge is given by the formula: ； In the formula, K N，k Indicates the first k The nitrogen absorption coefficient of a directed edge; K edge,k Indicates the first k The attenuation coefficient of a directed edge; V f,k Indicates the first k The denitrification rate per unit area of ​​a directed edge; W k Indicates the first k The average width of the water surface on each directed edge; Q k Indicates the first k The flow of a directed edge is less than or equal to Q edge,k ; Indicates the first k Localized comprehensive correction parameters for each directed edge; K N,k Indicates the first k The nitrogen absorption coefficient of a directed edge; No. k The phosphorus absorption coefficient for each directed edge is given by the following formula: ； In the formula, K p,k Indicates the first k The phosphorus absorption coefficient of a directed edge; K edge,k Indicates the first k The attenuation coefficient of a directed edge; α p,k and β p,k Indicates the first k The localization parameters of a directed edge; T r,k Indicates the first k Hydraulic residence time of a directed edge; No. k The salt transport coefficient of a directed edge is given by the formula: ； In the formula, K s,k Indicates the first k Salt transport coefficient of a directed edge; K edge,k Indicates the first k The attenuation coefficient of a directed edge; E k Indicates the first k The evaporation intensity of a directed edge of water surface; W k Indicates the first k The average width of the water surface on each directed edge; ρ Indicates the density of water; Q k Indicates the first k The flow of a directed edge; This represents the soil salt release coefficient on the slope.

4. The method for simulating the diffusion of non-point source pollution along the irrigation route in coastal saline-alkali land irrigation areas based on drainage ditch network topology according to claim 1, characterized in that, Step S5 includes the following steps: S5.1: Calculate the first x The initial nitrogen load for each farmland unit is calculated using the following formula: ； In the formula, L N,x Indicates the first x Initial nitrogen load of each farmland unit; A x Indicates the first x The area of ​​each farmland unit; F N,x Indicates the first x Nitrogen fertilizer application rate per unit area of ​​each farmland unit; μ N Represents the empirical coefficient; C season Indicates crop growth period or hydrological scenario correction factor; S5.2: Calculate the first x The initial phosphorus load for each farmland unit is calculated using the following formula: ； In the formula, A x Indicates the first x The area of ​​each farmland unit; F P,x Indicates the first x Phosphate fertilizer application rate per unit area of ​​each farmland unit; μ P Indicates the phosphorus loss rate; C clay This represents the correction factor for soil clay content; L P,x Indicates the first x Initial phosphorus load of each farmland unit; S5.3: Calculate the first x The initial salinity load of each farmland unit is calculated using the following formula: ； In the formula ，Q drain,x Indicates the first x Drainage of each farmland unit; C salt,x Indicates the first x Initial salinity concentration of each farmland unit; Indicates the salt concentration coefficient; S 0,x Indicates the first x Initial salinity load of each farmland unit; A x Indicates the first x The area of ​​each farmland unit; P Indicates rainfall amount; I rr,x Indicates the first x Irrigation volume for each farmland unit; Indicates the first x The runoff and drainage coefficients of each farmland unit.

5. The method for simulating the diffusion of non-point source pollution along the irrigation route in coastal saline-alkali land irrigation areas based on drainage ditch network topology according to claim 4, characterized in that, Step S6 includes the following steps: S6.1: Calculate the first... x The nitrogen inflow load of each farmland unit into the river is calculated using the following formula: ； In the formula, L N,out,x Indicates the first x Nitrogen load into the river in each farmland unit; L N,x Indicates the first x Initial nitrogen load of each farmland unit; D k This represents the length of the k-th directed edge; K N，k Indicates the first k The nitrogen absorption coefficient of a directed edge; Z x Indicates the first x The number of directed edges covered by the path of each farmland unit; S6.2: Calculate the first... x The phosphorus inflow load of each farmland unit into the river is calculated using the following formula: ； In the formula, L P,out,x Indicates the first x Phosphorus load into the river in each farmland unit; L P,x Indicates the first x Initial phosphorus load of each farmland unit; D k This represents the length of the k-th directed edge; K p,k Indicates the first k The phosphorus absorption coefficient of a directed edge; Z x Indicates the first x The number of directed edges covered by the path of each farmland unit; S6.3: Calculate the first x The formula for the salinity load into the river for each farmland unit is: ； In the formula, S out,x Indicates the first x The salt load of each farmland unit entering the river; S 0,x Indicates the first x Initial salinity load of each farmland unit; D k This represents the length of the k-th directed edge; K s,k Indicates the first k Salt transport coefficient of a directed edge; Z x Indicates the first x The number of directed edges covered by the path of a farmland unit.

6. The method for simulating the diffusion of non-point source pollution along the irrigation route in coastal saline-alkali land irrigation areas based on drainage ditch network topology according to claim 4, characterized in that, Step S7 includes the following steps: S7.1: Nitrogen inflow load into the river for all farmland units in the target irrigation district, calculated using the following formula: ； In the formula, L N,out,x Indicates the first x Nitrogen load into the river in each farmland unit; Ng Indicates the number of farmland units in the target irrigation district; TL N,out,x This represents the nitrogen inflow load into the river for all farmland units in the target irrigation district; S7.2: Phosphorus inflow load into the river for all farmland units in the target irrigation area, calculated using the following formula: ； In the formula, L P,out,x Indicates the first x Phosphorus load into the river in each farmland unit; Ng Indicates the number of farmland units in the target irrigation district; TL P,out,x This represents the phosphorus inflow load into the river for all farmland units in the target irrigation district; S7.3: Salt inflow load into the river for all farmland units in the target irrigation district, calculated using the following formula: ； In the formula, S out,x Indicates the first x The salt load of each farmland unit entering the river; TS out This represents the salt load into the river for all farmland units in the target irrigation area; Ng This indicates the number of farmland units in the target irrigation district.

7. The method for simulating the diffusion of non-point source pollution along the irrigation route in coastal saline-alkali land irrigation areas based on drainage ditch network topology according to claim 1, characterized in that, Step S1 includes the following steps: S1.1: Obtain basic data for the target irrigation area; S1.2: The DEM local correction preprocessing algorithm is used to rasterize the drainage ditch vector data to generate drainage ditch raster data with the same spatial resolution as the original DEM data. S1.3: Based on the rasterized drainage ditch data, obtain the "downcut depth" value of each ditch cell, and perform cell-by-cell correction on the DEM data of the target irrigation area. The formula is as follows: ； In the formula, H corrected ( i,j ) indicates position ( i,j The corrected elevation value at point ) H original ( i,j ) indicates position ( i,j The elevation value before correction at point ) D ( i,j ) indicates position ( i,j The incision depth of the ditch pixel at point ).

Images