WO2023168781A1 - Soil cadmium risk prediction method based on spatial-temporal interaction relationship - Google Patents

Abstract

Disclosed in the present invention is a soil cadmium (Cd) risk prediction method based on a spatial-temporal interaction relationship, comprising: firstly, obtaining a data set related to soil Cd pollution in a target prediction area in a plurality of different historical periods, and extracting and cross-verifying data to obtain spatial point location distribution data of enterprises subjected to the soil Cd pollution in the target prediction area in different historical periods, as well as a soil Cd concentration and a soil pH at each sampling point in the target prediction area in different historical periods; then, performing bivariate local spatial autocorrelation analysis on a soil Cd risk level and a density value of polluted enterprises in different periods to obtain a spatial-temporal interaction relationship between soil Cd and the polluted enterprises; and finally, identifying a soil Cd risk in the target prediction area in different periods, and predicting a future soil Cd risk category by using a patch-generating land use simulation model. The present invention expands the existing method and thought for soil Cd partitioning, and has important theoretical and practical significance for further management and prevention of soil Cd pollution.

Description

A soil cadmium risk prediction method based on spatiotemporal interaction Technical field

The invention belongs to the field of soil cadmium (Cd) risk prediction, and specifically relates to a soil Cd risk prediction method based on spatiotemporal interaction.

Background technique

Due to the source-sink theoretical relationship between polluting enterprises and agricultural land soil pollution, more and more researchers have begun to study the relationship between soil heavy metals and polluting enterprises, thereby providing an evidence chain for traceability of soil heavy metal pollution. However, traditional source analysis methods, such as multivariate and geostatistical analysis, principal component analysis, positive definite matrix decomposition, etc., all ignore the spatiotemporal pollution characteristics of soil heavy metals to a certain extent. In some areas, opposite results may be obtained, which is difficult to be effective. and precise guidance on soil heavy metal management. In addition, soil heavy metal pollution is often closely related to human activities, especially polluting enterprises. The use of traditional source analysis methods does not provide a reliable spatiotemporal coupling relationship between sources and sinks, and cannot well solve the spatiotemporal variability of pollution. This makes At present, the prevention and control of heavy metal pollution in soil caused by changes in polluting enterprises has become quite difficult.

At present, some technical regulations or standards on the zoning management of heavy metal pollution in soil have been introduced. However, these zoning management methods only take into account the current level of heavy metal pollution in soil and ignore the future risks of heavy metal pollution in soil, such as soil pollution caused by polluting enterprises. risk. In any case, soil pollution management based on future risks can be used as an important reference for future soil management, and simulating the dynamic change characteristics of soil Cd future risk probability zones is of great significance for effective decision-making in current and future soil Cd pollution management. However, how to predict the future risk of soil heavy metal pollution is an urgent technical problem that needs to be solved.

Contents of the invention

The purpose of the present invention is to solve the problems existing in the prior art and provide a soil Cd risk prediction method based on spatiotemporal interaction.

The specific technical solutions adopted by the present invention are as follows:

A soil cadmium risk prediction method based on spatiotemporal interaction, the steps are as follows:

S1. Obtain multiple data sets related to soil Cd pollution in different historical periods in the target prediction area, and through data extraction and cross-validation, obtain the spatial point distribution data of soil Cd pollution enterprises in the target prediction area in different historical periods, and Soil Cd concentration and soil pH at each sampling point in the target prediction area in different historical periods;

S2. For each historical period, use the kernel density method to conduct kernel density analysis on the spatial point distribution data of soil Cd pollution enterprises in the target prediction area of the historical period, and obtain the density value of polluting enterprises. At the same time, based on the target prediction area of the historical period, For the soil Cd concentration and soil pH at each sampling point, the spatial distribution of soil Cd concentration and soil pH in the target prediction area was obtained through interpolation method. Combined with the soil Cd concentration risk screening values corresponding to different pH ranges, different differences in the target prediction area were determined. The soil Cd risk level of the location;

S3. For each historical period, after gridding the target prediction area, count the soil Cd risk level and polluting enterprise density value in each grid, and then double-check the soil Cd risk level and polluting enterprise density value for the historical period. Variable local spatial autocorrelation analysis is used to obtain the spatiotemporal interaction relationship between soil Cd and polluting enterprises;

S4. For each historical period, according to the corresponding soil Cd risk level, receptor vulnerability and the spatio-temporal interaction relationship, identify different levels of soil Cd risk areas in the target prediction area according to the preset risk area grading standards. ;

S5. Obtain the influencing factors of soil Cd risk evolution corresponding to each historical period in the target prediction area, combine the identification results of soil Cd risk areas in each historical period, and use the patch-generating Land Use Simulation Model (PLUS) ) to achieve prediction of future soil Cd risk areas within the target prediction area.

Preferably, in S1, the data set related to soil Cd pollution includes soil-related enterprise data, POI data, historical industrial enterprise data and multi-period high-resolution remote sensing images related to soil Cd pollution, as well as soil in different periods. Sampling point data; among them, soil-related enterprise data includes the name of the enterprise related to soil Cd pollution and the city where it belongs; POI data includes the enterprise name and its longitude and latitude information; historical industrial enterprise data includes the enterprise name and the year of production activities; soil sampling Point data include soil Cd concentration and soil pH data at different sampling points.

Further, in the S1, the method for obtaining the spatial point distribution data of the soil Cd pollution enterprises is:

S11. Convert the unstructured data in the soil-related enterprise data related to soil Cd pollution into structured data, thereby obtaining the names of soil-related enterprises related to soil Cd pollution; The company names in the industrial enterprise data and POI data are separately processed by word segmentation, entities at different levels in the enterprise names are extracted, and then similarity matching is performed. The soil-related enterprises that completely match the entities at all levels in the enterprise names are regarded as polluting enterprises, and the historical industries are The year of production activities of the polluting enterprise in the enterprise data and the geographical location information in the POI data are associated with the polluting enterprise, and the data of polluting enterprises including the geographical location and year of production activities are obtained.

S12. Extract high-resolution remote sensing images of each polluting enterprise and its surrounding sites in each production activity year from multi-period high-resolution remote sensing image information, and then use an image semantic segmentation model based on deep learning to map each polluting enterprise's corresponding Building features are extracted from high-resolution remote sensing image data to determine whether there are buildings or corporate factories in the image area in each production activity year, to achieve remote sensing verification of the spatial point distribution of polluting enterprises in different years, and to eliminate pollution that has not passed remote sensing verification. For enterprise data, the remaining polluting enterprise data are divided according to the year of production to obtain the spatial point distribution data of soil Cd pollution enterprises in the target prediction area in different historical periods.

Furthermore, in S11, the company name is segmented using the word segmentation engine jieba, and is divided into entities at four levels: administrative division, font size, industry, and organizational form. Then similarity matching is performed by calculating edit distance. If two companies If the four levels of entities in the names match exactly, the two company names are deemed to match.

Further, in S12, the image semantic segmentation model based on deep learning is the U-Net convolutional neural network model; when extracting building features from high-resolution remote sensing images, the longitude and latitude of polluting enterprises in different historical periods are used and the information on the corresponding remote sensing images to obtain the corresponding remote sensing images around the point; perform data preprocessing on the acquired remote sensing images, generate corresponding annotation data, and use the U-Net convolutional neural network model for training and modeling; and then use The trained U-Net convolutional neural network model performs image segmentation on high-resolution remote sensing images. Each pixel in the returned result map has two categories: building and non-building; the returned result map belonging to the same production year is compared with Compare the data of polluting enterprises. If there are buildings in the returned result map at the location of the polluting enterprise, it means that the polluting enterprise at that point really exists; if there are only non-buildings in the returned result map at the location of the polluting enterprise, it means that there are only non-buildings in the returned result map. The polluting enterprise at the point does not exist, and the polluting enterprise information is reviewed twice to determine whether there is a polluting enterprise at the point.

As a preferred method, the specific method of S3 is as follows: for each historical period, after dividing the soil Cd risk level, according to the boundary of the target prediction area, obtain the range enclosed by its minimum circumscribed rectangle; then use the minimum Starting from a certain vertex of the circumscribed rectangle, a grid and corresponding grid points are generated; the grid points are used to extract the nuclear density and soil Cd risk classification values of polluting enterprises in different historical periods, which respectively represent the aggregation degree and concentration of polluting enterprises in different historical periods in the boundary. Soil Cd risk level; perform bivariate local Moran analysis on the grid points with soil Cd risk level and polluting enterprise density value, and obtain the spatio-temporal interaction relationship between the two. The specific formula is as follows:

in

and

Variables a and b respectively represent soil Cd risk levels and polluting enterprise densities in different historical periods on grid i and grid j. Variables a and b respectively represent soil Cd risk levels and polluting enterprise densities in different historical periods in the grid; W _ij is the spatial weight matrix of grid i and grid j, which is obtained based on the Euclidean distance weight between grid i and grid j; I ^ab represents the local invariance between the a attribute at grid i and the b attribute at grid j. blue index; when I ^ab is significantly positive, then the soil Cd risk level at grid i has a significant local positive correlation with the density of polluting enterprises at grid j; when I ^ab is significantly negative, then the grid The soil Cd risk level at grid i has a significant local negative correlation with the density of polluting enterprises at grid j; when I ^ab is not significant, it is considered that the soil Cd risk level at grid i is related to the density of polluting enterprises at grid j Density has no obvious spatiotemporal interaction.

Preferably, in S4, in the preset risk area grading standards, the soil Cd risk area is divided into three risk control areas: high, medium and low, and risk uncertainty areas.

Preferably, the receptor vulnerability is population density.

Preferably, in S5, the patch generation land use simulation model extracts the changing characteristics of soil Cd risk areas at each level between historical periods, and uses a random forest algorithm to mine the influencing factors of each soil Cd risk evolution one by one. Obtain the change occurrence probability of soil Cd risk areas at each level and the contribution of each soil Cd risk evolution influencing factor to the change characteristics. Finally, under the constraints of the change occurrence probability, based on the latest soil Cd risk evolution influencing factors, achieve the target prediction area Prediction of future soil Cd risk areas in China.

As an option, the specific steps for the patch generation land use simulation model to predict future soil Cd risk areas in the target prediction area are as follows: First, perform overlay analysis on the soil Cd risk areas of each level in the two historical periods, and extract Change characteristics of each grade of soil Cd risk area; secondly, based on the random forest model, establish the occurrence probability of change in each grade of soil Cd risk area, and calculate the contribution rate of each soil Cd risk evolution influencing factor to each land use change, The occurrence probability is used to estimate the dynamic change trend of each level of soil Cd risk area in the grid unit; finally, the future level of soil Cd risk area is predicted by integrating the time series of soil Cd risk area data and the combined probability generated by the model, where The following formula estimates the combined probability that grid cell p will be occupied by a soil Cd risk area of level k in the future:

in

Indicates the combined probability that grid unit p converts from the soil Cd risk area of the original level to the soil Cd risk area of the target level k at the iteration time t; P _p,k indicates the occurrence of the soil Cd risk area of the target level k on the grid unit p The probability;

Represents the neighborhood effect of the soil Cd risk area of the target k level on the grid unit p corresponding to the iteration time t;

represents the inertia coefficient of the soil Cd risk area of target k level at iteration time t; sc _c→k represents the conversion cost from the original c level soil Cd risk area to the target k level soil Cd risk area; patch generation land use The simulation model uses a roulette selection mechanism to determine which land use type will occupy a grid cell.

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

This invention is based on the public information of soil-related enterprises and at the same time matches the data of pollution-related enterprises investigated and disclosed by the local government to obtain the data of soil-related enterprises related to pollution. There is no need to classify the enterprise data and determine whether the enterprise is a polluting enterprise in the future. Data cleaning and modeling are directly performed based on the data of polluting enterprises disclosed by the government where the target prediction area is located; in addition, the present invention obtains the density value based on the kernel density analysis of polluting enterprises in different periods, and conducts bi-variable local modeling with the soil Cd risk level in the corresponding period. Lan analysis analyzes the spatiotemporal interaction between the two, so that the discrete Cd pollution point data and polluting enterprise point data can more accurately explore the spatiotemporal interaction; finally, through the soil Cd risk zoning characteristics and the influencing factors of the risk zoning evolution in different periods , using a patch generation land use simulation model to model soil Cd risk zoning in different periods at present, and predict future soil Cd risk zoning characteristics, expanding the existing soil Cd zoning methods and ideas, and providing further guidance for the further management and prevention of soil Cd. Pollution has important theoretical and practical significance, and has the value of promotion and application.

Description of the drawings

Figure 1 is a partial display of the results in the embodiment, including (a) soil Cd risk level zoning in 2002, (b) population density zoning in 2002, (c) polluting enterprise density distribution in 2002, (d) soil Cd in 2002 Local Moran spatial interactions with polluting companies;

Figure 2 is a partial display of the results in the embodiment, including (a) soil Cd risk level zoning in 2012, (b) population density zoning in 2012, (c) polluting enterprise density distribution in 2012, (d) soil Cd in 2012 Local Moran spatial interactions with polluting companies;

Figure 3 is a partial display of the results in the embodiment, including (a) schematic diagram of soil Cd risk management zoning in 2002, (b) schematic diagram of soil Cd risk management zoning in 2012, (c) patch generation land use simulation model predicted in 2012 Schematic diagram of soil Cd risk management zoning, (d) Schematic diagram of soil Cd risk management zoning in 2022 predicted by the patch generation land use simulation model;

Figure 4 is the model accuracy of the patch generation land use simulation model in the embodiment to predict the risk zoning in 2012 based on the 2002 risk zoning.

Detailed ways

The present invention will be further elaborated and described below in conjunction with the accompanying drawings and specific embodiments. The technical features of various embodiments of the present invention can be combined accordingly as long as they do not conflict with each other.

This invention uses a similarity matching algorithm based on multi-source data to spatially identify polluting enterprises. It mainly establishes a similarity calculation method for data matching based on multi-source enterprise data. At the same time, it uses the U-Net model to establish remote sensing verification for identifying polluting enterprises. , and conduct bivariate local Moran index analysis on the obtained polluting enterprise data and soil Cd pollution data in different periods to construct the spatio-temporal interaction relationship between polluting enterprises and soil Cd risk classification, and use source-sink theory and spatio-temporal interaction relationship to construct Determine the characteristics of soil Cd risk zoning in different periods, and finally use the patch generation land use simulation model to predict future soil Cd risk zoning, thereby playing a directional role in the research on soil Cd pollution management zoning.

In a preferred embodiment of the present invention, a soil Cd risk prediction method based on spatiotemporal interaction is provided. The specific implementation steps are as follows:

S1. Obtain multiple data sets related to soil Cd pollution in different historical periods in the target prediction area, and through data extraction and cross-validation, obtain the spatial point distribution data of soil Cd pollution enterprises in the target prediction area in different historical periods, and Soil Cd concentration and soil pH at each sampling point in the target prediction area in different historical periods.

In this embodiment, the specific implementation of the above S1 step is as follows:

1) Data acquisition:

Obtain data sets related to soil Cd pollution, including: soil-related enterprise data, POI data, historical industrial enterprise data and multi-period high-resolution remote sensing images related to soil Cd pollution, as well as soil sampling point data in different periods; among them, involving Soil enterprise data includes the names of enterprises related to soil Cd pollution and the cities where they belong. AutoNavi POI data includes enterprise names and their longitude and latitude information; historical industrial enterprise data includes enterprise names and years of production activities; soil sampling point data includes different sampling Soil Cd concentration and soil pH data at points.

In this embodiment, pollution-related soil-related enterprise data, POI data, historical industrial enterprise data and multi-period high-resolution remote sensing images can be obtained using the http protocol and a public web page. The POI data uses Amap POI data; soil Sampling points are obtained from field sampling, which can be queried from soil census data or obtained through other means. Each sampling point should contain the soil Cd concentration and soil pH at the location of the sampling point. The reason for introducing soil pH is that the Cd risk in the soil is closely related to pH, so soil Cd concentration and soil pH need to be combined to comprehensively determine the soil Cd risk area. In addition, AutoNavi POI data and historical industrial enterprise data should try to cover enterprises within the target prediction area, so as to facilitate subsequent matching of more polluting enterprise data as much as possible.

2) Data preprocessing:

Convert the unstructured data in the soil-related enterprise data related to soil Cd pollution into structured data, thereby obtaining the names of soil-related enterprises related to soil Cd pollution; the obtained soil-related enterprise names, as well as historical industrial enterprise data and The company names in the POI data are segmented separately, entities at different levels in the company names are extracted, and then similarity matching is performed. The soil-related companies that completely match the entities at all levels in the company names are regarded as polluting companies, and the entities in the historical industrial company data are The year of production activities of the polluting enterprise and the geographical location information in the POI data are associated with the polluting enterprise, and the data of polluting enterprises including the geographical location and year of production activities are obtained.

In this embodiment, information in pollution-related enterprises exists in unstructured document data PDF, word or picture format, and optical character recognition (Optical Character Recognition, OCR) technology is used to identify it as based on the Office Open XML standard compressed file format, thereby converting it into structured data.

In this embodiment, the enterprise name can be segmented using the word segmentation engine jieba and divided into entities at four levels: administrative division, font size, industry, and organizational form. Then similarity matching is performed by calculating the edit distance. If two enterprise names If the entities at the four levels are completely matched, the two enterprise names are deemed to match. Jieba word segmentation is a dictionary-based word segmentation method. In order to obtain better segmentation results, it is necessary to use a custom dictionary to better identify unregistered words. By adding administrative division names and general name characteristic words to the custom dictionary , making the segmentation of names more accurate. In this embodiment, similarity matching calculates the edit distance of two characters to obtain the difference between strings. The smaller the difference, the higher the similarity. Edit distance refers to the minimum number of edit operations required between two strings to convert one into the other. If their edit distance is smaller, it means that their differences are smaller and their similarity is higher. Character operations allowed by edit distance include: a. Delete a character, b. Insert a character, c. Modify a character. The core idea is to perform minimal editing operations on a single character inside a string to obtain another string. The algorithm is mathematically defined as Formula 2.2.

where lev _a,b (i,j) represents the edit distance between the first i characters of string a and the first j characters of string b. When min(i,j)=0, the edit distance between a and b is max(i,j); when min(i,j)≠0, it is the minimum value of the following three situations: 1.lev _{a, b} (i-1,j)+1,2.lev _a,b (i,j-1)+1, and 3.

They respectively represent deletion of a _i , insertion of b _j and replacement of b _j ,

is an indicator function, which means that when a _i = b _j, it takes 0; when a _i ≠ b _j , its value is 1.

Based on edit distance, the similarity calculation formula between two strings is as follows:

In the formula, |a| and |b| are the string lengths of strings a and b, max(|a|,|b|) is the larger string length of strings a and b, lev _a,b ( i,j) is the minimum number of words of editing operations required to convert a to b.

3) Remote sensing verification:

Extract high-resolution remote sensing images of each polluting enterprise and its surrounding sites in each production activity year from multi-period high-resolution remote sensing image information, and then use an image semantic segmentation model based on deep learning to classify the high-resolution images corresponding to each polluting enterprise. High-resolution remote sensing image data is used to extract building features to determine whether there are buildings or corporate factories in the image area in each production activity year, to achieve remote sensing verification of the spatial point distribution of polluting enterprises in different years, and to eliminate data from polluting enterprises that have not passed remote sensing verification. , divide the remaining polluting enterprise data according to the year of production, and obtain the spatial point distribution data of soil Cd pollution enterprises in the target prediction area in different historical periods.

In this embodiment, the image semantic segmentation model based on deep learning can adopt the U-Net convolutional neural network model. When extracting building features from high-resolution remote sensing images, the corresponding remote sensing images around the point are obtained based on the longitude and latitude of the location of polluting enterprises in different historical periods and the information on the corresponding remote sensing images; then the obtained remote sensing images are Data preprocessing, generate corresponding annotation data and use U-Net convolutional neural network model for training modeling; then use the trained U-Net convolutional neural network model to perform image segmentation on high-resolution remote sensing images, and return the results Each pixel in the figure has two categories: building and non-building; finally, the returned result map belonging to the same production year is compared with the polluting enterprise data. If the location of the polluting enterprise has a building in the returned result map, it means that the polluting enterprise has a building in the returned result map. The polluting enterprise at the point does exist; if the location of the polluting enterprise only contains non-buildings in the returned result map, it means that the polluting enterprise at the point does not exist, and the polluting enterprise information will be reviewed twice to determine the point. Whether there are polluting enterprises. The secondary review here can be achieved through manual review or other forms.

The U-Net convolutional neural network model belongs to the existing technology. The model mainly extracts features from the image through the steps of convolution and pooling, and restores the feature map to the size of the original image through the step of deconvolution. Convolution and The specific calculation of pooling is as follows:

Among them, it is assumed that k is the convolution layer of the model,

Represents the j-th feature map generated in the k-th layer. Compare the feature map of x _i ^k-1 with the pixels in the corresponding convolution kernel

Convolution calculations are performed and summed, b represents the corresponding bias. F(X) is the activation function of this layer, and ReLU is selected as the activation function here.

Among them, down() is used as a downsampling function. Under the average pooling step, all pixels in a fixed-size pixel area are added, and the size of the final feature map becomes 1/n of the original. The specific structure and training process of the U-Net convolutional neural network model will not be described again.

S2. For each historical period, use the kernel density method to conduct kernel density analysis on the spatial point distribution data of soil Cd pollution enterprises in the target prediction area of the historical period, and obtain the density value of polluting enterprises. At the same time, based on the target prediction area of the historical period, For the soil Cd concentration and soil pH at each sampling point, the spatial distribution of soil Cd concentration and soil pH in the target prediction area was obtained through interpolation method. Combined with the soil Cd concentration risk screening values corresponding to different pH ranges, different differences in the target prediction area were determined. Soil Cd risk level of the location.

The soil Cd concentration risk screening values corresponding to different pH ranges can be determined based on relevant standards and specifications or expert experience. In this example, "Soil Environmental Quality - Agricultural Land Soil Pollution Risk Management and Control Standards" GB 15618-2018 is used to determine the soil Cd concentration risk screening values corresponding to different pH ranges, and then based on the soil Cd concentration risk of each point The screening value is used as a benchmark to determine whether each point exceeds the Cd risk screening value, so as to achieve the classification of soil Cd risk levels based on the judgment results.

S3. For each historical period, after gridding the target prediction area, count the soil Cd risk level and polluting enterprise density value in each grid, and then double-check the soil Cd risk level and polluting enterprise density value for the historical period. Through local spatial autocorrelation analysis of variables, the spatiotemporal interaction relationship between soil Cd and polluting enterprises was obtained.

In this embodiment, the specific implementation method of step S3 is as follows: for each historical period, after completing the division of soil Cd risk levels, according to the boundary of the target prediction area, obtain the range enclosed by its minimum circumscribed rectangle; then use this Starting from a certain vertex of the minimum circumscribed rectangle, a grid and corresponding grid points are generated; the grid points are used to extract the nuclear density of polluting enterprises and soil Cd risk classification values in different historical periods, which respectively represent the aggregation degree of polluting enterprises in different historical periods in the boundary. and soil Cd risk level; perform bivariate local Moran analysis on the grid points with soil Cd risk level and polluting enterprise density value, and obtain the spatio-temporal interaction relationship between the two. The specific formula is as follows:

in

and

Variables a and b respectively represent soil Cd risk levels and polluting enterprise densities in different historical periods on grid i and grid j. Variables a and b respectively represent soil Cd risk levels and polluting enterprise densities in different historical periods in the grid; W _ij is the spatial weight matrix of grid i and grid j, which is obtained based on the Euclidean distance weight between grid i and grid j; I ^ab represents the local invariance between the a attribute at grid i and the b attribute at grid j. blue index; when I ^ab is significantly positive, then the soil Cd risk level at grid i has a significant local positive correlation with the density of polluting enterprises at grid j; when I ^ab is significantly negative, then the grid The soil Cd risk level at grid i has a significant local negative correlation with the density of polluting enterprises at grid j; when I ^ab is not significant, it is considered that the soil Cd risk level at grid i is related to the density of polluting enterprises at grid j Density has no obvious spatiotemporal interaction.

S4. For each historical period, according to the corresponding soil Cd risk level, receptor vulnerability and the spatio-temporal interaction relationship, identify different levels of soil Cd risk areas in the target prediction area according to the preset risk area grading standards. .

In this embodiment, the population density in different historical periods in the target prediction area is used as the receptor vulnerability index. When the population density is 10 times higher than my country's average population density in that year, it is considered to be high vulnerability, otherwise it is considered low vulnerability. The specific levels of soil Cd risk areas can be set according to actual needs, and the standards for each classification can be adjusted based on relevant guidelines or specifications or expert experience. In this embodiment, in the preset risk area grading standards, the soil Cd risk area is divided into three risk control areas: high, medium and low, and risk uncertainty areas. In this embodiment, the source-sink theory and the spatio-temporal interaction determined in the previous step are used to determine the levels of soil Cd risk areas in different periods. The risk area grading standards adopted are as follows:

Table 1 Risk zone grading standards

S5. Obtain the influencing factors of soil Cd risk evolution corresponding to each historical period in the target prediction area, combine the identification results of soil Cd risk areas in each historical period, and use the patch generation land use simulation model (PLUS) to realize the future prediction of the target prediction area. Prediction of soil Cd risk areas.

In the present invention, the patch generation land use simulation model PLUS belongs to the existing model software. The download link is: https://github.com/HPSCIL/Patch-generating_Land_Use_Simulation_Model model software can run independently on Windows Vista/7/8/X64 Environment, PLUS model software can run independently in Windows Vista/7/8/X64-bit environment. The PLUS model extracts the change characteristics of soil Cd risk areas at each level between historical periods, and uses the random forest algorithm to mine the influencing factors of each soil Cd risk evolution one by one to obtain the change occurrence probability of soil Cd risk areas at each level and each soil The contribution of Cd risk evolution influencing factors to the change characteristics is ultimately constrained by the change occurrence probability and based on the latest soil Cd risk evolution influencing factors, the prediction of future soil Cd risk areas in the target prediction area is achieved. As a land use change simulation model, the PLUS model integrates land expansion analysis strategies and a cellular automaton model based on multi-type random patch seeds. It is an effective method for simulating dynamic changes in land use types. Therefore, based on the source-sink theoretical relationship, considering the soil Cd pollution level and receptor vulnerability level, the bivariate local Moran index was used to establish the spatiotemporal interaction between soil Cd and polluting enterprises, and at the same time, the patch generation land use simulation model was used to Predicting future soil Cd risks is of urgent practical significance for scientific and reasonable zoning management of soil Cd pollution risks, and for guiding policymakers and stakeholders to manage and control soil Cd pollution.

In this embodiment, the specific steps for the patch generation land use simulation model to predict future soil Cd risk areas in the target prediction area are as follows: First, the soil Cd risk areas of each level in the two historical periods are superimposed and analyzed. Extract the change characteristics of each level of soil Cd risk area; secondly, establish the occurrence probability of change in each level of soil Cd risk area based on the random forest model, and calculate the contribution rate of each soil Cd risk evolution influencing factor to each land use change , the occurrence probability is used to estimate the dynamic change trend of each level of soil Cd risk area in the grid unit; finally, the future level of soil Cd risk area is predicted by integrating the time series of soil Cd risk area data and the combined probability generated by the model, where The combined probability that grid unit p will be occupied by a k-level soil Cd risk zone in the future is estimated by the following formula:

in

Indicates the combined probability that grid unit p converts from the soil Cd risk area of the original level to the soil Cd risk area of the target level k at the iteration time t; P _p,k indicates the occurrence of the soil Cd risk area of the target level k on the grid unit p The probability;

Represents the neighborhood effect of the soil Cd risk area of the target k level on the grid unit p corresponding to the iteration time t;

represents the inertia coefficient of the soil Cd risk area of target k level at iteration time t; sc _c→k represents the conversion cost from the original c level soil Cd risk area to the target k level soil Cd risk area; patch generation land use The simulation model uses a roulette selection mechanism to determine which land use type will occupy a grid cell.

There is no doubt that the dominant soil Cd risk area level with the highest combination probability on the grid unit is prioritized for allocation conversion, but other soil Cd risk area levels with relatively low combination probability still have a chance to be assigned, even if the chance is small. To achieve this, the patch generation land use simulation model uses a roulette selection mechanism to determine which soil Cd risk zone class will occupy the grid cell. The roulette selection mechanism means that the areas occupied by different types of roulette wheels represent distribution probabilities. The more areas a soil Cd risk zone level occupies, the greater the probability of allocation, and it will not deprive other areas of the opportunity and possibility of being allocated.

Below, the soil cadmium risk prediction method based on spatiotemporal interaction in the above embodiment is applied to a specific example to further demonstrate the technical effects that can be achieved by the present invention, so that those skilled in the art can better understand the essence of the present invention. The specific implementation step framework in the following examples is as mentioned above S1 to S5, and will not be completely repeated. The specific implementation details and technical effects are mainly shown below.

Example

Select an area along the southeastern coast of China as the research area, and use the method shown in the aforementioned embodiments S1 to S5 of the present invention to predict soil cadmium risk. The specific steps are as follows:

First, based on the http protocol and the public web page, the pollution-related enterprise data in the study area, AutoNavi POI data, historical industrial enterprise data and multi-period high-resolution remote sensing images were obtained; among them, the soil-related enterprise data includes polluting enterprise information, The Amap POI data includes the latitude and longitude information of the GCJ-02 coordinate system, which is downloaded from the Amap Map API. The source of historical industrial enterprise data is the Chinese Industrial Enterprise Historical Database. The high-resolution remote sensing image data comes from the Resource and Environmental Science and Data of the Chinese Academy of Sciences. Center (https://www.resdc.cn/). Soil Cd concentration and pH data were obtained from field sampling points in the study area.

Information in pollution-related enterprises that exists in unstructured document data in PDF, word, or image formats will be identified using OCR technology as a compressed file format based on the Office Open XML standard and converted into structured data. Based on the polluting enterprise data, historical industrial enterprise data and Amap POI data, the enterprise names in the above data were segmented using the word segmentation engine jieba, and segmented into "city" + "enterprise name" + "industry name" + " "Suffix" four items; if the four items of the company names in the two data completely match after word segmentation processing, then the information will be returned to the polluting enterprise data, and the polluting enterprise data including the geographical location and the year of starting production activities will be obtained; at the same time, Convert the GCJ-02 coordinate system of POI data to the WGS-84 geographic coordinate system.

Secondly, based on the remote sensing image information of each polluting enterprise and its corresponding point, the corresponding remote sensing images around the point are extracted. Preprocess the data through the acquired remote sensing images, and generate corresponding annotation data to train and model the U-Net convolutional neural network model. details as follows:

The obtained remote sensing images and corresponding annotation maps of the corresponding points are cut into 256*256 images, and data enhancement operations such as flipping and cropping are performed to increase the diversity of the data. The training sample size is 10 as the model. The training input is, the initial learning rate is 5e-5, and the loss function uses binary cross entropy. After reaching 200K iterations, the model basically reaches the local optimum and converges. In the subsequent model verification process, the overall accuracy, recall rate, and F1-score of the model segmentation results were used as the accuracy evaluation of model performance.

TP represents the number of pixels that are classified correctly and are classified as buildings; TN represents the number of pixels that are classified correctly and are classified as background; FN represents the number of pixels that are classified as background and are classified as buildings. Number; FP represents the number of pixels whose category is background and is divided into building categories.

Then, the trained U-Net convolutional neural network model is tested and verified in the research area. The model performs image segmentation on the image at the corresponding point. Each pixel in the returned result map has only two categories: building and non-building. If there are buildings in the returned result map, it means that the enterprise at that point really exists; if there are only non-buildings in the returned result map, it means that the enterprise at that point does not exist; compare the results with the data of polluting enterprises. , if the two are the same, the valid polluting enterprise data will be automatically filtered out; if the two are different, the enterprise information will be reviewed twice to determine whether there is an enterprise at that point.

Then, use the kernel density method to perform spatial interpolation on the previously obtained point data of polluting enterprises in different periods, and set the output resolution to 100m to obtain the spatial distribution density map of polluting enterprises in different periods in the study area (Figure 1c and Figure 2c) , that is, kernel density analysis is completed. The inverse distance weighting method was used to predict the spatial distribution of soil Cd concentration in the study area in different periods, and the "Soil Environmental Quality - Agricultural Land Soil Pollution Risk Management and Control Standards" was used to classify soil Cd risk levels (Figure 1a and Figure 2a). A specific 100m resolution grid and corresponding grid points are generated based on the research boundary of the study area, and the soil Cd risk levels and polluting enterprise density values in the grid in different periods are calculated based on the grid points. Bivariate local spatial autocorrelation analysis was performed on soil Cd risk levels and polluting enterprise density values in different periods. The selected spatial neighbor relationship was K-Nearest neighbors, and the spatiotemporal interaction relationship between soil Cd and polluting enterprises in different periods was obtained (Figure 1d and Figure 1 2d).

Finally, as shown in Figure 2, based on the soil Cd risk level, receptor vulnerability (taking population density as an example (Figure 1b and Figure 2b)) and spatiotemporal interaction in the study area at different periods, refer to Table 1 for soil Cd risk levels in the study area at different periods. Cd for risk identification (Figure 3a and Figure 3b). Subsequently, by collecting the influencing factors of soil Cd risk evolution in the study area in different periods (including digital elevation model, density of polluting enterprises, PM2.5, average annual precipitation, average annual temperature, population density and nighttime lighting), the 2002 risk management zoning and impact Factors predicted the risk management partition in 2012 and built a model (Figure 3c). The overall accuracy, kappa coefficient, and producer accuracy and user accuracy of different levels of risk control areas were used as evaluation indicators to obtain the accuracy evaluation results of the model (Figure 4) ; Finally, the patch generation land use simulation model was used to predict future soil Cd risk categories in the study area (Figure 3d) for subsequent monitoring and management.

The above-described embodiment is only a preferred solution of the present invention, but it is not intended to limit the present invention. Those of ordinary skill in the relevant technical fields can also make various changes and modifications without departing from the spirit and scope of the present invention. Therefore, any technical solution obtained by adopting equivalent substitution or equivalent transformation shall fall within the protection scope of the present invention.

Claims (10)

1. A soil cadmium risk prediction method based on spatiotemporal interaction, characterized by the following steps:

   S1. Obtain multiple data sets related to soil Cd pollution in different historical periods in the target prediction area, and through data extraction and cross-validation, obtain the spatial point distribution data of soil Cd pollution enterprises in the target prediction area in different historical periods, and Soil Cd concentration and soil pH at each sampling point in the target prediction area in different historical periods;

   S2. For each historical period, use the kernel density method to conduct kernel density analysis on the spatial point distribution data of soil Cd pollution enterprises in the target prediction area of the historical period, and obtain the density value of polluting enterprises. At the same time, based on the target prediction area of the historical period, For the soil Cd concentration and soil pH at each sampling point, the spatial distribution of soil Cd concentration and soil pH in the target prediction area was obtained through interpolation method. Combined with the soil Cd concentration risk screening values corresponding to different pH ranges, different differences in the target prediction area were determined. The soil Cd risk level of the location;

   S3. For each historical period, after gridding the target prediction area, count the soil Cd risk level and polluting enterprise density value in each grid, and then double-check the soil Cd risk level and polluting enterprise density value for the historical period. Variable local spatial autocorrelation analysis is used to obtain the spatiotemporal interaction relationship between soil Cd and polluting enterprises;

   S4. For each historical period, according to the corresponding soil Cd risk level, receptor vulnerability and the spatio-temporal interaction relationship, identify different levels of soil Cd risk areas in the target prediction area according to the preset risk area grading standards. ;

   S5. Obtain the influencing factors of soil Cd risk evolution corresponding to each historical period in the target prediction area, combine the identification results of soil Cd risk areas in each historical period, and use the patch generation land use simulation model (PLUS) to realize the future prediction of the target prediction area. Prediction of soil Cd risk areas.
2. The soil cadmium risk prediction method based on spatiotemporal interaction as claimed in claim 1, wherein in S1, the data set related to soil Cd pollution includes soil-related enterprise data and POI related to soil Cd pollution. Data, historical industrial enterprise data and multi-period high-resolution remote sensing images, as well as soil sampling point data in different periods; among them, soil-related enterprise data includes the names of enterprises related to soil Cd pollution and the cities where they belong, and POI data includes the names of enterprises and its latitude and longitude information; historical industrial enterprise data includes enterprise names and years of production activities; soil sampling point data includes soil Cd concentration and soil pH data at different sampling points.
3. The soil cadmium risk prediction method based on spatiotemporal interaction according to claim 2, characterized in that, in the S1, the method for obtaining the spatial point distribution data of the soil Cd pollution enterprises is:

   S11. Convert the unstructured data in the soil-related enterprise data related to soil Cd pollution into structured data, thereby obtaining the names of soil-related enterprises related to soil Cd pollution; The company names in the industrial enterprise data and POI data are separately processed by word segmentation, entities at different levels in the enterprise names are extracted, and then similarity matching is performed. The soil-related enterprises that completely match the entities at all levels in the enterprise names are regarded as polluting enterprises, and the historical industries are The year of production activities of the polluting enterprise in the enterprise data and the geographical location information in the POI data are associated with the polluting enterprise, and the data of polluting enterprises including the geographical location and year of production activities are obtained.

   S12. Extract high-resolution remote sensing images of each polluting enterprise and its surrounding sites in each production activity year from multi-period high-resolution remote sensing image information, and then use an image semantic segmentation model based on deep learning to map each polluting enterprise's corresponding Building features are extracted from high-resolution remote sensing image data to determine whether there are buildings or corporate factories in the image area in each production activity year, to achieve remote sensing verification of the spatial point distribution of polluting enterprises in different years, and to eliminate pollution that has not passed remote sensing verification. For enterprise data, the remaining polluting enterprise data are divided according to the year of production to obtain the spatial point distribution data of soil Cd pollution enterprises in the target prediction area in different historical periods.
4. The soil cadmium risk prediction method based on spatio-temporal interaction as claimed in claim 3, characterized in that in S11, the company name is segmented using the word segmentation engine jieba, and is divided into four categories: administrative division, trade size, industry, and organizational form. The entities at each level are then matched by calculating the edit distance for similarity. If the entities at the four levels of the two company names completely match, the two company names are deemed to match.
5. The soil cadmium risk prediction method based on spatiotemporal interaction as claimed in claim 3, characterized in that, in the S12, the image semantic segmentation model based on deep learning is a U-Net convolutional neural network model; When extracting building features from high-resolution remote sensing images, the corresponding remote sensing images around the point are obtained based on the longitude and latitude of polluting enterprises in different historical periods and the information on the corresponding remote sensing images; data preprocessing is performed on the obtained remote sensing images to generate the corresponding Label the data and use the U-Net convolutional neural network model for training modeling; then use the trained U-Net convolutional neural network model to segment the high-resolution remote sensing images, and return each pixel in the result image There are two categories: building and non-building; compare the returned result map belonging to the same production year with the polluting enterprise data. If there is a building in the returned result map at the location of the polluting enterprise, it means that the polluting enterprise at that point really exists; If there are only non-buildings in the returned result map for the location of the polluting enterprise, it means that the polluting enterprise does not exist at that point, and the polluting enterprise information will be reviewed twice to determine whether there is a polluting enterprise at that point.
6. A soil Cd risk prediction method based on spatiotemporal interaction as claimed in claim 1, characterized in that the specific method of S3 is as follows: for each historical period, after dividing the soil Cd risk levels, according to The boundary of the target prediction area is obtained to obtain the range enclosed by its minimum circumscribed rectangle; then starting from a certain vertex of the minimum circumscribed rectangle, a grid and corresponding grid points are generated; the grid points are used to extract the polluting enterprise core values in different historical periods. Density and soil Cd risk classification values represent the aggregation degree of polluting enterprises and soil Cd risk level in different historical periods in the boundary respectively; bivariate local Moran analysis is performed on the grid points with soil Cd risk level and density value of polluting enterprises. The spatio-temporal interaction relationship between the two is obtained. The specific formula is as follows:

   

   in

   

   and

   

   Variables a and b respectively represent soil Cd risk levels and polluting enterprise densities in different historical periods on grid i and grid j. Variables a and b respectively represent soil Cd risk levels and polluting enterprise densities in different historical periods in the grid; W _ij is the spatial weight matrix of grid i and grid j, which is obtained based on the Euclidean distance weight between grid i and grid j; I ^ab represents the local invariance between the a attribute at grid i and the b attribute at grid j. blue index; when I ^ab is significantly positive, then the soil Cd risk level at grid i has a significant local positive correlation with the density of polluting enterprises at grid j; when I ^ab is significantly negative, then the grid The soil Cd risk level at grid i has a significant local negative correlation with the density of polluting enterprises at grid j; when I ^ab is not significant, it is considered that the soil Cd risk level at grid i is related to the density of polluting enterprises at grid j Density has no obvious spatiotemporal interaction.
7. A soil Cd risk prediction method based on spatio-temporal interaction as claimed in claim 1, characterized in that, in said S4, in the preset risk zone grading standards, the soil Cd risk zone is divided into high, medium, There are three low risk control areas and risk uncertainty areas.
8. A soil Cd risk prediction method based on spatiotemporal interaction according to claim 1, characterized in that the receptor vulnerability is population density.
9. A soil Cd risk prediction method based on spatiotemporal interaction as claimed in claim 1, characterized in that, in the S5, the patch generation land use simulation model extracts soil Cd risk areas of various levels between historical periods. change characteristics, and use the random forest algorithm to mine the influencing factors of soil Cd risk evolution one by one to obtain the change occurrence probability of soil Cd risk areas at each level and the contribution of each soil Cd risk evolution influencing factor to the change characteristics, and finally in the change Under the constraints of occurrence probability, based on the latest soil Cd risk evolution influencing factors, the prediction of future soil Cd risk areas in the target prediction area is realized.
10. A soil Cd risk prediction method based on spatiotemporal interaction as claimed in claim 1, characterized in that the specific steps for the patch generation land use simulation model to predict future soil Cd risk areas in the target prediction area are as follows : First, the soil Cd risk areas of each level in the two historical periods were overlaid and analyzed to extract the change characteristics of each level of soil Cd risk area; secondly, the occurrence of changes in each level of soil Cd risk area was established based on the random forest model. probability, and calculate the contribution rate of each soil Cd risk evolution influencing factor to each land use change. The occurrence probability is used to estimate the dynamic change trend of each grade of soil Cd risk area in the grid unit; finally, by integrating the time series of soil The combined probability generated by Cd risk area data and the model is used to predict future soil Cd risk areas of each level, where the combined probability that grid unit p will be occupied by k level soil Cd risk areas in the future is estimated by the following formula:

    

    in

    

    Indicates the combined probability that grid unit p converts from the soil Cd risk area of the original level to the soil Cd risk area of the target level k at the iteration time t; P _p,k indicates the occurrence of the soil Cd risk area of the target level k on the grid unit p The probability;

    

    Represents the neighborhood effect of the soil Cd risk area of the target k level on the grid unit p corresponding to the iteration time t;

    

    represents the inertia coefficient of the soil Cd risk area of target k level at iteration time t; sc _c→k represents the conversion cost from the original c level soil Cd risk area to the target k level soil Cd risk area; patch generation land use The simulation model uses a roulette selection mechanism to determine which land use type will occupy a grid cell.

Images