Yellow River irrigation area non-point source pollution load space-time analysis method

Abstract

The invention provides a Yellow River water irrigation area non-point source pollution load space-time analysis method, and relates to the technical field of Yellow River water irrigation area non-point source pollution monitoring and analysis, and the method comprises the following steps: S1, collecting and preprocessing multi-source basic data; according to the Yellow River water irrigation area non-point source pollution load space-time analysis method, total nitrogen and ammonia nitrogen are accurately recognized as core pollutants through three-stage fuzzy comprehensive evaluation, the recognition accuracy is improved to 90% or above, and the problem that recognition is fuzzy in a traditional method is solved; by combining nitrogen and oxygen isotopes and a flora active chip, the manure sewage is determined as a main source of nitrogen for the first time, the traceability error is controlled within 10%, and the pollution source analysis blank is filled up; a core pollution area is defined through a gravity center migration model, the pollution load space-time distribution visualization degree is improved by 60%, and area differentiation treatment is achieved; a coupling migration model is constructed, the mechanism of reverse migration of pollutants caused by irrigation return water is clarified for the first time, and the interpretation precision reaches 85%.

Description

Technical Field

[0001] This invention relates to the field of non-point source pollution monitoring and analysis technology in the Yellow River irrigation area, and particularly to a spatiotemporal analysis method for non-point source pollution load in the Yellow River irrigation area. Background Technology

[0002] The Yellow River irrigation area is an important grain production base in northern my country, but it has long been affected by agricultural irrigation and fertilization activities. In the southern irrigation area, the risk of groundwater and drainage water pollution is significantly higher than that of surface water, with pollutants such as total nitrogen, ammonia nitrogen, and fluoride exceeding standards severely. Existing analytical methods rely on single monitoring data (such as water quality sampling) and do not incorporate technologies such as isotope tracing, microbial activity monitoring, and remote sensing dynamic monitoring, making it impossible to accurately identify specific pollution sources such as manure and sewage. Analysis of the reverse migration pattern of pollutants along the Yellow River runoff caused by irrigation backflow is lacking, making it difficult to explain the phenomenon of "higher pollutant concentrations downstream than upstream" after irrigation. The application of tools such as the Invest model and Nemerow index does not consider the unique "topographical runoff-irrigation-fertilization superposition" driving effect of the Yellow River irrigation area, resulting in low accuracy in simulating pollution loss (error exceeding 30%). Furthermore, the lack of correlation analysis between pollutant cluster types and loss driving factors fails to provide a basis for classification and control. For example, the shortcomings of current technologies include: pollution source analysis relies solely on conventional water quality indicators without incorporating specific markers such as nitrogen and oxygen isotopes and functional microbial activity, making it impossible to distinguish between different nitrogen sources such as manure, chemical fertilizers, and soil organic nitrogen; the analysis of the spatiotemporal distribution of pollution loads does not incorporate centroid migration models and UAV remote sensing data, resulting in vague definitions of core pollution areas (such as the eastern part of the irrigation area); and the simulation of pollutant migration does not construct a coupled model of "irrigation return water - soil infiltration - surface runoff," making it impossible to explain the causes of reverse migration. Summary of the Invention

[0003] Based on the technical problems existing in the background art, the present invention proposes a spatiotemporal analysis method for non-point source pollution load in the Yellow River irrigation area.

[0004] This invention proposes a spatiotemporal analysis method for non-point source pollution load in the Yellow River irrigation area, comprising the following steps: S1. Multi-source basic data acquisition and preprocessing: Collect basic data on agricultural non-point source pollution, water quality monitoring data of the Yellow River main stream and irrigation area, obtain dynamic data of surface cover in irrigation area using UAV multispectral technology, and collect soil pore water in layers and detect pollutant concentration using an in-situ soil pore water acquisition-detection integrated device. At the same time, collect meteorological data and DEM topographic data. S2. Identification of major pollutants and determination of priority control factors: The fuzzy comprehensive evaluation method is adopted to construct a three-level evaluation index system of "pollutant concentration-ecological risk-human health risk". Based on the water quality monitoring data obtained in S1, the membership degree and weight of each pollutant are calculated to identify the major pollutants and priority control factors of the irrigation area water body. S3. Source analysis and migration impact mechanism analysis of core pollutants: Taking autumn irrigation in the irrigation area as the research object, water samples were collected. The nitrogen-oxygen isotope ratio analysis method was combined with a real-time monitoring chip for the activity of denitrifying functional bacteria to trace the main sources of nitrogen pollution. The mechanism of enhanced activity of denitrifying functional bacteria after autumn irrigation on the aggravation of groundwater denitrification was analyzed. The driving mechanism of irrigation activities on nitrogen pollution migration was clarified through indoor simulation experiments. S4. Analysis of the spatiotemporal distribution characteristics of pollution load and definition of the core area: The pollution index of different water bodies in the irrigation area is calculated using the Nemerow comprehensive pollution index method to determine the pollution level and risk differences. Based on the pollution load calculation model and the centroid migration model, the centroid coordinates of the pollution load in different seasons of each sub-region are calculated to define the core pollution area and establish a correlation model between irrigation amount, fertilizer amount, rainfall amount and pollution load. S5. Pollutant migration patterns and causes analysis: The concentration of pollutants at different sections along the Yellow River runoff direction in the irrigation area was detected, and the changing trend of pollutants along the runoff direction was analyzed. Combined with the remote sensing data of S1 and the soil pore water detection data, a coupled migration model of "irrigation return water-soil seepage-surface runoff" was constructed to infer the main causes of the reverse migration of pollutants. S6. Pollution load clustering classification: The R-type clustering analysis method is adopted, using the concentration data of each pollutant obtained in S1-S5 as variables, setting a clustering threshold based on variance contribution rate > 85%, and calculating the similarity between variables through Euclidean distance to classify the pollution load clustering type. S7. Spatiotemporal Differentiation Quantitative Analysis of Pollutant Loss: Based on the Invest model, the DEM data, meteorological data, soil data, and crop planting data from S1 are used as inputs. A "topographic water catchment-irrigation and fertilization superimposed driving module" is added to simulate the loss of total nitrogen and total phosphorus in the southern irrigation area, analyze the temporal and spatial differentiation of the loss, and distinguish nitrogen and phosphorus migration paths through soil column experiments. S8. Analysis of Chemical Component Sources and Evaluation of Water Quality Suitability: A self-organizing neural network combined with positive definite matrix factor analysis was used to construct a "component-source" mapping model to analyze the main sources and contribution ratios of chemical components in canal water and drainage. An entropy weight water quality evaluation method was used to construct a "drinking suitability" evaluation index system, calculate the comprehensive water quality index, and determine the proportion of samples that are not suitable for drinking.

[0005] Preferably, the basic data on agricultural non-point source pollution in S1 includes crop planting structure, growth period of crops such as sunflower, fertilizer type and application amount, irrigation system and irrigation water volume, and irrigation area water bodies including groundwater, drainage water bodies, surface water, canal water and drainage. The detection of water quality monitoring data is carried out by high performance liquid chromatography and ion chromatography, with a monitoring frequency of once a month, and three parallel samples are collected each time. The relative deviation of parallel samples is controlled within 5%.

[0006] Preferably, the UAV multispectral technology in S1 uses a hexacopter UAV equipped with a 5-band hyperspectral imager, flying at an altitude of 100m with a spatial resolution of 0.5m, flying twice a month, with the flight path covering the entire irrigation area. The acquired remote sensing data is atmospherically and geometrically corrected using ENVI software, and then the NDVI index and soil moisture content are obtained through inversion. The soil pore water in-situ acquisition-detection integrated device includes a 120cm long stainless steel sampling tube, a 0.45μm mixed fiber filter membrane wrapped in the tube wall, and an in-situ ion sensor embedded in the tube. The stratified acquisition depth is 0-20cm, 20-60cm, and 60-100cm. Each sampling point is sampled 3 times, and the average value is taken as the detection result.

[0007] Preferably, the specific steps of the fuzzy comprehensive evaluation method in S2 are as follows: S21. Determine the evaluation factor set U = {total nitrogen (u1), ammonia nitrogen (u2), nitrite nitrogen (u3), nitrate nitrogen (u4), total phosphorus (u5), TDS (u6), chemical oxygen demand (u7), perfluorinated pollutants (u8), fluorides (u9)}; S22. Determine the evaluation level set V = {Clean (v1), Slightly Polluted (v2), Moderately Polluted (v3), Heavily Polluted (v4)}, where the clean level corresponds to the concentration of each pollutant being lower than the Class III standard of GB / T14848-2017 and the Class III standard of GB3838-2022, the lightly polluted level corresponds to the Class III-IV standard, the moderately polluted level corresponds to the Class IV-V standard, and the heavyly polluted level corresponds to the Class V standard. S23. Using the analytic hierarchy process (AHP), five experts in the field of environmental engineering were invited to construct a judgment matrix and calculate the weight vector W={w1,w2,...,w9} for each evaluation factor. The rationality of the weights was verified through a consistency test, where the weight of total nitrogen was w1=0.25, the weight of ammonia nitrogen was w2=0.20, and the weight of fluoride was w9=0.15. S24. Construct a fuzzy relation matrix R, where rij represents the membership degree of the i-th evaluation factor to the j-th evaluation level. The membership degree is calculated using a descending half-trapezoidal distribution function. When the factor concentration ≤ the cleanliness standard, the membership degree r... i1 =1、r i2 =r i3 =r i4 =0, when the concentration is between the clean and slightly polluted standards, r i1 =(Mild Standard - Concentration) / (Mild Standard - Cleanliness Standard), r i2 =1-r i1 r i3 =r i4 =0; S25. The weighted average method is used to calculate the fuzzy comprehensive evaluation result B=W×R, and the main pollutants and priority control factors are determined according to the principle of maximum membership.

[0008] Preferably, the specific steps of the nitrogen-oxygen isotope ratio analysis method in S3 are as follows: S31. Collect 500 mL of water sample, filter it through a 0.45 μm filter membrane, add 10 mL of concentrated sulfuric acid to acidify it to pH < 2, and use distillation to separate nitrogen from the sample and convert nitrogen into silver nitrate precipitate. S32. Determination of δ-carbon dioxide in silver nitrate precipitate using stable isotope ratio mass spectrometry. 15 N and δ 18 The ratio of O, where δ 15 The test accuracy for N is ±0.2‰, and δ 18 The testing accuracy for O is ±0.3‰, with reference to international standard materials AIR and VSMOW. S33. Establish a database of isotope ratios for known pollution sources: δ¹⁸O⁻ of manure and sewage. 15 N=7-15‰, δ 18 O=5-12‰, fertilizer δ 15 N=-5-5‰、δ 18 O = -10⁻⁵‰, atmospheric deposition δ 15 N=-10-0‰、δ 18 O=10-20‰, soil organic nitrogen δ 15 N=2-8‰, δ 18 O=8-15‰, the sample isotope ratios are compared with the database to determine the main sources of nitrogen pollution; The real-time monitoring chip for denitrifying microbial activity uses a silicon-based nanoelectrode array. The activity of nitratereductase enzyme is used as the detection index, with a detection limit of 0.01 U / L. The monitoring frequency is 1 day before autumn irrigation, 1 day after autumn irrigation, 3 days after autumn irrigation, and 7 days after autumn irrigation. Three parallel samples are collected for each monitoring. The activity data are analyzed for trends using Origin software.

[0009] Preferably, the formula for calculating the Nemerow Integrated Pollution Index in S4 is: PI = √[(PI avg 2 +PI max 2 ) / 2], where PI is the Nemerow Composite Pollution Index, PI avg PI is the arithmetic mean of the pollution indices of each pollutant. max The pollution index P represents the maximum value among all pollution indices. Ii =C i / S i C i S represents the measured concentration of the pollutant. iThe standard concentrations for pollutant evaluation are as follows: groundwater evaluation adopts Class III standard of the "Groundwater Quality Standard" (GB / T14848-2017), and surface water evaluation adopts Class III standard of the "Surface Water Environmental Quality Standard" (GB3838-2022); the formula for calculating the centroid coordinates of the centroid migration model is: X=∑(Q i ×X i ) / ∑Q i Y=∑(Q i ×Y i ) / ∑Q i Where X and Y are the abscissa and ordinate of the pollution load centroid, respectively, and Qi is the pollution load of the i-th evaluation unit. i Y i These are the x-coordinate and y-coordinate of the geometric center of the i-th evaluation unit, respectively. The evaluation units are divided into 1km×1km grids, with a total of 500 evaluation units. The correlation model was analyzed using Pearson linear correlation analysis. The correlation coefficients R between irrigation amount, fertilizer amount, rainfall amount and pollution load amount were calculated using SPSS software. The correlation coefficients R between irrigation amount and pollution load amount were 0.85, R between fertilizer amount and rainfall amount were 0.82, and R between rainfall amount and pollution load amount were 0.74. All of them passed the significance test of P < 0.01.

[0010] Preferably, the construction steps of the coupled migration model of "irrigation return water-soil infiltration-surface runoff" in S5 are as follows: S51. Construct a soil seepage sub-model based on the HYDRUS-1D model, input the soil data collected in S1, and simulate the soil moisture seepage rate under different irrigation amounts. S52. Construct a surface runoff sub-model based on the SWAT model. Input the DEM data, land use data (70% cultivated land, 20% grassland, and 10% water area) and meteorological data collected in S1, and calculate the surface runoff flow using the SCS curve number method. S53. Construct an irrigation return water sub-model, input irrigation water volume, ditch slope, and soil infiltration rate, and calculate irrigation return water flow rate through the water balance equation Q return water = Q irrigation - Q infiltration - Q evaporation. S54. A coupling interface was developed using C++ to couple the output data of the soil infiltration sub-model, surface runoff sub-model, and irrigation return water sub-model. The migration coefficient of pollutants in infiltration water was set to 0.8, in surface runoff to 0.9, and in irrigation return water to 0.95. The migration paths and concentration changes of pollutants in the water-soil system were simulated. The model was validated using the measured concentration data collected in S1, and the validation error was controlled within 15%.

[0011] Preferably, the R-type clustering analysis method in S6 adopts the Ward clustering method, implemented using SPSS 26.0 software. The specific steps are as follows: S61. Standardize the concentration data of the nine pollutants obtained in S1-S5 (standardization formula: x'=(xx)). min ) / (x max -x min (Eliminating the influence of dimensions). S62. Calculate the Euclidean distance between standardized variables (distance formula d). ij =√[∑(x ik -x jk Construct the distance matrix; S63. Use Ward's method for clustering, which minimizes the increase in the sum of squares within each cluster when merging two clusters. Set the clustering threshold to a variance contribution rate > 85%. Calculate the variance contribution rate of each principal component through principal component analysis. Stop clustering when the cumulative variance contribution rate reaches 85%. S64. The silhouette coefficient is used to verify the validity of the clustering results. The formula for calculating the silhouette coefficient is s. i =(b i -a i ) / max(a i ,b i ), where a i Let b be the average distance between the i-th sample and other samples of the same class. i Let be the average distance between the i-th sample and the nearest outlier sample. The clustering result is considered valid when the average silhouette coefficient is greater than 0.7. The final four pollution load clustering types are as follows: Clusters of reducing nitrogen, total phosphorus-TDS-chemical oxygen demand, perfluorinated pollutants, and total nitrogen.

[0012] Preferably, the parameter optimization of the Invest model in S7 includes: The field water holding capacity and saturated hydraulic conductivity in the soil data were adjusted to the measured values ​​based on S1 in-situ detection. The growth period in the crop parameters was adjusted to the actual growth period of sunflowers in the Yellow River irrigation area. The crop coefficients were adjusted to: seedling stage 0.3, jointing stage 0.7, budding stage 1.1, grain filling stage 1.0, and maturity stage 0.6. The "Topographic Water Catchment-Irrigation and Fertilization Overlay Driving Module" is an independently developed module written in Python. The input parameters include topographic slope, irrigation amount, nitrogen application amount, and average monthly rainfall. The module calculates the pollution loss driving intensity under different parameter combinations through overlay analysis. Soil column experiments were conducted using transparent plexiglass columns filled with undisturbed soil. Three treatment groups were set up, with each group being replicated three times. Soil solutions at different depths were collected through sampling holes to detect nitrogen and phosphorus concentrations. The migration rates of nitrogen and phosphorus were calculated, and it was found that nitrogen mobility was significantly higher than that of phosphorus.

[0013] Preferably, the specific steps of the self-organizing neural network combined with the positive definite matrix factor analysis method in S8 are as follows: S81. A self-organizing neural network was used to reduce the dimensionality of the chemical composition data of canal water and drainage samples. The SOM network was set to a 4×4 output layer and 11 input layer neurons. The initial learning rate was 0.05 and the number of iterations was 1000. After dimensionality reduction, 6 feature vectors were obtained. S82. Input the 6 eigenvectors into the positive definite matrix factor analysis model, set the number of factors to 3, and determine the characteristic components of each source by rotating the factor loading matrix: the characteristic components of evaporite dissolution are TDS, Cl-, SO4²-, the characteristic components of silicate weathering are Ca²+, Mg²+, HCO3-, and the characteristic components of agricultural input are total nitrogen, total phosphorus, and fluoride. S83. Calculate the contribution percentage of each source: the contribution percentage of evaporite dissolution is 35-45%, the contribution percentage of silicate weathering is 20-30%, and the contribution percentage of agricultural input is 30-40%. The error in the contribution percentage is calculated using the bootstrap method, with the error controlled within 5%. The specific steps of the entropy weight water quality assessment method are as follows: S84. Standardize water quality indicators. The standardization formula for positive indicators is x'=(xx). min ) / (x max -x min The standardized formula for the negative index is x'=(x max -x) / (x max -x min ); S85. Calculate the entropy value e of each indicator. j =-k×∑(p ij ×lnp ij ), where k=1 / lnm (m is the number of samples, in this example m=100), p ij =x' ij / ∑x' ij (If x' ij =0, let p ij =1 / m); S86. Calculate the weight w of each indicator. j =(1-e j ) / ∑(1-e j ); S87, Calculate the comprehensive water quality index FI = ∑(w j ×Fij ), where F ij FI is the score of the j-th indicator for the i-th sample. When FI < 0.6, it is judged as suitable for drinking; when FI ≥ 0.6, it is judged as unsuitable for drinking.

[0014] The spatiotemporal analysis method for non-point source pollution load in the Yellow River irrigation area proposed in this invention has the following beneficial effects: Through a three-level fuzzy comprehensive evaluation, total nitrogen and ammonia nitrogen are accurately identified as core pollutants, with an identification accuracy rate exceeding 90%, solving the problem of fuzzy identification in traditional methods; combining nitrogen and oxygen isotopes with microbial activity chips, manure wastewater is identified as the main source of nitrogen for the first time, with the source tracing error controlled within 10%, filling a gap in pollution source analysis; the core pollution area is defined through a centroid migration model, improving the visualization of the spatiotemporal distribution of pollution load by 60%, enabling differentiated regional governance; and a coupled migration model is constructed, elucidating for the first time that irrigation backflow leads to the reverse migration of pollutants. The mechanism explains the accuracy of 85%; it classifies four types of pollution load clusters to enable the formulation of classified control strategies (such as focusing on denitrifying bacteria for reducing nitrogen control); it optimizes the Invest model, reducing the simulation error of nitrogen and phosphorus loss to below 15%, which is significantly better than the traditional model (error 30%+); the entropy weight method combined with SOM-PMF improves the accuracy of water quality drinking suitability determination to 88%, and the accuracy of chemical component source contribution analysis to 90%; it integrates eight types of technologies and six models to form a closed loop of "data-analysis-application", which can provide comprehensive technology and support the precise treatment of non-point source pollution in the Yellow River irrigation area, and can reduce the pollution load emission in the irrigation area by 20-30%. Attached Figure Description

[0015] Figure 1 This is a flowchart of a spatiotemporal analysis method for non-point source pollution load in the Yellow River irrigation area proposed in this invention. Detailed Implementation

[0016] Reference Figure 1 This invention proposes a spatiotemporal analysis method for non-point source pollution load in the Yellow River irrigation area, comprising the following steps: S1. Multi-source basic data acquisition and preprocessing 1.1 Basic Pollution Data Collection For the Yellow River irrigation area (taking the Yellow River South Irrigation Area as a typical study area, including the Jianshe Irrigation Area and the Dugui Hangjin Irrigation Area), the following basic data were collected: Basic data on agricultural non-point source pollution: Through a combination of field surveys and farmer questionnaires, the following data were obtained: crop planting structure (sunflower 65%, corn 25%, other crops 10%), growth stages of sunflower and other crops (emergence early to mid-April, jointing stage mid to late May, budding stage mid to late June, grain filling stage July to late August, maturity stage mid to late September), fertilizer types (urea: N content 46%; diammonium phosphate: N content 18%, P2O5 content 46%) and application rates (40% basal fertilizer, applied at the emergence stage; 60% topdressing, applied three times at the jointing, budding, and grain filling stages, with a total nitrogen application rate of 300-375 kg). hm -2 Irrigation system (spring irrigation in early April, autumn irrigation in early September, winter irrigation in early November, irrigation method is furrow irrigation, irrigation water volume 900-1350m³) 3 hm -2 ); Water quality monitoring data: Water quality indicators were collected from the Yellow River main stream and irrigation areas (groundwater: collected from monitoring wells drilled to a depth of 30m; drainage water: cross-sections of drainage canals at the end of the irrigation area; surface water: water bodies in tributaries surrounding the irrigation area; canal water: cross-sections of main and branch irrigation canals; drainage: drainage outlets of underground pipes in farmland), including pH value (detection range 0-14, accuracy ±0.01), dissolved oxygen (DO, unit mg / L, detection limit 0.01 mg / L), total nitrogen (TN, unit mg / L), and ammonia nitrogen (NH4). + -N, unit mg / L), nitrite nitrogen (NO2) - -N, unit mg / L), nitrate nitrogen (NO3) - -N (mg / L), total phosphorus (TP, mg / L), total dissolved solids (TDS, mg / L), chemical oxygen demand (COD, mg / L), perfluorinated pollutants (PFAS, μg / L), and fluorides (F) - (unit: mg / L) Detection instruments and methods: TN was detected using alkaline potassium persulfate digestion-ultraviolet spectrophotometry (instrument: Shimadzu UV-2600 ultraviolet spectrophotometer, detection limit 0.05 mg / L); NH4 + -N was detected using Nessler's reagent spectrophotometry (instrument: Hach DR6000 UV-Vis spectrophotometer, detection limit 0.025 mg / L); NO3 - -N, NO2 --N was analyzed using ion chromatography (instrument: Dionex ICS-2100 ion chromatograph, detection limit 0.001 mg / L); TP was analyzed using ammonium molybdate spectrophotometry (instrument: same as above, detection limit 0.01 mg / L); COD was analyzed using potassium dichromate method (instrument: Hach COD digester DRB200, detection limit 5 mg / L); PFAS was analyzed using ultra-high performance liquid chromatography-tandem mass spectrometry (instrument: Waters ACQUITYUPLC-XevoTQ-S, detection limit 0.01 μg / L); F - Ion-selective electrode method was used (instrument: Leici PXSJ-216 ion meter, detection limit 0.01 mg / L). Sample processing requirements: Sampling should be conducted once a month, with 3 parallel samples collected each time. The sampling volume is 500 mL (PTFE bottles for PFAS sampling, and brown glass bottles for other indicators). Samples should be stored at 4℃ after collection and tested within 24 hours. The relative deviation of parallel samples should be ≤5%. If it exceeds 5%, resampling should be carried out. Outliers should be removed using the Grubbs method (significance level α=0.05).

[0017] 1.2 Remote Sensing Dynamic Monitoring Drone configuration: A DJI Matrice 350RTK hexacopter drone equipped with a DJI Zenmuse P1 hyperspectral imager (5 bands: 450nm±5nm, 550nm±5nm, 650nm±5nm, 750nm±5nm, 850nm±5nm, spectral resolution 5nm), flight altitude 100m (ensuring ground sampling distance GSD=0.5m, spatial resolution 0.5m), flight speed 5m / s, forward overlap 80%, lateral overlap 60%; Flight and Data Processing: Flights were conducted twice a month (3 days before irrigation and 7 days after irrigation, capturing surface changes before and after irrigation). The flight paths covered the entire Yellow River South Irrigation Area (approximately 500 km²). Before each flight, ground control points (GCPs, 10 points deployed, RTK positioning, accuracy ±2 cm) were used for calibration. Remote sensing data preprocessing steps included: ① Atmospheric correction (using the FLAASH model, inputting atmospheric profile data and aerosol optical thickness); ② Geometric correction (based on GCPs, correction error ≤ 1 pixel); ③ Radiometric calibration (converting DN values ​​to reflectance); ④ Cropping (cropping according to irrigation area boundaries, removing non-study areas). Inversion parameters and methods: NDVI (Normalized Difference Vegetation Index) is calculated using the formula NDVI=(ρ 850 -ρ 650 ) / (ρ 850 +ρ 650NDVI is used to characterize crop growth (NDVI > 0.6 indicates vigorous growth, NDVI < 0.3 indicates slow growth); soil moisture content is determined using a single-band inversion model (based on 750nm band reflectance, formula: soil moisture content = -0.23 × ρ). 750 +0.35, R²=0.82, inversion error <5%), and inversion was performed using ENVI5.6 software.

[0018] 1.3 In-situ detection of soil pore water Device Structure: Utilizing a self-developed "In-situ Soil Pore Water Acquisition-Detection Integrated Device," the device comprises: ① Sampling Unit: A 120cm long, 5cm outer diameter 304 stainless steel sampling tube. Sampling holes (2mm diameter, 12 holes per group, spirally distributed) are opened every 20cm along the tube wall. Each sampling hole is wrapped with a 0.45μm mixed fiber filter membrane (material: cellulose acetate-cellulose nitrate, acid and alkali resistant, 0.45μm particle size cutoff, preventing soil particles from entering). The filter membrane is fixed by a double-layer stainless steel mesh (0.1mm pore size) wrapped around it, laser-welded to the outer wall of the sampling tube to prevent membrane detachment (existing devices often use glue, which is easily corroded by soil); ② Detection Unit: Three in-situ ion sensors are embedded inside the tube (corresponding to soil layers of 0-20cm, 20-60cm, and 60-100cm respectively; sensor model: Sentek GS3; detection parameter: NO3). - -N, NH4 + -N, response time <30s, detection range 0.1-100mg / L, accuracy ±0.05mg / L), sensor calibration procedure: use 10mg / L, 50mg / L, and 100mg / L NO3 - -N, NH4 + -N standard solution, perform two-point calibration of Sentek GS3 sensor under constant temperature of 25℃, calibration error must be ≤2%, field calibration once every 30 days; ③ Control unit: built-in data acquisition unit (sampling interval 1h, storage capacity 16GB) and lithium battery (30-day battery life), supports 4G remote data transmission, data transmission verification: control unit supports 4G Cat.1 transmission, data sampling interval can be adjusted by remote command (1h-24h), automatically generates data integrity report daily (alarm triggered when missing rate >5%). Sampling point layout and operation: In the Yellow River South Irrigation Area, 30 sampling points were laid out in a gradient of "upstream-midstream-downstream" (upstream: 8 points with a slope of 5-8°; midstream: 14 points with a slope of 2-5°; downstream: 8 points with a slope of 0-2°), with a spacing of ≥5km between sampling points, avoiding interference areas such as ditches and roads; the sampling tubes were vertically inserted into the soil to a depth of 100cm using a hydraulic drilling rig, with the filter membrane in close contact with the soil, and testing began after standing for 24 hours. The testing frequency was synchronized with the water quality monitoring (once a month), with each test lasting 24 hours, and the average value was taken as the pore water pollutant concentration of that soil layer.

[0019] 1.4 Auxiliary Data Acquisition Meteorological data: Daily meteorological data for the past five years (e.g., 2020-2024) were collected from three national meteorological stations surrounding the Yellow River South Irrigation Area (e.g., Bayannur, Wuhai, and Ordos). This included daily precipitation (mm, accuracy 0.1mm), daily average temperature (°C, accuracy 0.1°C), daily average wind speed (m / s, accuracy 0.1m / s), and daily average relative humidity (%, accuracy 1%). The data were sourced from the China Meteorological Data Network, and missing data were supplemented using linear interpolation. Topographic data: 30m resolution DEM (Digital Elevation Model) data were acquired from the Geospatial Data Cloud and processed using ArcGIS 10.8 software to extract topographic parameters such as slope (accuracy 0.1°), aspect, and catchment area, which were then used as topographic factor inputs for subsequent models.

[0020] S2. Identification of Major Pollutants and Determination of Priority Control Factors The three-level fuzzy comprehensive evaluation method of "pollutant concentration-ecological risk-human health risk" is adopted, and the specific steps are as follows: 2.1 Construction of Evaluation Factor Set The evaluation factor set U = {total nitrogen (u1), ammonia nitrogen (u2), nitrite nitrogen (u3), nitrate nitrogen (u4), total phosphorus (u5), TDS (u6), chemical oxygen demand (u7), perfluorinated pollutants (u8), fluorides (u9)} is set to cover all characteristic pollutants monitored in the irrigation area, ensuring that no key factors are missed.

[0021] 2.2 Setting of Evaluation Level Set The evaluation level set V = {Clean (v1), Slightly Polluted (v2), Moderately Polluted (v3), Heavyly Polluted (v4)} is defined, and the level standards strictly refer to the current national standards. Groundwater: Based on the "Groundwater Quality Standard" (GB / T14848-2017), Class III (clean v1), Class IV (slightly polluted v2), Class V (moderately polluted v3), and above Class V (heavily polluted v4), the specific values ​​are as follows: u1≤10mg / L (v1), 10<u1≤20mg / L (v2), 20<u1≤30mg / L (v3), u1>30mg / L (v4); u2≤0.5mg / L (v1), 0.5<u2≤1.0mg / L (v2), 1.0<u2≤2.0mg / L (v3), u2>2.0mg / L (v4); u9≤1.0mg / L (v1), 1.0<u9≤2.0mg / L (v2), 2.0<u9≤4.0mg / L (v3), u9>4.0mg / L (v4). Surface water: Based on the "Surface Water Environmental Quality Standard" (GB3838-2022), Class III (clean v1), Class IV (lightly polluted v2), Class V (moderately polluted v3), and above Class V (heavily polluted v4), the specific values ​​are: u1≤1.0mg / L (v1), 1.0<u1≤2.0mg / L (v2), 2.0<u1≤5.0mg / L (v3), u1>5.0 mg / L (v4); u2≤0.5mg / L (v1), 0.5<u2≤1.0mg / L (v2), 1.0<u2≤2.0mg / L (v3), u2>2.0mg / L (v4); u5≤0.2mg / L (v1), 0.2<u5≤0.4mg / L (v2), 0.4<u5≤0.6mg / L (v3), u5>0.6mg / L (v4); Other factors: TDS (u6) refers to GB / T14848-2017: ≤1000mg / L (v1), 1000<u6≤2000mg / L (v2), 2000<u6≤3000mg / L (v3), u6>3000mg / L (v4); COD (u7) refers to GB3838-2022: ≤20mg / L (v1), 20<u7≤30mg / L (v2), 30<u7≤40mg / L (v3), u7>40mg / L (v4); PFAS (u8) refers to the "Standards for Drinking Water Quality" (GB5749-2022): ≤0.001mg / L (v1), 0.001<u8≤0.005mg / L (v2), 0.005<u8≤0.01mg / L (v3), u8>0.01mg / L (v4).

[0022] 2.3 Weight Calculation The weights of each factor are calculated using the Analytic Hierarchy Process (AHP). The steps are as follows: The hierarchical structure is constructed as follows: target layer (evaluation of water pollution in the Yellow River irrigation area), criterion layer (pollutant concentration, ecological risk, and human health risk, with weights of 0.4, 0.3, and 0.3, respectively), and indicator layer (i.e., evaluation factor set U). Invite experts to score: Invite 5 experts with associate senior or higher professional titles in the field of environmental engineering to construct a judgment matrix using a 1-9 scale (1 = equally important, 3 = slightly important, 5 = significantly important, 7 = strongly important, 9 = extremely important, 2, 4, 6, and 8 are intermediate values). For example, the judgment matrix for the criterion layer: |Criterion Level|Pollutant Concentration|Ecological Risk|Human Health Risk| |----------|------------|----------|--------------| |Pollutant concentration|1|3 / 1|4 / 1| |Ecological Risk|1 / 3|1|2 / 1| |Human Health Risks|1 / 4|1 / 2|1| Consistency check: Calculate the largest eigenvalue λ of the judgment matrix. max Consistency index CI=(λ) max -n) / (n-1) (n is the matrix order), random consistency index RI (RI=0.58 when n=3), consistency ratio CR=CI / RI, CR<0.1 is required, if not satisfied, the judgment matrix is ​​adjusted; Weight determination: The weights of each factor were calculated using the eigenvector method. The final weight vector is W={w1=0.25,w2=0.20,w3=0.08,w4=0.09,w5=0.12,w6=0.06,w7=0.07,w8=0.05,w9=0.15}, where total nitrogen (w1) and ammonia nitrogen (w2) account for 45% of the weights, which meets the requirement of "not less than 40%".

[0023] 2.4 Membership Degree and Comprehensive Evaluation The membership degree of each factor to different pollution levels is calculated using a descending half-trapezoidal distribution function. The specific formula is as follows: The membership degree r of cleanliness (v1) ij1 : When C i ≤S i1 When (v1 standard value), r ij1 =1; When S i1 <C i ≤S i2 At (v2 standard value), r ij1=(S i2 -C i ) / (S i2 -S i1 ); When C i >S i2 At that time, r ij1 =0; Membership r of mild pollution (v2) ij2 : When C i ≤S i1 At that time, r ij2 =0; When S i1 <C i ≤S i2 At that time, r ij2 =(C i -S i1 ) / (S i2 -S i1 ); When S i2 <C i ≤S i3 (v3 standard value) r ij2 =(S i3 -C i ) / (S i3 -S i2 ); When C i >S i3 At that time, r ij2 =0; When S i3 <C i ≤S i4 At that time, r ij4 =(C i -S i3 ) / (S i4 -S i3 ); When C i >S i4 At that time, r ij4 =1; Membership r of heavily polluted (v4) ij4 : When C i ≤S i3 (v3 standard value) r ij4 =0; When S i3 <C i ≤ Si4 (v4 standard value) r ij4 =(Ci -S i3 ) / (S i4 -S i3 ); When C i >S i4 At that time, r ij4 =1; Example: If the measured concentration C of ammonia nitrogen (u2) in a groundwater sample is... i =2.5mg / L, corresponding to S i3 =2.0mg / L, S i4 =5.0mg / L, then r ij4 =(2.5-2.0) / (5.0-2.0)=0.167.

[0024] Among them, C i S represents the measured concentration of the i-th factor (unit: mg / L or μg / L; PFAS units must be uniformly expressed as mg / L). i1 S i2 S i3 S i4 These are the standard values ​​for levels v1, v2, v3, and v4 corresponding to the i-th factor, respectively; Construct a fuzzy relation matrix R: R is a 9×4 matrix, where each row corresponds to one evaluation factor, each column corresponds to one pollution level, and the elements are the membership degrees r. ij ; Comprehensive evaluation calculation: The weighted average method is used to calculate the comprehensive evaluation vector B = W × R (matrix multiplication), that is, B = [b1, b2, b3, b4], where b j =∑(w i ×r ij (i=1-9, j=1-4); Result determination: Based on the "maximum membership principle", select b. j The level corresponding to the maximum value is taken as the pollution level of the water body, and the top 3 factors with the highest weight (total nitrogen u1, ammonia nitrogen u2, fluoride u9) are extracted as priority control factors.

[0025] S3. Source apportionment and migration impact mechanism analysis of core pollutants 3.1 Nitrogen Source Traceability Autumn irrigation in the Yellow River South Irrigation Area (early September, irrigation volume 1125m³) 3 hm -2 Taking ( ) as the research object, the specific steps are as follows: Sample collection: Assuming 10 sampling points are set up in each of the Jianshe Irrigation Area and the Dugui Hangjin Irrigation Area, 20 samples (500 mL each) of groundwater (30m depth of monitoring well), drain water (middle layer of cross section), and canal water (middle layer of irrigation branch canal) are collected. The samples are collected slowly using a sampler to avoid disturbance. PFAS test samples are collected in polytetrafluoroethylene bottles, and other samples are collected in brown glass bottles. The samples are transported under refrigeration at 4℃. Isotope pretreatment: Take 200 mL of water sample and filter it through a 0.45 μm mixed fiber filter membrane to remove suspended particulate matter; Add 10 mL of concentrated sulfuric acid (analytical grade, 98% concentration) to acidify to pH < 2 to inhibit microbial activity; A semi-micro Kjeldahl distillation apparatus (instrument: Shanghai Fiber Inspection KDN-103F) was used, with excess magnesium oxide (superior grade) added (magnesium oxide dosage: 0.5 g of superior grade magnesium oxide (100 mesh particle size, pre-dried at 105℃ for 2 h) per 200 mL water sample, to ensure that the pH is stable at 9-10 during the distillation process), the distillation temperature was 105℃, the distillation time was 30 min, and the distillate was collected in a receiving flask containing 20 mL of 2% boric acid solution; Add 0.1 mol / L silver nitrate solution (analytical grade) to the distillate (end point of silver nitrate addition: add 0.1 mol / L silver nitrate solution until a white precipitate appears in the solution and does not disappear within 30 seconds, record the amount used (usually 8-12 mL), to ensure that nitrogen is completely converted to AgNO3), until no more precipitate is formed, let stand for 24 h, filter and collect the silver nitrate precipitate (AgNO3), wash 3 times with ultrapure water, and dry at 80℃ to constant weight (purity > 99.5%). Isotope determination: δ¹² Isotope concentrations in AgNO₃ were determined using a stable isotope ratio mass spectrometer (IRMS, Thermo Delta VAdvantage). 15 N and δ 18 O, Measurement conditions: Ion source temperature 1300℃, ionization voltage 10kV, resolution 500, measurement accuracy: δ 15 N±0.2‰, δ 18 O±0.3‰, using the international standard material AIR (atmospheric nitrogen, δ) 15 N=0‰) is δ 15 N reference, VSMOW (Vienna Standard Mean Seawater, δ) 18 O=0‰) is δ 18 O reference, isotope ratio calculation formula: δ(‰) = (Rsample / Rstandard - 1) × 1000, where R is the ratio of heavy isotope to light isotope (e.g., ...). 15 N / 14 N、 18 O / 16 O); Microbial activity monitoring: The "Real-time Monitoring Chip for Denitrifying Microbial Activity" is used. Chip structure: ① Substrate: Silicon wafer (2cm × 2cm); ② Nanoelectrode array: 50 × 50 gold nanoelectrodes (50nm diameter, 10μm spacing), modified with nitrate reductase-specific antibody; ③ Signal processing unit: Built-in microcurrent amplifier (detection range 10... -12 -10 -6 A) With wireless transmission module; During monitoring, take 10 mL of water sample (or soil pore water), add it to the chip reaction cell, react at 37℃ for 30 min, capture nitrate reductase secreted by bacteria through electrodes, and detect enzyme activity (unit U / L, 1U is defined as the amount of enzyme that catalyzes the reduction of 1 μmol of nitrate per minute at 37℃), detection limit 0.01U / L; Source identification: Establish an isotope database of known nitrogen sources: ① Manure wastewater (originating from livestock farms and domestic sewage in irrigation areas, δ¹²⁺) 15 N=7-15‰, δ 18 O=5-12‰, nitratereductase activity 0.5-1.2U / L); ② Fertilizers (urea and diammonium phosphate, δ 15 N=-5-5‰、δ 18 O=-10⁻⁵‰, nitratereductase activity <0.1U / L); ③ Soil organic nitrogen (topsoil in irrigated areas, δ¹⁰ 15 N=2-8‰, δ 18 O=8-15‰, nitratereductase activity 0.2-0.5U / L); ④ Atmospheric deposition (local rainfall, δ 15 N=-10-0‰、δ 18 O=10-20‰, nitratereductase activity <0.05U / L); δ of the sample 15 N, δ 18 By comparing O with the database and combining it with enzyme activity values, it was determined that the main source of nitrogen was manure wastewater (sample δ). 15 N=8-12‰, δ 18 O=6-10‰, enzyme activity 0.6-1.0U / L, matching degree >90%.

[0026] 3.2 Analysis of the impact on microbial community activity Monitoring Plan: Five typical sampling points (midstream agricultural intensive areas) were selected in the irrigation area. Groundwater samples were collected one day before autumn irrigation (T0), one day after autumn irrigation (T1), three days after autumn irrigation (T2), and seven days after autumn irrigation (T3). The nitrate reductase activity of Pseudomonas and Paracoccus was detected using the aforementioned "real-time monitoring chip for denitrification functional bacterial community activity". At the same time, high-throughput sequencing (instrument: Illumina MiSeq) was used to analyze the abundance of bacterial communities (primers: 338F5'-ACTCCTACGGGAGGCAGCAG-3', 806R5'-GGACTACHVGGGTWTCTAAT-3'). Monitoring results: At T0, the bacterial activity was 0.4-0.6 U / L, and the abundance was 1.2 × 10⁻⁶. 6 -1.5×10 6 Copies / g soil; at T1, activity increased to 0.6-0.8 U / L (an increase of 30-50%), and abundance increased to 1.8 × 10⁻⁶. 6 -2.2×10 6 Copies / g soil; activity remained at 0.7-0.9 U / L during T2-T3, and abundance remained stable at 2.0 × 10⁻⁶. 6 -2.3×10 6 copies / g soil; Indoor verification experiment: Take soil from the topsoil layer of the irrigation area (0-20cm, bulk density 1.3g / cm³), pass it through a 2mm sieve, put it into a 500mL serum bottle (soil volume 300g), add groundwater (simulating field moisture, moisture content 25%), adjust pH=7.0, purge nitrogen for 30min (anaerobic environment), and place it in a 25℃ constant temperature incubator. Three treatment groups were set up: ① control group (inoculated without aseptic colonies); ② low activity group (inoculated with Pseudomonas, activity 0.5 U / L); ③ high activity group (inoculated with Pseudomonas + Paracoccus, activity 0.8 U / L), with each group repeated 3 times. NO3 was detected after incubation for 0, 1, 3, 5, and 7 days. - -N concentration (ion chromatography), calculate denitrification rate (unit: mg / (kg)) d)); Results: The denitrification rate in the control group was 0.8-1.0 mg / (kg) d); Low activity group 1.0-1.2 mg / (kg) d); High-activity group: 1.2-1.35 mg / (kg) d) indicates that when the activity of the microbial community increases by 30-50%, the denitrification rate increases by 25-35%, exacerbating NO3. - -N is converted into gaseous nitrogen such as N2 and N2O, reducing NO3 in groundwater. - -N concentration, but increases the risk of greenhouse gas emissions.

[0027] 3.3 Verification of Irrigation Driving Mechanism Indoor simulation experiment design: An organic glass soil column (20cm in diameter, 100cm in height, filled with irrigated soil in layers: 0-20cm topsoil layer, 20-60cm plow pan layer, and 60-100cm parent material layer, with bulk densities of 1.3, 1.4, and 1.5g / cm³, respectively, hereinafter referred to as soil column) was used. A drainage outlet was set at the bottom of the soil column and connected to a water collection bottle. Treatment settings: 3 sets of irrigation volumes (W1=900m³) 3 hm -2 W2=1125m 3 hm -2 W3=1350m 3 hm -2 Converted to soil column irrigation volume: W1=1.8L, W2=2.25L, W3=2.7L, nitrogen application rate for 2 groups (N1=300kg) hm -2 N2 = 375 kg hm -2 The nitrogen application rates for the soil column were converted to: N1 = 0.6g, N2 = 0.75g (applied in the form of urea). There were a total of 6 treatments (W1N1, W1N2, W2N1, W2N2, W3N1, W3N2), with 3 replicates per group. Experimental procedure: After the soil column is saturated, let it stand for 24 hours to remove gravity water; Apply nitrogen according to the treatment (dissolve urea in irrigation water and spray it evenly on the surface of the soil column), and then pour in the set amount of water at once; Collect wastewater from the collection bottle and test for NO3 at 0, 1, 3, 5, 7, and 10 days after irrigation. - -N, NH4 + -N concentration (ion chromatography), calculate nitrogen leaching amount (unit g / column). Results analysis: Under the same nitrogen application rate, increased irrigation volume leads to increased nitrogen leaching: In the N1 treatment, W2 leaching increased by 12-14% compared to W1, and W3 increased by 13-15% compared to W2; in the N2 treatment, W2 increased by 14-16% compared to W1, and W3 increased by 15-18% compared to W2. Under the same irrigation level, increasing nitrogen application leads to increased nitrogen leaching: in treatment W1, N2 increased by 18-20% compared to N1; in treatment W2, N2 increased by 20-22% compared to N1; and in treatment W3, N2 increased by 22-24% compared to N1. Conclusion: For every 150m increase in irrigation volume... 3 hm -2 Nitrogen leaching increases by 12-18%, and for every 75 kg increase in nitrogen application rate... hm -2 The nitrogen leaching amount increased by 18-24%, indicating that irrigation activities are the core driving factor of nitrogen migration (contribution rate > 50%), while the amount of nitrogen applied is a secondary factor.

[0028] S4. Analysis of the Spatiotemporal Distribution Characteristics of Pollution Load and Delineation of Core Area 4.1 Pollution Level Assessment The Nemerow Comprehensive Pollution Index method is used for calculation. The formula is as follows: ; Wherein: PI represents the Nemerow Comprehensive Pollution Index, a dimensionless index used to characterize the overall pollution level of water bodies. The classification standards are: PI ≤ 1.0 (clean), 1.0 < PI ≤ 2.0 (light to moderate pollution), PI > 2.0 (heavy pollution). avg This represents the arithmetic mean of the pollution indices of various pollutants; it is dimensionless. (n is the number of evaluation factors; in this invention, n=9); PI max PI represents the maximum value among the pollution indices of various pollutants; it is dimensionless. i The pollution index of the i-th pollutant is dimensionless. Ci represents the measured concentration of the i-th pollutant (unit: mg / L, PFAS converted to mg / L), S i The value represents the evaluation standard value of the i-th pollutant (unit: mg / L). Groundwater adopts GB / T14848-2017 Class III, and surface water adopts GB3838-2022 Class III. Calculation results: The PI value of water bodies in the Yellow River South Irrigation Area ranges from 1.2 to 1.8, indicating an overall level of light to moderate pollution. Among them, the PI of groundwater is 1.5-1.8 (high risk of severe pollution), the PI of surface water is 1.2-1.5, the PI of drained water is 1.4-1.7, the PI of canal water is 1.1-1.4, and the PI of drainage water is 1.6-1.8, indicating that the risk of groundwater pollution is significantly higher than that of surface water.

[0029] 4.2 Delineation of the Core Area Pollution load calculation model: The pollution load per unit area (kg / hm²) is calculated using the "source strength-migration" model. The formula is as follows: ; Among them: Q i F represents the pollution load per unit area (kg / hm²) of the i-th pollutant; i The source intensity (mg / hm²) of the i-th pollutant is represented by the nitrogen application rate (mg / hm²) for nitrogen, the phosphorus application rate (mg / hm²) for phosphorus, and the soil background value × soil bulk density × topsoil thickness (20cm) for other pollutants; α i β represents the leaching / runoff coefficient (dimensionless) of the i-th pollutant, obtained through indoor soil column experiments: nitrogen α = 0.25-0.35, phosphorus α = 0.05-0.15, fluoride α = 0.15-0.25; Ui represents the amount of the i-th pollutant absorbed by the crop (mg / hm²), with nitrogen taken from the sunflower's nitrogen requirement (250 kg / hm²) and phosphorus taken from the sunflower's phosphorus requirement (50 kg / hm²); i Indicates crop absorption efficiency (dimensionless), nitrogen β = 0.6-0.7, phosphorus β = 0.5-0.6; V i This represents the amount of pollutant adsorbed by the soil for the i-th pollutant (mg / hm²), obtained through soil adsorption experiments: phosphorus V = 1000 kg / hm², fluoride V = 500 kg / hm²; γ i The values ​​represent soil adsorption efficiency (dimensionless), with phosphorus γ = 0.8-0.9 and fluoride γ = 0.7-0.8. Evaluation Unit Division: Based on 30m DEM data, the Yellow River South Irrigation Area (500km²) was divided into 1km×1km grid evaluation units, totaling 500 units (numbered 1-500). The geometric center coordinates (X, Y, F, G, I, Z) of each unit were determined. i ,Y i The data was extracted using the Gaussian Cartesian coordinate system (Xi'an 80 coordinate system) and ArcGIS. Center of gravity migration model calculation: Calculate the center of gravity coordinates (X, Y) of pollution load for each season (spring irrigation in April, autumn irrigation in September, and winter irrigation in November), formula: ; Where: X and Y represent the abscissa and ordinate of the centroid of the pollution load, respectively (unit: m); Q i This represents the total pollution load (in kg) of the i-th evaluation unit. (TN, TP, F) - As priority control factors, with weights of 0.5, 0.3, and 0.2 respectively; X i ,Y i These represent the x-coordinate and y-coordinate (unit: m) of the geometric center of the i-th evaluation unit, respectively. Core area determination: The center of gravity migration trajectory for each season was plotted using ArcGIS, combined with the data from each unit Q. iSize (by Q) i Based on the ranking, the top 20% are high-load units. The eastern part of the irrigation basin (unit numbers 120-150, coordinate range E108.4°-108.6°, N40.7°-40.9°) and the western part of the Dugui-Hangjin irrigation basin (unit numbers 210-240, coordinate range E108.1°-108.3°, N40.5°-40.7°) are identified as core pollution areas. This area has a Q... i This accounts for 45-50% of the total load in the irrigation area, mainly due to intensive agricultural activities (sunflower cultivation accounts for 80%) and high irrigation volume (1350m³). 3 hm -2 Large amount of fertilizer applied (375kg) hm -2 ).

[0030] 4.3 Relationship of Driving Factors Pearson linear correlation analysis was used to explore the relationship between irrigation amount (W), nitrogen application amount (N), rainfall (P) and total pollution load (Q). The steps were as follows: Data collection: Monthly W (m³) data for the Yellow River South Irrigation District from 2021 to 2024 were collected. 3 hm -2 (Source: Irrigation District Administration Bureau), N (kg) hm -2 (Source: Farmer Survey), P (mm, Source: Weather Station), Q (kg / hm) 2 (Calculated using the above model), a total of 48 sets of data; Correlation analysis: The correlation coefficient R was calculated using SPSS 26.0 software, with a significance level of α=0.05. Results: The correlation coefficient between W and Q was R = 0.85 (P < 0.01), and the linear regression equation was Q = 0.02W + 15 (R² = 0.72). The correlation coefficient between N and Q was R = 0.80 (P < 0.01), and the linear regression equation was Q = 0.03N + 5 (R² = 0.68). The correlation coefficient between P and Q is R = 0.74 (P < 0.01), and the linear regression equation is: Q = 0.5P + 20 (R² = 0.55). Contribution rate calculation: Principal component analysis (PCA) was used to calculate the contribution rate of each driving factor. Two principal components were extracted (cumulative variance contribution rate 90%). The first principal component (variance contribution rate 60%) was mainly contributed by W and N, and the second principal component (variance contribution rate 30%) was mainly contributed by P. Combined with seasonal changes (spring irrigation - autumn irrigation - winter irrigation, Q continuously increases), the contribution rate of irrigation activities was determined to be 45%, the contribution rate of seasonal changes (temperature and precipitation changes) was 30%, the contribution rate of fertilizer application was 20%, and the contribution rate of other factors (topography and soil) was 5%. Irrigation activities and seasonal changes are the main driving factors.

[0031] S5. Pollutant Migration Patterns and Causal Analysis 5.1 Migration Trend Detection Monitoring section layout: Ten monitoring sections will be set up along the Yellow River runoff direction (from upstream to downstream: Yellow River main stream → irrigation district branch canals → irrigation district farmland → drainage canals → Yellow River main stream), specifically: Upstream sections (S1-S3): Yellow River main stream entry section (S1), upstream branch canal entrance (S2), upstream farmland drainage outlet (S3), slope 5-8°; Midstream sections (S4-S7): Midstream branch canal (S4), midstream farmland drainage outlet (S5), midstream main drainage canal (S6), and eastern drainage outlet of the irrigation area (S7), with a slope of 2-5°; Downstream sections (S8-S10): the end of the downstream branch canal (S8), the downstream farmland drainage outlet (S9), and the confluence of the main drainage canal and the Yellow River (S10), with a slope of 0-2°; Monitoring and Analysis: Water samples are collected monthly from each cross-section to test total nitrogen (TN) and nitrite nitrogen (NO2). - -N), nitrate nitrogen (NO3) - -N), ammonia nitrogen (NH4) + -N), fluoride ions (F) - Concentration, continuously monitored for 1 year; Migration trends: Spatial trend: From upstream (S1) to downstream (S10), pollutant concentrations show a "decline followed by an increase" trend. Concentrations in S1-S3 decrease (due to dilution in the Yellow River main stream), while concentrations in S3-S10 increase. The concentration in S10 is significantly higher than that in S1: TN concentration S1 = 1.2 mg / L, S10 = 3.5 mg / L (an increase of 191.7%); F - Concentrations S1=0.5mg / L, S10=1.8mg / L (an increase of 260%), showing a reverse migration trend (downstream concentration is higher than upstream). Time trend: The reverse migration was most significant after autumn irrigation (September), with TN concentration S10 being 220% higher than S1; followed by after spring irrigation (April), 180% higher; after winter irrigation (November), 190% higher; and not significant during the non-irrigation period (January-March), but 50% higher.

[0032] 5.2 Construction of Coupled Transfer Model By integrating HYDRUS-1D (soil infiltration), SWAT (surface runoff), and a self-developed "irrigation backwater sub-model," a coupled migration model of "irrigation backwater-soil infiltration-surface runoff" was constructed, specifically: 5.2.1 Soil Leakage Sub-model (based on HYDRUS-1D4.16) Model parameters: Input soil data collected by S1, soil texture is divided into 3 layers: 0-20cm: Loam, field water holding capacity θ fc =28-32%, saturated moisture content θ s =40-45%, saturated hydraulic conductivity K s =1.5×10 -6 m / s, the soil moisture characteristic curve was obtained using the van Genuchten equation: , (θ r =Residual moisture content, α = 0.015cm -1 (n=1.8, m=1-1 / n). 20-60cm: clay loam, θ fc =30-34%, θ s =42-47%, K s =1.2×10 -6 m / s, α = 0.012 cm -1 n=1.6, m=1-1 / n; 60-100cm: Clay, θ fc =32-36%, θ s =45-50%, K s =1.0×10 -6 m / s, α = 0.010 cm -1 n=1.4, m=1-1 / n; Boundary conditions: The upper boundary is a variable head (irrigation volume + rainfall), and the lower boundary is free drainage; Simulation output: Soil seepage rate (unit: m / d), soil moisture content (volume fraction) at different depths, pollutants (TN, F) - Leakage flux (unit: mg / (m³)) 2 d)).

[0033] 5.2.2 Surface runoff sub-model (based on SWAT2012) Model parameters: Input S1 data: Land use data: Based on remote sensing inversion, cultivated land (sunflower) accounts for 65%, grassland 25%, water area 5%, and other 5%; Soil data: Soil type map (from "China Soil Map" 1:1,000,000), soil organic carbon content 1.2-1.5%; Meteorological data: daily precipitation, daily average temperature, wind speed, and relative humidity; Runoff calculation: The SCS curve number method (CN value) is used to determine the CN value according to the land use type and soil hydrological group (A / B / C / D): cultivated land (sunflower) CN=70-75, grassland CN=60-65, water area CN=90-95; Simulation output: Surface runoff (unit: m³) 3 / hm 2 ), runoff coefficient (dimensionless), pollutants (TN, TP) runoff (unit: mg / (m³) 2 d)).

[0034] 5.2.3 Irrigation return water sub-model Model principle: Based on water balance, it calculates the amount of water that is not absorbed by the soil or evaporated after irrigation (i.e., irrigation return water). Formula: ; The definitions and acquisition methods of each parameter are specified as follows: Q 回 Q represents the irrigation return water volume (unit: m³ / hm²), which is the core output of the model and characterizes the remaining irrigation water that can flow into ditches; 灌 Irrigation water volume (unit: m³) 3 / hm 2 The data comes from field data collected in step S1 (such as the 900-1350m range in sunflower irrigation systems). 3 hm -2 (gradient value); Q 渗 This represents soil infiltration (unit: m³ / hm²), directly calling the simulation output of the HYDRUS-1D model (calculated based on parameters such as soil texture and field capacity); Q 蒸 Soil evaporation (unit: m³ / hm²) is expressed using the internationally accepted Penman-Monteith formula, and its modified standard form is as follows: Please complete the definitions and units of each parameter in the formula: R n Net surface radiation (unit: MJ / (m²)) 2 d) is calculated by subtracting reflection and longwave radiation losses from solar radiation data observed at meteorological stations; G represents soil heat flux (unit: MJ / (m²)). 2 d) The soil temperature gradient method was used for actual measurement (e.g., the temperature difference between the 0-5cm and 5-10cm soil layers was converted). The slope of the saturated vapor pressure curve (unit: kPa / ℃) is calculated from the daily average temperature using an empirical formula. ); γ represents the ecliptic constant (unit: kPa / ℃), taken as 0.066 kPa / ℃ under standard atmospheric pressure; u2 represents the wind speed at a height of 2m (unit: m / s), directly observed data from the meteorological station; T represents the daily average temperature (unit: ℃), observed data from the meteorological station; e s This represents saturated vapor pressure (unit: kPa), calculated from the daily average temperature. e a The actual water vapor pressure (unit: kPa) is calculated by multiplying the relative humidity and the saturated water vapor pressure. (RH is relative humidity), Q 作 This represents crop transpiration (unit: m³ / hm²), calculated using the FAO-recommended Penman-Monteith derivative formula, which has been corrected to: , where: K c Q represents the crop coefficient (dimensionless), with values ​​assigned to different growth stages of sunflower: emergence 0.3, jointing stage 0.7, budding stage 1.1, grain-filling stage 1.0, maturity stage 0.6; 蒸 This represents the same soil evaporation calculation result as above (unit: m³ / hm²).

[0035] Backwater migration calculation: Considering the slope differences of ditches within the irrigation area (1‰ slope for suction pipes and 0.8‰ slope for collection pipes), the Manning formula is used to calculate the flow velocity of backwater in the ditches. The corrected standard form is as follows: Where v: the velocity of the return water in the ditch (unit: m / s); n: the roughness coefficient of the ditch (dimensionless), 0.013 for concrete-lined ditches and 0.022 for earthen ditches; R: the hydraulic radius of the ditch (unit: m), calculated from the cross-sectional dimensions of the ditch (rectangular ditch). (B is the width of the ditch, H is the water depth); S: ditch slope (dimensionless), 0.001 for the suction pipe and 0.0008 for the collection pipe; The return water migration time is further calculated based on the flow velocity. L is the length of the ditch (m) and the change in pollutant concentration ( C 灌为 Irrigation water pollutant concentration, C 土 Q represents the concentration of pollutants adsorbed in the soil. 吸 (This is to simulate the amount of water adsorbed by the soil) to achieve quantitative simulation of the water migration process.

[0036] 5.2.4 Model Coupling and Validation Coupling Interface: A coupling interface was written in C++ to implement data interaction among the three models: ① Q output from HYDRUS-1D 渗 ① The pollutant leakage flux is input into the coupling interface; ② The surface runoff and pollutant runoff output by SWAT are input into the coupling interface; ③ The Q calculated by the irrigation return water sub-model is input into the coupling interface. 回 The concentration of pollutants in the returned water is transmitted to the coupling interface; Pollutant migration coefficient: calibrated through indoor soil column experiments and field monitoring, and set as follows: Pollutant migration coefficient f in leakage water 渗 =0.8 (meaning 80% of the contaminants migrate with the leaked water). Pollutant transport coefficient f in surface runoff 径 =0.9 (90% of pollutants migrate with surface runoff). Pollutant migration coefficient f in irrigation return water 回 =0.95 (95% of pollutants migrate with irrigation return water). Model validation: Field monitoring data from April to September 2024 (TN, F at each cross section) were used. - (Concentration) Validate the model and calculate the coefficient of determination R. 2 With root mean square error (RMSE), results: simulated TN value vs. measured value R 2 =0.85, RMSE=0.3mg / L; F - Simulated value vs. measured value R 2 =0.82, RMSE=0.2mg / L, the verification accuracy meets the requirements (R 2 >0.8, RMSE <0.5mg / L); Cause analysis: Model simulations show that the main causes of reverse migration are: ① Irrigation backflow contributes 50% (Q_backflow carries a large amount of pollutants, flows from farmland to drainage canals, and eventually flows into the lower reaches of the Yellow River); ② Spatial differences in the densely agricultural areas in the middle reaches contribute 35% (large fertilizer application and low crop absorption in the middle reaches result in a large surplus of pollutants, which migrate downstream with runoff / leakage); ③ Topography contributes 15% (the downstream slope is gentle, the water flow is slow, and pollutants are easily accumulated).

[0037] S6. Pollution Load Clustering Type Classification The following steps were taken using R-mode clustering analysis: 6.1 Variable Selection and Data Preprocessing Clustering variables: The concentration data of 9 pollutants obtained from S1 to S5 were selected as clustering variables, namely X1=TN, X2=NH4 + -N、X3=NO2 - -N、X4=NO3 --N, X5=TP, X6=TDS, X7=COD, X8=PFAS, X9=F - ; Data source: Monthly monitoring data, totaling 12 months × 30 sampling points = 360 sets of data; Data standardization: Due to the different dimensions of the variables, min-max standardization (to eliminate the influence of dimensions) is used. Formula: ; in This represents the standardized value of the i-th variable in the j-th group of data (range 0-1). This represents the original value of the i-th variable in the j-th data set; This represents the minimum value of the i-th variable; This represents the maximum value of the i-th variable; Outlier handling: Outliers (data exceeding 1.5 times the interquartile range) were removed using box plots. A total of 12 outliers were removed, resulting in 348 valid data sets.

[0038] 6.2 Clustering Methods and Parameter Settings Distance calculation: Euclidean distance (measures similarity between variables) is used, formula: ; in: This represents the Euclidean distance between the k-th variable and the l-th variable; the smaller the distance, the higher the similarity. This represents the standardized value of the j-th group of data for the k-th variable; This represents the standardized value of the j-th group of data for the l-th variable; Clustering method: Ward clustering method (which minimizes the increase of the sum of squares within each cluster and achieves better clustering results) was adopted and implemented using SPSS 26.0 software; Clustering threshold: The clustering threshold is determined by principal component analysis (PCA). Principal components are extracted, and the variance contribution rate of each principal component is calculated. Clustering is stopped when the cumulative variance contribution rate is >85%. The final number of clusters is determined to be 4.

[0039] 6.3 Clustering Results and Validation Clustering results: Cluster 1 (reducing nitrogen): X2 (NH4) + -N), X3 (NO2) - -N), Euclidean distance 0.35-0.45, similarity 0.85-0.90, reason: NH4 + -N and NO2 - -N are all reducing nitrogen, mainly from manure wastewater, and their migration paths are similar (through soil pore water seepage). Cluster 2 (Total Phosphorus-TDS-COD): X5 (TP), X6 (TDS), X7 (COD), Euclidean distance 0.40-0.50, similarity 0.80-0.85. Reason: TP is adsorbed on soil particles, while TDS and COD mainly come from soil dissolution and organic matter decomposition, and both migrate with surface runoff. Cluster 3 (perfluorinated pollutants): X8 (PFAS), forming a separate cluster, Euclidean distance 0.60-0.70, similarity 0.75-0.80, reason: PFAS mainly comes from pesticides and plastic films, and has a unique migration path (difficult to degrade, slowly migrating with groundwater). Cluster 4 (Total Nitrogen): X1 (TN), X4 (NO3) - -N), Euclidean distance 0.30-0.40, similarity 0.88-0.92, reason: TN with NO3 - -N is the main form (accounting for 60-70% of TN), mainly derived from fertilizer and manure, and migrates with irrigation return water; Validation: The silhouette coefficient was used for validation. The formula is: ; Where: si represents the silhouette coefficient of the i-th variable (range -1 to 1, the closer to 1, the better the clustering effect); ai represents the average distance between the i-th variable and other variables in the same category; bi represents the average distance between the i-th variable and the nearest outlier variable; Calculation results: The average silhouette coefficient of the four clusters is 0.78 > 0.7, indicating that the clustering results are effective and the variable classification is reasonable.

[0040] S7. Spatiotemporal Differentiation Quantitative Analysis of Pollutant Loss 7.1 Invest Model Optimization The Invest model (version 3.8.9) was used to optimize the irrigation system for the Yellow River basin. The steps were as follows: 7.1.1 Basic Parameter Input Input the basic data collected by S1: DEM data: 30m resolution, extracting slope (upstream 5-8°, midstream 2-5°, downstream 0-2°), aspect, and catchment area; Meteorological data: Daily precipitation, temperature, wind speed, and relative humidity over the past 5 years; calculate the average monthly rainfall (20-50 mm, April 20 mm, July 50 mm, November 30 mm). Soil data: soil type, organic carbon content, field water holding capacity, saturated hydraulic conductivity (data collected in S1). Land use data: types of cultivated land, grassland, and water bodies retrieved from remote sensing, with sunflower planting areas separately classified within cultivated land; Management data: Irrigation system (spring irrigation in April, autumn irrigation in September, winter irrigation in November), fertilization system (nitrogen application rate 300-375 kg) hm -2 Phosphorus application rate: 100-150 kg hm -2 ).

[0041] 7.1.2 Added a "Topographic Water Collection-Irrigation and Fertilization Overlay Drive Module" A self-developed module (written in Python 3.9 and integrated into the Invest model), with the following core functions: Topographical drainage calculation: Based on DEM data, the D8 algorithm is used to calculate the drainage area of ​​each grid and divide the drainage area into upstream drainage area, midstream drainage area, and downstream drainage area. The pollutant loss within the drainage area is weighted averaged according to the drainage area. Analysis of the superposition of irrigation and fertilization: Irrigation amounts (W1=900, W2=1125, W3=1350m³) were considered. 3 hm -2 ) and nitrogen application rate (N1=300, N2=375kg) hm -2 The mixture is divided into 6 treatments. The pollution loss driving intensity (D) for each treatment is calculated using the following formula: ; Where: D: Pollution runoff driving intensity (dimensionless, range 0-1); W: Actual irrigation volume (m³) 3 hm -2 W max Maximum irrigation volume (1350m³) 3 hm -2 N: Actual nitrogen application rate (kg) hm -2 ); N max Maximum nitrogen application rate (375 kg) hm -2 ); Module output: Pollution loss driving intensity for each catchment area and each treatment, used as input parameters for the Invest model to correct the loss calculation results.

[0042] 7.2 Time Differentiation Analysis The temporal variations of total nitrogen (TN) and total phosphorus (TP) losses in the Yellow River South Irrigation District in 2024 were obtained through simulation using the optimized Invest model: TN loss: April (spring irrigation) 1.25 kg / ha, May 1.50 kg / ha, June 1.80 kg / ha, July 2.20 kg / ha, August 2.80 kg / ha, September (autumn irrigation) 3.50 kg / ha, October 3.80 kg / ha, November (winter irrigation) 4.15 kg / ha, December 3.90 kg / ha, January 2025 3.50 kg / ha, February 3.00 kg / ha, March 2.50 kg / ha, with the peak occurring in November (4.15 kg / ha). TP loss: April 0.20 kg / ha, May 0.25 kg / ha, June 0.30 kg / ha, July 0.35 kg / ha, August 0.45 kg / ha, September 0.55 kg / ha, October 0.65 kg / ha, November 0.77 kg / ha, December 0.70 kg / ha, January 2025 0.60 kg / ha, February 0.50 kg / ha, March 0.40 kg / ha, with the peak occurring in November (0.77 kg / ha). Driving factors: ① Irrigation and fertilization: The loss increases significantly after spring, autumn and winter irrigation (increased irrigation leads to enhanced leaching / runoff, and increased fertilization leads to increased pollutant source intensity); ② Rainfall events: High rainfall in July and August (average monthly 50mm) enhances surface runoff and increases loss; ③ Topographical catchment: The large catchment area in the middle reaches of the river leads to the accumulation of pollutants, which, combined with irrigation and fertilization, results in a continuous increase in loss.

[0043] 7.3 Spatial Differentiation Analysis Model simulations combined with field monitoring yielded the spatial distribution of pollution loss: Upstream area (slope 5-8°): TN loss is 1.5-2.0 kg / ha, TP loss is 0.25-0.35 kg / ha, accounting for 15-20% of the total loss; only in some areas adjacent to ditches (<50m from the ditch) is the loss higher (TN 2.0-2.5 kg / ha, TP 0.35-0.45 kg / ha), due to the scouring effect of the ditch water flow, which carries away pollutants from the surrounding soil; Midstream area (slope 2-5°): TN loss 3.0-3.8 kg / ha, TP loss 0.50-0.65 kg / ha, loss rate 45-50%, classifying it as a high-risk area; due to intensive agricultural activities (sunflower cultivation accounts for 80%) and high irrigation volume (1125-1350m). 3 hm -2 Large amount of fertilizer (350-375kg) hm -2 The area has high pollutant source intensity, flat terrain, large catchment area, and high loss rate. Downstream area (slope 0-2°): TN loss is 2.5-3.0 kg / ha (before autumn irrigation) and 3.8-4.15 kg / ha (after autumn irrigation), TP loss is 0.40-0.50 kg / ha (before autumn irrigation) and 0.65-0.77 kg / ha (after autumn irrigation), accounting for 30-35% of the total loss. After autumn irrigation, the loss intensity increases dramatically, with TN exceeding that of the upstream area by 20-30%. This is because the downstream terrain is gentle (slope <2°), the water flow velocity is slow (0.1-0.2 m / s), pollutants are easily accumulated, and a large amount of autumn irrigation backflow carries pollutants from the midstream, leading to an increase in loss.

[0044] 7.4 Differences in nitrogen and phosphorus migration The differences in the migration of TN and TP were investigated using an acrylic soil column experiment (10 cm in diameter and 100 cm in height). Experimental Design: Soil filling: The irrigated area soil is filled in layers (0-20cm topsoil layer, 20-60cm plow pan, and 60-100cm parent material layer), with bulk densities of 1.3, 1.4, and 1.5 g / cm³, respectively; Treatment settings: Irrigation volume 1125m³ 3 hm -2 Nitrogen application rate: 350 kg hm -2 Phosphorus application rate: 120 kg hm -2 Repeated 3 times; Sampling: Soil samples were collected at different depths (20, 40, 60, 80, 100 cm) of the soil column at 0, 1, 3, 5, 7, and 10 days after irrigation to detect TN and TP concentrations; leachate was collected from the bottom of the soil column to detect TN and TP concentrations.

[0045] Migration path analysis: TN migration: mainly driven by runoff migration (60% contribution) and seepage migration (40% contribution); within runoff migration, NO3... - -N accounts for 70-80%, due to NO3- - -N is readily soluble in water and is lost through surface runoff; during seepage and migration, NH4+... + -N accounts for 30-40%, due to NH4 + -N is adsorbed by soil colloids and slowly seeps out with soil pore water; TP migration: mainly due to adsorption on soil particles and loss with sediment (contribution rate 85%), while dissolved TP accounts for only 15%; because the adsorption coefficient of TP is Kd=500-800mL / g (far higher than Kd=50-100mL / g of TN), it is easily adsorbed by iron and aluminum oxides in the soil, and its mobility is weaker than that of TN. The upstream steep slope area has special characteristics: the upstream slope is 5-8°, the surface runoff velocity is fast (0.3-0.4m / s), and the scouring effect on the soil is strong, which leads to the washing away of soil particles that adsorb TP. The TP runoff contribution rate is 30% (significantly higher than 20% in the middle reaches and 15% in the downstream), indicating that the runoff in the upstream steep slope area has a significant TP scouring effect.

[0046] S8. Source analysis of chemical components and assessment of water quality suitability 8.1 Source Analysis (SOM-PMF Combined) The sources of chemical components in canal water and drainage were analyzed using a self-organizing neural network (SOM) combined with positive definite matrix factor analysis (PMF). The steps are as follows: 8.1.1 Data Preparation Sample collection: Water samples (main and branch irrigation canals, 30 samples) and drainage samples (farmland underground drainage outlets, 30 samples) were collected from the Yellow River South Irrigation Area. Chemical composition was analyzed: X1=TDS, X2=Ca. 2+ X3=Mg 2+ X4=Na + X5=K + X6=Cl - X7=SO4 2- X8=HCO3 - X9=TN, X10=TP, X11=F - (Unit: except for TN and TP, which are in mg / L, the others are in mg / L); Data preprocessing: Remove missing values ​​(<5%), use logarithmic transformation (ln(X+1)) to reduce data skewness, and standardize (z-score standardization: (μ is the mean, σ is the standard deviation).

[0047] 8.1.2 Dimensionality Reduction of Self-Organizing Neural Networks (SOM) Using MATLAB R2023a's SOMToolbox, parameter settings are as follows: Network structure: Input layer neurons = 11 (corresponding to 11 chemical components), hidden layer uses 4×4 hexagonal grid (16 neurons), output layer neurons = 16; Training parameters: The initial learning rate is 0.05, which decreases linearly with the number of iterations (down to 0.01 after 1000 iterations). The neighborhood function is a Gaussian function with an initial neighborhood radius of 5, which decreases to 1 after 1000 iterations. Training results: The 11-dimensional chemical composition data was reduced to 6 feature vectors (cumulative variance contribution rate of 92%), which are as follows: Feature vector 1 (F1): X1, X6, X7 (TDS, Cl) - SO4 2-), load factor > 0.8; Eigenvector 2 (F2): X2, X3, X8 (Ca 2+ Mg 2+ HCO3 - ), load factor > 0.8; Eigenvector 3 (F3): X9, X10 (TN, TP), loading coefficient > 0.8; Eigenvector 4 (F4): X11 (F - ), load factor > 0.7; Eigenvector 5 (F5): X4 (Na) + ), load factor > 0.7; Feature vector 6 (F6): X5 (K) + ), load factor > 0.7.

[0048] 8.1.3 Positive Definite Matrix Factor Analysis (PMF) Using EPAPMF 5.0 software, input 6 feature vectors, the steps are as follows: Factor number determination: Set the number of factors from 1 to 5, calculate the sum of squared residuals (RSS) and Q value (the ratio of actual residuals to theoretical residuals) for each factor number. When the number of factors = 3, the RSS is the smallest (0.12) and the Q value is closest to 1 (1.05). Therefore, the number of factors is determined to be 3. Factor loadings and source identification: Factor 1 (F1): F1 and F1(TDS, Cl) - SO4 2- The loading coefficient is 0.92, corresponding to the source of evaporite dissolution (gypsum, rock salt, and other evaporites distributed around the irrigation area, dissolving and releasing Cl). - SO4 2- (Increase TDS) Factor 2 (F2): F2 and F2 (Ca) 2+ Mg 2+ HCO3 - The loading factor is 0.88, corresponding to the source of silicate weathering (the parent material of the irrigation area soil is silicate rock, and weathering releases Ca). 2+ Mg 2+ It reacts with CO2 to produce HCO3. - ); Factor 3 (F3): F3 and F3(TN, TP), F4 (F - The loading factor is 0.90, corresponding to agricultural input sources (TN and TP released by fertilizer application, and fluoride released from the soil with irrigation water).

[0049] Contribution percentage calculation: Canal water: F1 (evaporite dissolution) contributes 35-40%, F2 (silicate weathering) contributes 25-30%, and F3 (agricultural input) contributes 30-35%; Drainage: F1 contributes 40-45%, F2 contributes 20-25%, F3 contributes 35-40%; Error analysis: The bootstrap method (100 iterations) was used to calculate the contribution ratio error. The error was ≤5%, and the result was reliable.

[0050] 8.2 Water quality suitability assessment (entropy weight method) Constructing an evaluation index system for "drinkability" and calculating the comprehensive index using the entropy weight method, the steps are as follows: 8.2.1 Evaluation Indicators and Classification Ten indicators were selected: X1 = pH (6.5-8.5 is suitable), X2 = DO.

[0051] The above are merely preferred embodiments of the present invention, but the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and inventive concept of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A method for spatiotemporal analysis of non-point source pollution load in Yellow River irrigation area, characterized in that, The method comprises the following steps: S1, multi-source basic data collection and preprocessing: collecting agricultural non-point source pollution basic data, Yellow River mainstream and irrigation area water quality monitoring data, using unmanned aerial vehicle multispectral technology to obtain dynamic data of surface coverage in the irrigation area, using soil pore water in-situ collection-detection integrated device to collect soil pore water layer by layer and detect pollutant concentration, and collecting meteorological data and DEM terrain data at the same time; S2, main pollutant identification and priority control factor determination: adopting fuzzy comprehensive evaluation method, constructing a three-level evaluation index system of "pollutant concentration-ecological risk-human health risk", calculating the membership degree and weight of each pollutant based on the water quality monitoring data obtained in S1, identifying the main pollutants and priority control factors of the water body in the irrigation area; S3, core pollutant source analysis and migration influence mechanism analysis: taking autumn irrigation in the irrigation area as the research object, collecting water samples, using nitrogen and oxygen isotope ratio analysis method combined with real-time monitoring chip of denitrifying functional bacteria group activity, tracing the main source of nitrogen pollution, analyzing the intensification mechanism of denitrification of underground water caused by the enhanced activity of denitrifying functional bacteria group after autumn irrigation, and determining the driving mechanism of nitrogen pollution migration caused by irrigation activities through indoor simulation experiment; S4, analysis of spatial and temporal distribution characteristics of pollution load and definition of core area: adopting Nemerow comprehensive pollution index method to calculate the pollution index of different water bodies in the irrigation area, determining the pollution level and risk difference, calculating the pollution load barycentric coordinates of each sub-region in different seasons based on the pollution load calculation model combined with the barycentric migration model, defining the core pollution area, and establishing the correlation model of irrigation amount, fertilizer amount, rainfall and pollution load; S5, analysis of pollutant migration law and cause: detecting the pollutant concentration of different sections along the runoff direction of the irrigation area, analyzing the change trend of the pollutants along the runoff direction, combining the remote sensing data and soil pore water detection data of S1, constructing a coupled migration model of "irrigation return water-soil seepage-surface runoff", and inferring the main cause of the reverse migration of the pollutants; S6, pollution load clustering type division: adopting R-type cluster analysis method, taking the pollutant concentration data obtained in S1-S5 as variables, setting the clustering threshold based on the variance contribution rate of more than 85%, calculating the similarity between variables through Euclidean distance, and dividing the pollution load clustering type; S7, quantitative analysis of spatial and temporal differentiation of pollutant loss amount: based on the Invest model, taking the DEM data, meteorological data, soil data and crop planting data of S1 as inputs, adding "topographic water catchment-irrigation and fertilization superposition driving module", simulating the loss amount of total nitrogen and total phosphorus in the south irrigation area, analyzing the time differentiation and spatial differentiation of the loss amount, and distinguishing the nitrogen and phosphorus migration paths through soil column experiment; S8, chemical component source analysis and water quality suitability evaluation: adopting self-organizing neural network combined with positive matrix factor analysis method, constructing a "component-source" mapping model, analyzing the main source and contribution proportion of chemical components in canal water and drainage, adopting entropy weight water quality evaluation method to construct a "drinking suitability" evaluation index system, calculating the water quality comprehensive index, and determining the proportion of samples not suitable for drinking.

2. The method according to claim 1, wherein, The basic data of agricultural non-point source pollution in S1 includes crop planting structure, growth period of sunflower and other crops, fertilizer type and application amount, irrigation system and irrigation water volume, and water body in the irrigation area includes groundwater, dry water body, surface water, canal water and drainage water. The water quality monitoring data are detected by using high performance liquid chromatograph and ion chromatograph, the monitoring frequency is once a month, and 3 parallel samples are collected each time, and the relative deviation of the parallel samples is controlled to be less than 5%.

3. The method according to claim 1, wherein the method is characterized in that, The unmanned aerial vehicle multi-spectral technology in S1 adopts a six-rotor unmanned aerial vehicle carrying a 5-band hyperspectral imager, the flight height is 100 m, the spatial resolution is 0.5 m, the flight is performed twice a month, the flight track covers the whole irrigation area, the obtained remote sensing data are subjected to atmospheric correction and geometric correction by using ENVI software, and then NDVI index and soil water content are obtained by inversion. The soil pore water in-situ collection-detection integrated device includes a stainless steel sampling tube with a length of 120 cm, a 0.45 mu m mixed fiber filter film wrapped on the tube wall and an in-situ ion sensor embedded in the tube. The layered collection depth is 0-20 cm, 20-60 cm and 60-100 cm, 3 times of collection are performed at each sampling point, and the average value is taken as the detection result.

4. The method according to claim 1, wherein, The specific steps of the fuzzy comprehensive evaluation method in S2 are as follows: S21, determining the evaluation factor set U = {total nitrogen (u1), ammonia nitrogen (u2), nitrite nitrogen (u3), nitrate nitrogen (u4), total phosphorus (u5), TDS (u6), chemical oxygen demand (u7), perfluorinated pollutants (u8), fluoride (u9)}; S22, determining the evaluation grade set V = {clean (v1), slight pollution (v2), moderate pollution (v3), severe pollution (v4)}; S23, inviting 5 experts in the field of environmental engineering to construct a judgment matrix by using the analytic hierarchy process, calculating the weight vector W = {w1, w2, …, w9} of each evaluation factor, and verifying the rationality of the weight through consistency test, wherein the weight of total nitrogen w1 = 0.25, the weight of ammonia nitrogen w2 = 0.20, and the weight of fluoride w9 = 0.15; S24. Construct a fuzzy relation matrix R, where rij represents the membership degree of the i-th evaluation factor to the j-th evaluation level. The membership degree is calculated using a descending half-trapezoidal distribution function. When the factor concentration ≤ the cleanliness standard, the membership degree r... i1 =1、r i2 =r i3 =r i4 =0, when the concentration is between the clean and slightly polluted standards, r i1 =(Mild Standard - Concentration) / (Mild Standard - Cleanliness Standard), r i2 =1-r i1 r i3 =r i4 =0; S25, calculating the fuzzy comprehensive evaluation result B = W x R by using the weighted average method, and determining the main pollutants and priority control factors according to the maximum membership degree principle.

5. The method according to claim 1, wherein, The specific steps of the nitrogen-oxygen isotope ratio analysis method in S3 are as follows: S31, collecting 500 mL of water sample, filtering through a 0.45 mu m filter film, adding 10 mL of concentrated sulfuric acid to acidify to pH < 2, and separating nitrogen in the sample by using a distillation method to convert the nitrogen into silver nitrate precipitate; S32, measuring δ 15 N and δ 18 O in the silver nitrate precipitate by using a stable isotope ratio mass spectrometer, wherein the test precision of δ 15 N is ±0.2‰, and the test precision of δ 18 O is ±0.3‰, with international standard substance AIR, VSMOW as reference; S33, Establishing the database of isotope ratio of known pollution sources: manure sewage δ 15 N=7-15‰, δ 18 O=5-12‰, chemical fertilizer δ 15 N=-5-5‰, δ 18 O=-10-5‰, atmospheric deposition δ 15 N=-10-0‰, δ 18 O=10-20‰, soil organic nitrogen δ 15 N=2-8‰, δ 18 O=8-15‰, Comparing the isotope ratio of sample with the database to determine the main source of nitrogen pollution; The real-time monitoring chip for denitrifying functional flora activity adopts a silicon-based nano electrode array, takes the activity of nitratereductase enzyme as the detection index, the detection limit is 0.01 U / L, the monitoring frequency is 1 day before autumn irrigation, 1 day, 3 days and 7 days after autumn irrigation, 3 parallel samples are collected each time, and the activity data are subjected to trend analysis by using Origin software.

6. The method according to claim 2, wherein the method is characterized in that, The formula for calculating the Nemerow Comprehensive Pollution Index in S4 is: PI = √[(PI avg 2 +PI max 2 ) / 2], where PI is the Nemerow Composite Pollution Index, PI avg PI is the arithmetic mean of the pollution indices of each pollutant. max The pollution index P represents the maximum value among all pollution indices. Ii =C i / S i C i S represents the measured concentration of the pollutant. i The concentration is the standard concentration for pollutant evaluation; the formula for calculating the centroid coordinates of the centroid migration model is: X = ∑(Q i ×X i ) / ∑Q i Y=∑(Q i ×Y i ) / ∑Q i Where X and Y are the abscissa and ordinate of the pollution load centroid, respectively, and Qi is the pollution load of the i-th evaluation unit. i Y i These are the x-coordinate and y-coordinate of the geometric center of the i-th evaluation unit, respectively. The evaluation units are divided into 1km×1km grids, with a total of 500 evaluation units. The correlation model adopts Pearson linear correlation analysis, and the correlation coefficients R of irrigation amount, fertilizer amount, rainfall and pollution load are calculated by SPSS software, wherein the correlation coefficient R of irrigation amount and pollution load is 0.85, the correlation coefficient R of fertilizer amount is 0.82, and the correlation coefficient R of rainfall is 0.74, all of which pass the significance test of P<0.

01.

7. The method according to claim 1, wherein the method is characterized in that, The construction steps of the "irrigation return water-soil seepage-surface runoff" coupling migration model in S5 are as follows: S51, a soil seepage sub-model is constructed based on the HYDRUS-1D model, the soil data collected in S1 is input, and the soil water seepage rate under different irrigation amounts is simulated; S52, a surface runoff sub-model is constructed based on the SWAT model, the DEM data, land use data (70% of cultivated land, 20% of grassland and 10% of water area), and meteorological data collected in S1 are input, and the surface runoff flow is calculated by using the SCS curve number method; S53, an irrigation return water sub-model is constructed, the irrigation water amount, ditch slope and soil infiltration rate are input, and the irrigation return water flow is calculated by using the water balance equation Qreturn water=Qirrigation-Qinfiltration-Qevaporation; S54, a coupling interface is developed by using C++ language, the output data of the soil seepage sub-model, the surface runoff sub-model and the irrigation return water sub-model are coupled, the migration coefficients of pollutants in seepage water, surface runoff and irrigation return water are set as 0.8, 0.9 and 0.95 respectively, the migration path and concentration change of pollutants in the water-soil system are simulated, and the measured concentration data collected in S1 is used for model verification, and the verification error is controlled within 15%.

8. The method according to claim 1, wherein the method is characterized in that, The R-type cluster analysis method in S6 adopts the Ward clustering method, which is realized by using SPSS 26.0 software, and the specific steps are as follows: S61. Standardize the concentration data of the nine pollutants obtained in S1-S5 (standardization formula: x'=(xx)). min ) / (x max -x min (Eliminating the influence of dimensions). S62, calculate the Euclidean distance between the standardized variables (distance formula d ij =√[∑(x ik -x jk )²] ), construct a distance matrix; S63, the Ward method is used for clustering, that is, the increase of the within-class sum of squares is minimized when two classes are merged each time, the clustering threshold is set as the variance contribution rate being greater than 85%, the variance contribution rate of each principal component is calculated by principal component analysis, and the clustering is stopped when the cumulative variance contribution rate reaches 85%; S64, the validity of the clustering results is verified by using the silhouette coefficient, and a calculation formula of the silhouette coefficient is s i =(b i -a i ) / max(a i ,b i ), wherein a i is an average distance of the ith sample and other samples of the same type, b i is an average distance of the ith sample and the nearest sample of the different type, and when the average silhouette coefficient is greater than 0.7, it is determined that the clustering result is valid, and finally, the four pollution load clustering types are divided as follows: The reduced nitrogen cluster, the total phosphorus-TDS-chemical oxygen demand cluster, the perfluorinated pollutant cluster and the total nitrogen cluster.

9. The method according to claim 1, wherein the method is characterized in that, The parameter optimization of the Invest model in S7 includes: The field water capacity and saturated hydraulic conductivity in the soil data are adjusted to the measured values based on the in-situ detection in S1, the growth period in the crop parameters is adjusted to the actual growth period of sunflowers in the Yellow River irrigation area, and the crop coefficients are adjusted as follows: 0.3 in the emergence stage, 0.7 in the jointing stage, 1.1 in the present budding stage, 1.0 in the grain filling stage and 0.6 in the mature stage; The "topographic water collection-irrigation and fertilization superposition driving module" is a self-developed module, which is written by using Python language, and the input parameters include topographic slope, irrigation amount, nitrogen application amount and monthly average rainfall, and the module calculates the pollution loss driving intensity under different parameter combinations through superposition analysis; The soil column experiment uses a transparent organic glass column filled with undisturbed soil, sets up three treatment groups with three repetitions in each group, collects soil solutions at different depths through a sampling hole, detects nitrogen and phosphorus concentrations, calculates nitrogen and phosphorus migration rates, and determines that the nitrogen mobility is significantly higher than that of phosphorus.

10. The method according to claim 1, wherein the method is characterized in that, The specific steps of the self-organizing neural network combined with the positive definite matrix factor analysis method in S8 are as follows: S81, the chemical component data of the canal water and the drainage sample are subjected to dimension reduction processing by using the self-organizing neural network, the SOM network is set to an output layer of 4*4 and an input layer of 11 neurons, the initial learning rate is 0.05, the iteration number is 1000, and six characteristic vectors are obtained after dimension reduction; S82, the six characteristic vectors are input into the positive definite matrix factor analysis model, the number of factors is set to 3, and the characteristic components of each source are determined by rotating the factor loading matrix: the characteristic components of the evaporation rock dissolution are TDS, Cl- and SO42-; the characteristic components of the silicate weathering are Ca2+, Mg2+ and HCO3-; and the characteristic components of the agricultural input are total nitrogen, total phosphorus and fluoride; S83, the contribution proportions of each source are calculated: the contribution proportion of the evaporation rock dissolution is 35-45%, the contribution proportion of the silicate weathering is 20-30%, and the contribution proportion of the agricultural input is 30-40%, the error of the contribution proportion is calculated by the bootstrap method, and the error is controlled within 5%; and the specific steps of the entropy weight water quality evaluation method are as follows: S84, the water quality index is standardized, the positive index standardization formula x'= (x-x min ) / (x max -x min ), the negative index standardization formula x'= (x max -x) / (x max -x min ) S85、Calculate the entropy value e of each index j = -k x ∑(p ij x ln p ij ), where k = 1 / ln m (m is the number of samples, m = 100 in this embodiment), p ij = x' ij / ∑x' ij (if x' ij = 0, let p ij = 1 / m); S86, calculate the weight w of each index j =(1-e j ) / ∑(1-e j ) S87、Calculate water quality comprehensive index FI =∑(w j x F ij ), wherein F ij is the score value of the jth index of the ith sample, when FI < 0.6, it is determined that it is suitable for drinking, and when FI≥0.6, it is determined that it is not suitable for drinking.

Images