A non-point source pollution intensity evaluation method based on multi-source data fusion

Abstract

The application discloses a kind of based on multi-source data fusion's non-point source pollution intensity evaluation method, comprising the following steps: S1.Through air space ground coordination network acquisition target watershed multi-source environmental data;S2.Coupling model is constructed to environmental resistance-pollution source intensity, and pollution migration index and dynamic load are calculated;S3.Non-point source pollution intensity index is calculated, and non-point source pollution risk classification is carried out.The application is through air space ground multi-source data cooperation, utilizes space-based satellite remote sensing system, air-based unmanned aerial vehicle monitoring system, ground-based internet of things sensing system acquisition target watershed environmental data, for the non-point source pollution intensity calculation and non-point source pollution risk evaluation of target watershed.

Description

Technical Field

[0001] This invention relates to the field of non-point source pollution assessment technology, specifically to a method for evaluating the intensity of non-point source pollution based on multi-source data fusion. Background Technology

[0002] Non-point source pollution (NPP) presents significant challenges to accurate monitoring and intensity quantification in the environmental field due to its dispersed, random, and time-dependent characteristics. Current mainstream technologies suffer from three major bottlenecks: First, fragmented data acquisition: traditional ground monitoring stations are sparse, making it difficult to cover complex underlying surface areas (such as farmland and mountains), resulting in widespread pollution "blind spots"; manual sampling has long cycles and cannot capture the dynamics of sudden pollution events such as rainfall and runoff. Second, insufficient model evaluation accuracy: existing models (such as the output coefficient method) rely on static parameters and do not consider the spatial heterogeneity of environmental resistance (such as surface roughness, slope, and vegetation cover) and its hindering effect on pollutant migration, leading to biased load estimations. Third, lack of multi-source data fusion: satellite, UAV, and ground sensor data are processed independently, lacking spatiotemporal alignment and coupling analysis, making it difficult to construct a full-chain dynamic assessment system encompassing "pollution source-migration path-receptor." Therefore, in the current process of evaluating the intensity of non-point source pollution, although there are ways to build corresponding evaluation models based on different scales, the parameters used in the existing technology are often very limited. For example, relying on sensors can only obtain indicators such as total phosphorus and total nitrogen in water bodies, while relying on drones can obtain watershed image data. However, these data from different sources have significant differences in both dimensions and spatiotemporal scales. These parameters will have different degrees of impact on non-point source pollution. The existing technology lacks coupling of these multi-source heterogeneous data, so it is impossible to effectively obtain more scientific evaluation results of the intensity of non-point source pollution.

[0003] To address these issues, this invention proposes a method for evaluating area source pollution based on multi-source data fusion technology from air, space, and ground, in order to solve the aforementioned problems. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention aims to provide a high-precision, real-time, and adaptive method for evaluating the intensity of non-point source pollution. By constructing a space-air-ground collaborative monitoring network, it achieves full coverage data collection across multiple scales (watershed-plot) and dimensions (water quality-soil-meteorology). It also develops a multi-parameter dynamic coupling model and optimizes the accuracy of pollution load quantification through environmental resistance correction and machine learning.

[0005] A method for assessing the intensity of non-point source pollution based on multi-source data fusion, characterized in that the assessment method includes the following steps:

[0006] S1. Collect multi-source environmental data of the target watershed through a space-air-ground integrated network;

[0007] In step S1, the space-air-ground collaborative network includes: a space-based satellite remote sensing system, an air-based unmanned aerial vehicle (UAV) monitoring system, and a ground-based Internet of Things (IoT) sensing system;

[0008] S2. Construct an environmental resistance-pollution source coupling model to calculate the pollution migration index and dynamic load;

[0009] S3. Calculate the non-point source pollution intensity index and classify the non-point source pollution risk.

[0010] The multi-source environmental data includes watershed land use type, vegetation cover, surface parameters, topographic data, soil physicochemical indicators, water body physicochemical indicators, soil permeability, and real-time rainfall intensity; the space-based satellite remote sensing system is used to acquire watershed land use type and vegetation cover, the air-based UAV monitoring system is used to collect surface parameters and topographic data, and the ground-based Internet of Things sensing system is used to monitor soil physicochemical indicators, water body physicochemical indicators, soil permeability, and real-time rainfall intensity in real time.

[0011] The airborne UAV monitoring system is equipped with multispectral sensors and lidar.

[0012] The surface parameters include surface roughness and water flow path distance.

[0013] The soil physicochemical indicators include total phosphorus, total nitrogen, and soil organic matter content.

[0014] The physicochemical indicators of the water body include total phosphorus, total nitrogen, dissolved oxygen, and turbidity.

[0015] The specific steps involved in obtaining the watershed land use type are as follows:

[0016] S11. Select a satellite data source and preprocess the data, including radiometric calibration, atmospheric correction, geometric fine correction, image fusion, image mosaicking and cropping;

[0017] S12. Perform feature extraction, the features including spectral features, vegetation index, water index, building index, impermeable surface index, texture features, spatial features, morphological features, time series features, terrain features, and contextual features;

[0018] S13. Classify land use types. Based on the 3D-CNN model, divide the target watershed into four categories: agricultural area, forestry area, industrial area and residential area, and output the area and proportion of each type of land.

[0019] Step S2 includes the following specific steps:

[0020] S21. Establish a unified spatiotemporal coordinate system for multi-source environmental data. The spatiotemporal coordinate system is established as follows:

[0021]

[0022] In the formula: —A set of spatiotemporal coordinate systems for the target watershed;

[0023] — Column coordinates of the grid, ranging from 1 to ;

[0024] — The row coordinates of the grid, ranging from 1 to ;

[0025] —The total number of grid cells in the x-direction;

[0026] —The total number of grid cells in the y direction;

[0027] —Time variable;

[0028] —Time window range, defined as ,in For the start time, End time;

[0029] Dimensionless processing was performed on multi-source environmental data. The formula for dimensionless processing of vegetation cover is as follows:

[0030]

[0031] In the formula: —Grid Normalized vegetation cover, with values ​​ranging from 0 to 1;

[0032] —Grid Original vegetation cover;

[0033] —Minimum vegetation coverage for the target watershed;

[0034] —Maximum vegetation coverage in the target watershed;

[0035] The standardized formula for calculating surface roughness is as follows:

[0036]

[0037] In the formula: —Grid Standardized surface roughness value;

[0038] —Grid The original surface roughness value;

[0039] —Minimum surface roughness of the watershed;

[0040] —Maximum surface roughness of the watershed;

[0041] The formula for calculating normalized soil permeability is as follows:

[0042]

[0043] In the formula: —Grid Normalized soil permeability;

[0044] —Grid Original soil permeability;

[0045] —Minimum soil permeability in the watershed;

[0046] —Maximum soil permeability in the watershed;

[0047] The formula for dimensionless treatment of slope is as follows:

[0048]

[0049] In the formula: —Grid Standardized slope value;

[0050] —Grid Original slope value, in units: ;

[0051] —Maximum slope value, in units of: ;

[0052] S22. Based on the land use types obtained in step S1, construct the pollution equivalent matrix C.

[0053]

[0054] In the formula: —The pollution equivalent coefficient for agricultural areas is 0.85;

[0055] —The pollution equivalent coefficient for forestry areas is 0.15;

[0056] —The pollution equivalent coefficient for the industrial zone is 0.6;

[0057] —The pollution equivalent coefficient for residential areas is 0.7;

[0058] The formula for calculating pollution source intensity is as follows:

[0059]

[0060] In the formula: —Grid The intensity of the pollution source at time t;

[0061] —Grid The area proportion of land of type k, where k=1 is agricultural area, k=2 is forestry area, k=3 is industrial area, and k=4 is residential area;

[0062] —Pollution equivalent coefficient, k=1, 2, 3, 4;

[0063] —Yearly accumulated days, valued as 1-365;

[0064] —Seasonal offset;

[0065] —Seasonal factor intensity coefficient, taken as 0.5;

[0066] —Grid The intensity of rainfall;

[0067] —The maximum rainfall intensity in the region;

[0068] —Rainfall erosivity index, taken as 1.2;

[0069] S23. Construction of the environmental resistance model, including the construction of the resistance factor matrix and spatial principal component analysis, wherein the resistance factor matrix is:

[0070]

[0071] In the formula: —Grid The drag factor vector;

[0072] The spatial principal component analysis includes constructing a local neighborhood data matrix for the central grid point. Take it Neighborhood:

[0073] ;

[0074] in, No more than , The smaller value in ;

[0075] Calculate the local covariance matrix:

[0076]

[0077] In the formula: —Local covariance matrix;

[0078] —Neighborhood mean vector ,in For the first Feature vectors of neighboring points;

[0079] Eigenvalue decomposition is performed on the local covariance matrix to obtain eigenvalues ​​and eigenvectors;

[0080] Perform drag coefficient synthesis:

[0081]

[0082] In the formula: —Grid Environmental resistance coefficient;

[0083] —The weight of the j-th principal component, , Local covariance matrix The j-th eigenvalue is given, and the eigenvalues ​​are arranged in ascending order. Local covariance matrix The sum of all eigenvalues;

[0084] —The eigenvector corresponding to the j-th eigenvalue of the local covariance matrix;

[0085] S24. Calculation of Pollution Migration Index:

[0086]

[0087] In the formula: —The pollution migration index of grid (x,y) at time t;

[0088] —The set of upstream grids flowing towards grid (x,y);

[0089] —The pollution source intensity of the upstream grid (u,v) at time t. ,in For grid The area proportion of land of type k, where k=1 is agricultural area, k=2 is forestry area, k=3 is industrial area, and k=4 is residential area;

[0090] —Drag influence coefficient, taken as 0.8;

[0091] —Distance attenuation coefficient, taken as 0.2;

[0092] —The distance of the water flow path from grid (u,v) to (x,y), in kilometers;

[0093] S25. Perform dynamic load forecasting, wherein the dynamic load forecasting uses the following formula:

[0094]

[0095] In the formula: —Grid(x,y) in future time Pollution load;

[0096] —Migration conversion factor, ranging from 0.8 to 1.2;

[0097] —The pollution migration index of grid (u,v) at time t;

[0098] —Base weight, ranging from 0.1 to 0.3;

[0099] —The environmental background of the grid (x,y);

[0100] —Dynamic correction item;

[0101] The formula for calculating the dynamic correction term is:

[0102]

[0103] In the formula: — The calibration parameter is 0.05;

[0104] —Real-time rainfall intensity, mm / h;

[0105] —Baseline rainfall, which is the monthly average;

[0106] —The pollution migration index of grid (x,y) at time t;

[0107] The formula for calculating the environmental background of the grid (x,y) is as follows:

[0108]

[0109] In the formula: —Soil weight, ranging from 0.3 to 0.7;

[0110] —Water body weight, taken as 0.2-0.8;

[0111]

[0112] In the formula: —Soil total phosphorus weighting coefficient;

[0113] —Soil total nitrogen weighting coefficient;

[0114] —Soil organic matter content weighting coefficient ;

[0115] —Total phosphorus in soil;

[0116] —Total nitrogen in soil;

[0117] —Soil organic matter content;

[0118]

[0119] In the formula: —Total phosphorus weighting coefficient in water bodies;

[0120] —Total nitrogen weighting coefficient in water bodies;

[0121] —Dissolved oxygen weighting coefficient in water body;

[0122] —Water turbidity weighting coefficient + + + =1;

[0123] —Total phosphorus in water;

[0124] —Total nitrogen in water;

[0125] —Dissolved oxygen in water;

[0126] —Water turbidity;

[0127] In step S3, the formula for calculating the intensity of non-point source pollution is:

[0128]

[0129] In the formula: —The area source pollution intensity index of the grid (x,y);

[0130] —Historical minimum pollution load value;

[0131] —The highest historical pollution load value.

[0132] In step S3, the specific method for classifying the risk of non-point source pollution is as follows: when The risk level of non-point source pollution is low; when The risk level of non-point source pollution is medium; when The risk level of non-point source pollution is high.

[0133] Compared with the prior art, the beneficial effects of the present invention are:

[0134] This invention utilizes multi-source data collaboration across air, space, and ground. It employs space-based satellite remote sensing systems, airborne unmanned aerial vehicle (UAV) monitoring systems, and ground-based Internet of Things (IoT) sensing systems to collect environmental data from target watersheds. This data is then used to calculate the intensity of non-point source pollution (NPP) in the target watersheds, achieving comprehensive utilization of various types of data related to NPP, including remote sensing data, point cloud data, and sensor data. Through the construction and data acquisition of a three-dimensional air-space-ground network, it is applicable to NPP assessment in different watersheds nationwide.

[0135] Data fusion maps satellite remote sensing data, UAV point cloud data, and ground sensor data onto a unified watershed geographic grid, eliminating scale differences and temporal shifts. Simultaneously, it performs dimensionless processing on multi-source heterogeneous data, achieving environmental data unification and overcoming the technical bottleneck of coupling multi-source heterogeneous parameters. This allows a series of parameters affecting non-point source pollution, such as land use type, vegetation cover, soil physicochemical indicators, water body physicochemical indicators, and real-time rainfall intensity, to participate in the assessment of non-point source pollution intensity.

[0136] The watershed land use types were taken into account, as different land use types contribute differently to non-point source pollution. The pollution source intensity was calculated by combining the land use types, which ensured the accuracy of the pollution source intensity calculation results.

[0137] In terms of coupled model construction, a calculation model for pollution source intensity was first constructed based on parameters such as land type and rainfall intensity. Furthermore, an environmental resistance model was constructed based on parameters such as vegetation cover, surface roughness, soil permeability, and slope. The environmental resistance coefficient of the grid was calculated based on the environmental resistance model. Then, the pollution migration index was calculated based on the pollution source intensity and water flow path distance. Finally, a dynamic pollution load calculation model was constructed by coupling parameters such as pollution source intensity, environmental resistance, environmental background, and rainfall intensity. This model fully considers the impact of migration costs, environmental background, and rainfall on the pollution load, achieving accurate calculation of the dynamic pollution load. The non-point source pollution intensity was obtained based on the pollution load. Attached Figure Description

[0138] Figure 1 This is a flowchart illustrating the overall evaluation method of the present invention.

[0139] Figure 2 This is a flowchart for obtaining watershed land use types according to the present invention. Detailed Implementation

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

[0141] This invention provides a method for evaluating the intensity of area source pollution based on multi-source data fusion. The evaluation method includes the following steps:

[0142] S1. Collect multi-source environmental data of the target watershed through a space-air-ground integrated network;

[0143] In step S1, the space-air-ground collaborative network includes: a space-based satellite remote sensing system, an air-based unmanned aerial vehicle (UAV) monitoring system, and a ground-based Internet of Things (IoT) sensing system;

[0144] S2. Construct an environmental resistance-pollution source coupling model to calculate the pollution migration index and dynamic load;

[0145] S3. Calculate the non-point source pollution intensity index and classify the non-point source pollution risk.

[0146] The multi-source environmental data includes watershed land use type, vegetation cover, surface parameters, topographic data, soil physicochemical indicators, and water body physicochemical indicators; the space-based satellite remote sensing system is used to acquire watershed land use type and vegetation cover, the airborne UAV monitoring system is used to collect surface parameters and topographic data, and the ground-based Internet of Things sensing system is used to monitor soil physicochemical indicators, water body physicochemical indicators, soil permeability, and real-time rainfall intensity in real time.

[0147] The airborne UAV monitoring system is equipped with multispectral sensors and lidar.

[0148] The surface parameters include surface roughness and water flow path distance.

[0149] The soil physicochemical indicators include total phosphorus, total nitrogen, and soil organic matter content.

[0150] The physicochemical indicators of the water body include total phosphorus, total nitrogen, dissolved oxygen, and turbidity.

[0151] Because different land features possess unique remote sensing characteristics, such as the high reflectivity of vegetation (crops) in agricultural areas in the near-infrared band and low reflectivity in the red band, the vegetation index in agricultural and forestry areas exhibits strong periodic changes with the seasons (sowing-growth-maturation-harvest). Industrial areas, on the other hand, are typically clustered in large areas with clear boundaries, and factories often exhibit distinct high-temperature zones in the thermal infrared band. Residential areas display highly structured textures, exhibiting dense, fragmented grid-like or honeycomb patterns, and their temporal characteristics change relatively slowly. Forestry areas are spatially concentrated in non-flat areas such as mountains and hills, and are strongly correlated with topography.

[0152] The specific steps involved in obtaining the watershed land use type are as follows:

[0153] S11. Select a satellite data source and preprocess the data, including radiometric calibration, atmospheric correction, geometric fine correction, image fusion, image mosaicking and cropping;

[0154] Radiometric calibration converts the original brightness values ​​of pixels into radiance or apparent reflectance. Atmospheric correction eliminates atmospheric effects to obtain the true surface reflectance. Geometric fine correction ensures accurate spatial location. Image fusion improves spatial resolution. Image mosaicking and cropping cover the target watershed area.

[0155] S12. Perform feature extraction, the features including spectral features, vegetation index, water index, building index, impermeable surface index, texture features, spatial features, morphological features, time series features, terrain features, and contextual features;

[0156] The spectral features include the original band reflectance values; vegetation indices include the Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EDI), and Soil-Adjusted Vegetation Index (SDI); water indices include the Normalized Difference Water Index (NDDI) to identify rivers, lakes, and irrigation ditches; building indices and impermeable surface indices are used to identify man-made building areas; texture features are used to distinguish between residential areas (high texture complexity) and farmland (low texture complexity); spatial and morphological features include plot aspect ratio, compactness, boundary density, and road density; time series features are used to calculate the vegetation index time series curve for the entire growing season (or many years), using Sentinel 2 or Landsat dense time series; terrain features include slope and aspect; contextual features mainly consider the category of pixels surrounding the target pixel, such as regular rectangular plots near rivers are likely farmland, while large building areas near main roads and railways may be industrial.

[0157] S13. Classify land use types. Based on the 3D-CNN model, divide the target watershed into four categories: agricultural area, forestry area, industrial area and residential area, and output the area and proportion of each type of land.

[0158] The 3D-CNN model architecture includes a dual-branch input processing. The spatiotemporal feature branch uses a 3×5×5 three-dimensional convolutional kernel (time×height×width) to capture seasonal changes and spatial texture. The first layer extracts macroscopic phenological features (such as crop growth cycles), retaining the temporal dimension through 1×2×2 pooling while compressing only the spatial dimension. The spectral feature branch uses a 1×1×1 convolutional kernel for intelligent band weighting to adaptively optimize multispectral combinations (such as enhancing the spectral differences between vegetation and buildings). Feature fusion concatenates the spatiotemporal feature map (32 channels) and the spectral feature map (32 channels) along the channel dimension to form a 64-channel fused feature. Deep feature extraction uses three levels of 3×3×3 convolutional layers (channel count 64→128→256), progressively performing 2×2×2 spatiotemporal pooling with a compression rate of 50%. Each layer includes batch normalization (BN) and ReLU activation. In terms of classification mechanism, the 256-channel time × height × width feature map is compressed into a 256-dimensional vector, and four land use types—agricultural area, forestry area, industrial area, and residential area—are output through global average pooling and fully connected layers.

[0159] Step S2 includes the following specific steps:

[0160] S21. Establish a unified spatiotemporal coordinate system for multi-source environmental data. The spatiotemporal coordinate system is established as follows:

[0161]

[0162] In the formula: —A set of spatiotemporal coordinate systems for the target watershed;

[0163] — Column coordinates of the grid, ranging from 1 to ;

[0164] — The row coordinates of the grid, ranging from 1 to ;

[0165] —The total number of grid cells in the x-direction;

[0166] —The total number of grid cells in the y direction;

[0167] —Time variable;

[0168] —Time window range, defined as ,in For the start time, End time;

[0169] Dimensionless processing was performed on multi-source environmental data. The formula for dimensionless processing of vegetation cover is as follows:

[0170]

[0171] In the formula: —Grid Normalized vegetation cover, with values ​​ranging from 0 to 1;

[0172] —Grid Original vegetation cover;

[0173] —Minimum vegetation coverage for the target watershed;

[0174] —Maximum vegetation coverage in the target watershed;

[0175] The standardized formula for calculating surface roughness is as follows:

[0176]

[0177] In the formula: —Grid Standardized surface roughness value;

[0178] —Grid The original surface roughness value;

[0179] —Minimum surface roughness of the watershed;

[0180] —Maximum surface roughness of the watershed;

[0181] The formula for calculating normalized soil permeability is as follows:

[0182]

[0183] In the formula: —Grid Normalized soil permeability;

[0184] —Grid Original soil permeability;

[0185] —Minimum soil permeability in the watershed;

[0186] —Maximum soil permeability in the watershed;

[0187] The formula for dimensionless treatment of slope is as follows:

[0188]

[0189] In the formula: —Grid Standardized slope value;

[0190] —Grid Original slope value, in units: ;

[0191] —Maximum slope value, in units of: ;

[0192] S22. Based on the land use types obtained in step S1, construct the pollution equivalent matrix C.

[0193]

[0194] In the formula: —The pollution equivalent coefficient for agricultural areas is 0.85;

[0195] —The pollution equivalent coefficient for forestry areas is 0.15;

[0196] —The pollution equivalent coefficient for the industrial zone is 0.6;

[0197] —The pollution equivalent coefficient for residential areas is 0.7;

[0198] Pollution emission factors (emissions per unit area) for each land use type may vary due to geographical location and actual environment. In this embodiment, a pollution equivalent matrix is ​​constructed based on the main non-point source pollutants of different land types. This indicates that agricultural areas contribute 85% to total nitrogen and total phosphorus emissions. This indicates that the forestry area is concerned about NO x Its contribution to emissions is 15%. This indicates that industrial zones contribute 60% to the emissions of heavy metals and COD (chemical oxygen demand). This indicates that residential areas contribute 70% to organic matter emissions. Other pollutants, such as carbon dioxide and PM2.5, also exist in agricultural, forestry, industrial, and residential areas, but they are not the main non-point source pollutants, so the pollution equivalent matrix in this embodiment does not consider such pollutants.

[0199] The formula for calculating pollution source intensity is as follows:

[0200]

[0201] In the formula: —Grid The intensity of the pollution source at time t;

[0202] —Grid The area proportion of land of type k, where k=1 is agricultural area, k=2 is forestry area, k=3 is industrial area, and k=4 is residential area;

[0203] —Pollution equivalent coefficient, k=1, 2, 3, 4;

[0204] —Yearly accumulated days, valued as 1-365;

[0205] —Seasonal offset;

[0206] —Seasonal factor intensity coefficient, taken as 0.5;

[0207] —Grid The intensity of rainfall;

[0208] —The maximum rainfall intensity in the region;

[0209] —Rainfall erosivity index, taken as 1.2;

[0210] S23. Construction of the environmental resistance model, including the construction of the resistance factor matrix and spatial principal component analysis, wherein the resistance factor matrix is:

[0211]

[0212] In the formula: —Grid The drag factor vector;

[0213] The spatial principal component analysis includes constructing a local neighborhood data matrix for the central grid point. Take it Neighborhood:

[0214] ;

[0215] in, No more than , The smaller value in ;

[0216] Calculate the local covariance matrix:

[0217]

[0218] In the formula: —Local covariance matrix;

[0219] —Neighborhood mean vector ,in For the first Feature vectors of neighboring points;

[0220] Eigenvalue decomposition is performed on the local covariance matrix to obtain eigenvalues ​​and eigenvectors;

[0221] Perform drag coefficient synthesis:

[0222]

[0223] In the formula: —Grid Environmental resistance coefficient;

[0224] —The weight of the j-th principal component, , Local covariance matrix The j-th eigenvalue is given, and the eigenvalues ​​are arranged in ascending order. Local covariance matrix The sum of all eigenvalues;

[0225] —The eigenvector corresponding to the j-th eigenvalue of the local covariance matrix;

[0226] S24. Calculation of Pollution Migration Index:

[0227]

[0228] In the formula: —The pollution migration index of grid (x,y) at time t;

[0229] —The set of upstream grids flowing towards grid (x,y);

[0230] —The pollution source intensity of the upstream grid (u,v) at time t. ,in For grid The area proportion of land of type k, where k=1 is agricultural area, k=2 is forestry area, k=3 is industrial area, and k=4 is residential area;

[0231] —Drag influence coefficient, taken as 0.8;

[0232] —Distance attenuation coefficient, taken as 0.2;

[0233] —The distance of the water flow path from grid (u,v) to (x,y), in kilometers;

[0234] S25. Perform dynamic load forecasting, wherein the dynamic load forecasting uses the following formula:

[0235]

[0236] In the formula: —Grid(x,y) in future time Pollution load;

[0237] —Migration conversion factor, ranging from 0.8 to 1.2;

[0238] —The pollution migration index of grid (u,v) at time t;

[0239] —Base weight, ranging from 0.1 to 0.3;

[0240] —The environmental background of the grid (x,y);

[0241] —Dynamic correction item;

[0242] The formula for calculating the dynamic correction term is:

[0243]

[0244] In the formula: — The calibration parameter is 0.05;

[0245] —Real-time rainfall intensity, mm / h;

[0246] —Baseline rainfall, which is the monthly average;

[0247] —The pollution migration index of grid (x,y) at time t;

[0248] The formula for calculating the environmental background of the grid (x,y) is as follows:

[0249]

[0250] In the formula: —Soil weight, ranging from 0.3 to 0.7;

[0251] —Water body weight, taken as 0.2-0.8;

[0252]

[0253] In the formula: —Soil total phosphorus weighting coefficient;

[0254] —Soil total nitrogen weighting coefficient;

[0255] —Soil organic matter content weighting coefficient ;

[0256] —Total phosphorus in soil;

[0257] —Total nitrogen in soil;

[0258] —Soil organic matter content;

[0259]

[0260] In the formula: —Total phosphorus weighting coefficient in water bodies;

[0261] —Total nitrogen weighting coefficient in water bodies;

[0262] —Dissolved oxygen weighting coefficient in water body;

[0263] —Water turbidity weighting coefficient + + + =1;

[0264] —Total phosphorus in water;

[0265] —Total nitrogen in water;

[0266] —Dissolved oxygen in water;

[0267] —Water turbidity;

[0268] In step S3, the formula for calculating the intensity of non-point source pollution is:

[0269]

[0270] In the formula: —The area source pollution intensity index of the grid (x,y);

[0271] —Historical minimum pollution load value;

[0272] —The highest historical pollution load value.

[0273] In step S3, the specific method for classifying the risk of non-point source pollution is as follows: when The risk level of non-point source pollution is low; when The risk level of non-point source pollution is medium; when The risk level of non-point source pollution is high.

[0274] Following the above method, taking a small watershed of a tributary of River A in a certain city as an example, a 1km×1km grid was used for modeling, and the downstream grid was selected for calculation. Within this grid, the agricultural area accounts for 50%, the forestry area for 20%, the industrial area for 10%, and the residential area for 20%. The rainfall intensity of this grid was calculated. Maximum rainfall intensity in the region Years accumulated =150, Seasonal Offset The two upstream grid cells of this grid and The relevant information is shown in Table 1.

[0275] Table 1 Upstream Grid and Calculation results

[0276] ;

[0277] The relevant parameters of the environmental background are as follows:

[0278] , , , , , ;

[0279] , , , , , , , .

[0280] Soil weight Water weight ;

[0281] Baseline rainfall ;

[0282] migration conversion coefficient ;

[0283] Base weight ;

[0284] Historical minimum pollution load value ;

[0285] Historical maximum pollution load value ;

[0286] but

[0287] ;

[0288] ;

[0289] ;

[0290] Environmental background 0.59 + 0.5 0.45 = 0.52

[0291] Dynamic correction item ;

[0292]

[0293]

[0294] Therefore, the risk level of the grid area source pollution is classified as low risk.

[0295] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the invention can be implemented in other specific forms without departing from its spirit or essential characteristics. Although this specification describes embodiments, not every embodiment contains only one technical solution. This method of description is merely for clarity. Those skilled in the art should consider the specification as a whole, and the technical solutions in each embodiment can be appropriately combined to form other embodiments that can be understood by those skilled in the art.

Claims

1. A method for evaluating the intensity of non-point source pollution based on multi-source data fusion, characterized in that, The evaluation method includes the following steps: S1. Collecting multi-source environmental data of the target watershed through an air-space-ground collaborative network; in step S1, the air-space-ground collaborative network includes: a space-based satellite remote sensing system, an air-based unmanned aerial vehicle monitoring system, and a ground-based Internet of Things sensing system; S2. Constructing an environmental resistance-pollution source coupling model and calculating the pollution migration index and dynamic load; S3. Calculating the non-point source pollution intensity index and classifying the non-point source pollution risk. In step S2, the pollution migration index is calculated as follows: In the formula: —The pollution migration index of grid (x,y) at time t; —The set of upstream grids flowing towards grid (x,y); —The pollution source intensity of the upstream grid (u,v) at time t. ,in For grid The area proportion of land of type k, where k=1 is agricultural area, k=2 is forestry area, k=3 is industrial area, and k=4 is residential area; —Drag influence coefficient, taken as 0.8; —Grid Environmental resistance coefficient; —Distance attenuation coefficient, taken as 0.2; —The water flow path distance from grid (u,v) to (x,y), in kilometers; step S2 also includes dynamic load prediction, which uses the following formula: In the formula: —Grid(x,y) in future time Pollution load; —Migration conversion factor, ranging from 0.8 to 1.2; —The pollution migration index of grid (u,v) at time t; —Base weight, ranging from 0.1 to 0.3; —The environmental background of the grid (x,y); —Dynamic correction term; the calculation formula for the dynamic correction term is: In the formula: — The calibration parameter is 0.05; —Real-time rainfall intensity, mm / h; —Baseline rainfall, which is the monthly average; —The pollution migration index of grid (x,y) at time t; the environmental background of grid (x,y) is calculated using the following formula: In the formula: —Soil weight, ranging from 0.3 to 0.7; —Water body weight, taken as 0.2-0.8; In the formula: —Soil total phosphorus weighting coefficient; —Soil total nitrogen weighting coefficient; —Soil organic matter content weighting coefficient ; —Total phosphorus in soil; —Total nitrogen in soil; —Soil organic matter content; In the formula: —Total phosphorus weighting coefficient in water bodies; —Total nitrogen weighting coefficient in water bodies; —Dissolved oxygen weighting coefficient in water body; —Water turbidity weighting coefficient + + + =1; —Total phosphorus in water; —Total nitrogen in water; —Dissolved oxygen in water; —Water turbidity; In step S3, the formula for calculating the intensity of non-point source pollution is: In the formula: —The area source pollution intensity index of the grid (x,y); —Historical minimum pollution load value; —The highest historical pollution load value.

2. The method for evaluating the intensity of area source pollution based on multi-source data fusion according to claim 1, characterized in that, The multi-source environmental data includes watershed land use type, vegetation cover, surface parameters, topographic data, soil physicochemical indicators, water body physicochemical indicators, soil permeability, and real-time rainfall intensity; the space-based satellite remote sensing system is used to acquire watershed land use type and vegetation cover, the air-based UAV monitoring system is used to collect surface parameters and topographic data, and the ground-based Internet of Things sensing system is used to monitor soil physicochemical indicators, water body physicochemical indicators, soil permeability, and real-time rainfall intensity in real time.

3. The method for evaluating the intensity of area source pollution based on multi-source data fusion according to claim 1, characterized in that, The airborne UAV monitoring system is equipped with multispectral sensors and lidar.

4. The method for evaluating the intensity of area source pollution based on multi-source data fusion according to claim 2, characterized in that, The surface parameters include surface roughness and water flow path distance.

5. The method for evaluating the intensity of area source pollution based on multi-source data fusion according to claim 4, characterized in that, The soil physicochemical indicators include total phosphorus, total nitrogen, and soil organic matter content.

6. The method for evaluating the intensity of area source pollution based on multi-source data fusion according to claim 5, characterized in that, The physicochemical indicators of the water body include total phosphorus, total nitrogen, dissolved oxygen, and turbidity.

7. The method for evaluating the intensity of area source pollution based on multi-source data fusion according to claim 2, characterized in that, The specific steps for obtaining watershed land use types include: S11. Selecting satellite data sources and preprocessing the data, including radiometric calibration, atmospheric correction, geometric fine correction, image fusion, image mosaicking and cropping; S12. Extracting features, including spectral features, vegetation index, water body index, building index, impermeable surface index, texture features, spatial features, morphological features, time series features, topographic features, and contextual features; S13. Classifying land use types, dividing the target watershed into four categories—agricultural areas, forestry areas, industrial areas, and residential areas—based on a 3D-CNN model, and outputting the area and proportion of each type of land.

8. The method for evaluating the intensity of area source pollution based on multi-source data fusion according to claim 6, characterized in that, Step S2 includes the following steps: S21. Establish a unified spatiotemporal coordinate system for multi-source environmental data. The spatiotemporal coordinate system is established as follows: In the formula: —A set of spatiotemporal coordinate systems for the target watershed; — Column coordinates of the grid, ranging from 1 to ; — The row coordinates of the grid, ranging from 1 to ; —The total number of grid cells in the x-direction; —The total number of grid cells in the y direction; —Time variable; —Time window range, defined as ,in For the start time, The end time is specified; dimensionless processing is performed on the multi-source environmental data, and the formula for dimensionless processing of vegetation cover is as follows: In the formula: —Grid Normalized vegetation cover, with values ​​ranging from 0 to 1; —Grid Original vegetation cover; —Minimum vegetation coverage for the target watershed; —Maximum vegetation cover for the target watershed; Standardized surface roughness calculation formula is as follows: In the formula: —Grid Standardized surface roughness value; —Grid The original surface roughness value; —Minimum surface roughness of the watershed; —Maximum surface roughness of the watershed; the normalized soil permeability calculation formula is as follows: In the formula: —Grid Normalized soil permeability; —Grid Original soil permeability; —Minimum soil permeability in the watershed; —Maximum soil permeability in the watershed; the dimensionless formula for slope is as follows: In the formula: —Grid Standardized slope value; —Grid Original slope value, in units: ; —Maximum slope value, in units of: S22. Based on the land use types obtained in step S1, construct the pollution equivalent matrix C: In the formula: —The pollution equivalent coefficient for agricultural areas is 0.85; —The pollution equivalent coefficient for forestry areas is 0.15; —The pollution equivalent coefficient for the industrial zone is 0.6; —The pollution equivalent coefficient for residential areas is 0.7; the formula for calculating pollution source intensity is: In the formula: —Grid The intensity of the pollution source at time t; —Grid The area proportion of land of type k, where k=1 is agricultural area, k=2 is forestry area, k=3 is industrial area, and k=4 is residential area; —Pollution equivalent coefficient, k=1, 2, 3, 4; —Yearly accumulated days, valued as 1-365; —Seasonal offset; —Seasonal factor intensity coefficient, taken as 0.5; —Grid The intensity of rainfall; —The maximum rainfall intensity in the region; —Rainfall erosivity index, taken as 1.2; S23. Construction of environmental resistance model, including construction of resistance factor matrix and spatial principal component analysis, wherein the resistance factor matrix is: In the formula: —Grid The drag factor vector; the spatial principal component analysis includes constructing a local neighborhood data matrix, for the central grid point Take it Neighborhood: ;in, No more than , The smaller value in ; Calculate the local covariance matrix: In the formula: —Local covariance matrix; —Neighborhood mean vector ,in For the first Eigenvectors of neighborhood points; eigenvalue decomposition of the local covariance matrix to obtain eigenvalues ​​and eigenvectors; synthesis of drag coefficients: In the formula: —The weight of the j-th principal component, , Local covariance matrix The j-th eigenvalue is given, and the eigenvalues ​​are arranged in ascending order. Local covariance matrix The sum of all eigenvalues; —The eigenvector corresponding to the j-th eigenvalue of the local covariance matrix.

9. The method for evaluating the intensity of area source pollution based on multi-source data fusion according to claim 8, characterized in that, In step S3, the specific method for classifying the risk of non-point source pollution is as follows: when The risk level of non-point source pollution is low; when The risk level of non-point source pollution is medium; when The risk level of non-point source pollution is high.

Images