Spatial prediction of element content in black rock series high background soil and division of impact area

Abstract

The present invention belongs to the field of soil environment technology, and specifically relates to a method for spatial prediction of element content in black rock series high background soil and division of influence area. The method comprises the following steps: spatial prediction of element content is performed based on the digital terrain model of the study area, polygonal vector data of the black rock series and regional soil geochemical survey data of the study area; data in the study area is retained in the element spatial prediction result, and the background threshold is determined according to the probability distribution and the reference soil environment standard, and then the influence range of the black rock series is determined. The method of the present invention solves the problems of inaccurate spatial prediction of soil element content in black rock series high background areas and lack of basis for division of high background areas.

Description

Technical Field

[0001] The present invention relates to the field of soil environment technology, and in particular to a method for spatial prediction of element content in black rock series high-background soil and division of impact areas. Background Art

[0002] The description of the background technology in the present invention belongs to the related technology related to the present invention and is only used to illustrate and facilitate the understanding of the invention content of the present invention. It should not be understood that the applicant explicitly believes or presumes that the applicant believes that it is the prior art of the present invention on the filing date of the first application;

[0003] Soil space is continuous, but current survey or monitoring methods based on sampling tests can only obtain the soil element content corresponding to "discrete" "points". Existing technologies use spatial interpolation methods based on geographic information system technology to obtain content information of unknown points. Commonly used interpolation methods include: inverse distance interpolation, Kriging interpolation, radial basis function interpolation, etc. These methods are "broad-spectrum" methods, lack of specificity for the migration characteristics of elements released by weathering of black rock series, and their spatial prediction effects are extremely limited.

[0004] At present, there is a lack of a "theoretical" method for the spatial division of the black rock system high background area, or the black rock system influence range, and its boundary range is often determined subjectively by the black rock system area shown in the geological map;

[0005] The content of trace elements such as cadmium in the bedrock of black rock series is very high, and they are mobile during the weathering process, resulting in high background content of these trace elements in related soils. At present, there are two major problems in the risk control of soil heavy metal elements and the development of beneficial elements such as selenium in such high background areas: 1) The spatial variation of trace elements in soil is complex, which makes it difficult to obtain more objective and accurate spatial distribution information of elements through test data of limited samples even in the case of high sampling density; and existing interpolation methods, such as inverse distance interpolation and Kriging interpolation, lack targeted characterization of the spatial migration characteristics of trace elements in such high background areas, and their spatial prediction capabilities are limited; 2) Due to the spatial mobility of elements, the black rock series area given by the geological map does not overlap with the high background area of ​​soil elements. At present, there is a lack of effective spatial division method of high background areas, and there is no basis for determining the high background boundary; Summary of the invention

[0006] The purpose of the embodiment of the present invention is to provide a method for spatial prediction of element content in black rock series high background soil and division of impact areas, which solves the problem that there is a lack of basis for the division of high background areas in black rock series.

[0007] A method for spatial prediction of element content in high-background soil of black rock series and division of impact areas comprises the following steps: performing element spatial prediction based on a digital terrain model of a study area, polygonal vector data of black rock series and regional soil geochemical survey data of the study area; retaining data in the study area in the element spatial prediction results, determining a background threshold based on probability distribution and reference environmental soil standards, and further determining the impact range of the black rock series.

[0008] Furthermore, the black rock series polygonal vector data is extracted from a geological map or obtained through field survey.

[0009] Furthermore, the regional soil geochemical survey data of the study area were obtained from soil sampling tests and sorted into point vector data.

[0010] Furthermore, performing element space prediction specifically includes the following steps:

[0011] (1) The black rock distribution polygon vector file BlkShlVec is rasterized according to the DEM resolution and coordinate information. The grid value corresponding to the black rock is set to 1, and the others are set to 0 or background to form the raster data BlkShale;

[0012] (2) Based on the GIS hydrological analysis function, the DEM is used to fill the depressions and remove the sunken areas, thus forming the file S1DEMfill;

[0013] (3) Calculate the water flow direction of each grid in S1DEMfill. It is recommended to use the D8 algorithm and use the code to give the water flow direction of each grid value to form the file S2FlowDir;

[0014] (4) Taking BlkShale as the weight factor, use S2FlowDir to extract the water flow, that is, calculate the total number of corresponding water grids for each grid to form the grid data S3FlowAcc;

[0015] (5) S3FlowAcc is reclassified and binarized, with 0 set as the background and other values ​​greater than 0 set as 1, resulting in the formation of raster data S4BlkRgn;

[0016] (6) Taking S4BlkRgn as the "source" data and the original DEM to provide surface information, the "surface" path distance from each grid to the nearest "source" is calculated to form a grid file S5Path0; during the calculation, the input DEM elevation can be multiplied by a coefficient to highlight or suppress the influence of the terrain; when calculating the path, attention should be paid to the unity of the horizontal coordinate unit and the elevation unit of the input data, such as using meters. If the input data coordinates do not meet the above requirements, projection and coordinate conversion should be performed; then S5Path0 is subjected to a grid operation S5Path=S5Path0+n, where n>0, to obtain the grid file S5Path after removing the 0 value;

[0017] (7) Overlay the soil geochemical survey point data pnts with S5Path to obtain the path distance attribute Path;

[0018] (8) Based on the element content attribute data in pnts, the X element content and the Path attribute are logarithmically transformed to obtain logX and logPath; the regression relationship between logX and logPath is established: logX = a × logPath + b, and the regression coefficients a and b are calculated;

[0019] (9) Using the raster data S5Path to calculate the logarithm, the raster data S5logPath is obtained as the logPath variable in the above regression model. The above regression coefficients a and b are used to perform raster calculations, and the predicted results of the obtained logX values ​​form the raster data S6logpX;

[0020] (10) Superimpose the point data pnts and S6logpX, and obtain the grid data value as the attribute of pnts, the attribute name is S6logpX; subtract the logX of pnts from S6logpX to obtain the attribute reslogX, that is, reslogX = logX - S6logpX;

[0021] (11) Perform simple Kriging interpolation on the reslogX attribute in the pnts point data, and the interpolated raster data is recorded as S7reslogX;

[0022] (12) The raster data S6logpX and S7reslogX are operated: S8logX=S6logpX+S7reslogX to obtain the raster data S8logX, which is then subjected to an antilogarithmic transformation to obtain the raster data S8X; S8X is the spatial prediction result of the soil X content.

[0023] Furthermore, the background threshold is determined according to the probability distribution and the reference environmental soil standard, and then the influence range of the black rock series is determined, which specifically includes the following steps:

[0024] The grid data S8X was cropped to retain only the data in the study area; basic statistical analysis such as histogram was performed on the S8X grid values, and the background threshold was determined based on its probability distribution and in reference to relevant soil environmental quality standards;

[0025] Using a determined threshold, the area corresponding to the value greater than the threshold in S8X is converted into polygonal vector data, which is recorded as S9Rgn;

[0026] Select the polygon that intersects with BlkShlVecs in S9Rgn and save it as S10Bk as the influence range of the black rock system;

[0027] Furthermore, the X element is cadmium, selenium or other elements.

[0028] Furthermore, in the step of performing element space prediction, after step (6), S5Path or S5logPath can be used as an auxiliary variable to perform regression kriging analysis on the soil sampling test value X or logX.

[0029] The embodiments of the present invention have the following beneficial effects:

[0030] This application provides a targeted spatial prediction method based on the high background characteristics of black rock series, which significantly improves the accuracy of spatial prediction compared with conventional methods, and its interpolation results can more accurately characterize the spatial distribution state of elements; it also provides a prediction method for the influence range of black rock series, which solves the problem of lack of basis for the current division of high background areas of black rock series. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] Figure 1 It is a flow chart of the element space prediction (interpolation) technology in an embodiment of the present invention;

[0032] Figure 2 The black rock system influence area (high background area) is divided in the embodiment of the present invention;

[0033] Figure 3 The spatial distribution prediction of soil Cd given by the present invention in the embodiment of the present invention;

[0034] Figure 4 is the ordinary Kriging interpolation result in the embodiment of the present invention;

[0035] Figure 5 It is the scatter point relationship between the measured value (logCd) and the predicted value (S8logCd of the present invention, ordinary Kriging logCdordK) of soil cadmium in the embodiment of the present invention;

[0036] Figure 6 The range of the black rock series soil cadmium high background area divided by the method in the embodiment of the present invention;

[0037] Figure 7The BlkShale image of Example 1 in the embodiment of the present invention;

[0038] Figure 8 It is the log(S5Path) image of Example 1 in the embodiments of the present invention. DETAILED DESCRIPTION

[0039] The present application is further described below in conjunction with embodiments;

[0040] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, in the following description, different "one embodiment" or "embodiment" does not necessarily refer to the same embodiment; different embodiments may be replaced or combined, and for those skilled in the art, other implementation methods may be obtained based on these embodiments without creative work;

[0041] A method for spatial prediction of element content in high-background soil of black rock series and division of impact areas comprises the following steps: performing element spatial prediction based on a digital terrain model of a study area, polygonal vector data of black rock series and regional soil geochemical survey data of the study area; retaining data within the study area in the element spatial prediction results, determining a background threshold based on probability distribution and reference environmental soil standards, and further determining the impact range of the black rock series.

[0042] In some embodiments of the present invention, the black rock series polygonal vector data is extracted from a geological map or obtained through field survey.

[0043] In some embodiments of the present invention, regional soil geochemical survey data of the study area is obtained by soil sampling and testing and is organized into point vector data.

[0044] In some embodiments of the present invention, performing element space prediction specifically includes the following steps:

[0045] (1) The black rock distribution polygon vector file BlkShlVec is rasterized according to the DEM resolution and coordinate information. The grid value corresponding to the black rock is set to 1, and the others are set to 0 or background to form the raster data BlkShale;

[0046] (2) Based on the GIS hydrological analysis function, the DEM is used to fill the depressions and remove the sunken areas, thus forming the file S1DEMfill;

[0047] (3) Calculate the water flow direction of each grid in S1DEMfill. It is recommended to use the D8 algorithm and use the code to give the water flow direction of each grid value to form the file S2FlowDir;

[0048] (4) Taking BlkShale as the weight factor, use S2FlowDir to extract the water flow, that is, calculate the total number of corresponding water grids for each grid to form the grid data S3FlowAcc;

[0049] (5) S3FlowAcc is reclassified and binarized, with 0 set as the background and other values ​​greater than 0 set as 1, resulting in the formation of raster data S4BlkRgn;

[0050] (6) Taking S4BlkRgn as the "source" data and the original DEM to provide surface information, the "surface" path distance from each grid to the nearest "source" is calculated to form a grid file S5Path0; during the calculation, the input DEM elevation can be multiplied by a coefficient to highlight or suppress the influence of the terrain; when calculating the path, attention should be paid to the unity of the horizontal coordinate unit and the elevation unit of the input data, such as using meters. If the input data coordinates do not meet the above requirements, projection and coordinate conversion should be performed; then S5Path0 is subjected to a grid operation S5Path=S5Path0+n, where n>0, to obtain the grid file S5Path after removing the 0 value;

[0051] (7) Overlay the soil geochemical survey point data pnts with S5Path to obtain the path distance attribute Path;

[0052] (8) Based on the element content attribute data in pnts, the X element content and the Path attribute are logarithmically transformed to obtain logX and logPath; the regression relationship between logX and logPath is established: logX = a × logPath + b, and the regression coefficients a and b are calculated;

[0053] (9) Using the raster data S5Path to calculate the logarithm, the raster data S5logPath is obtained as the logPath variable in the above regression model. The above regression coefficients a and b are used to perform raster calculations, and the predicted results of the obtained logX values ​​form the raster data S6logpX;

[0054] (10) Superimpose the point data pnts and S6logpX, and obtain the grid data value as the attribute of pnts, the attribute name is S6logpX; subtract the logX of pnts from S6logpX to obtain the attribute reslogX, that is, reslogX = logX - S6logpX;

[0055] (11) Perform simple Kriging interpolation on the reslogX attribute in the pnts point data, and the interpolated raster data is recorded as S7reslogX;

[0056] (12) The raster data S6logpX and S7reslogX are operated: S8logX=S6logpX+S7reslogX to obtain the raster data S8logX, which is then subjected to an antilogarithmic transformation to obtain the raster data S8X; S8X is the spatial prediction result of the soil X content.

[0057] In some embodiments of the present invention, determining the background threshold according to the probability distribution and the reference environmental soil standard, and then determining the influence range of the black rock series specifically includes the following steps:

[0058] The grid data S8X was cropped to retain only the data in the study area; basic statistical analysis such as histogram was performed on the S8X grid values, and the background threshold was determined based on its probability distribution and in reference to relevant soil environmental quality standards;

[0059] Using a determined threshold, the area corresponding to the value greater than the threshold in S8X is converted into polygonal vector data, which is recorded as S9Rgn;

[0060] Select the polygon that intersects with BlkShlVecs in S9Rgn and save it as S10Bk as the influence range of the black rock system;

[0061] In some embodiments of the present invention, the X element is cadmium, selenium or other elements.

[0062] In some embodiments of the present invention, in the step of performing element space prediction, after step (6), S5Path or S5logPath can be used as an auxiliary variable to perform regression kriging analysis on the soil sampling test value X or logX.

[0063] The present invention utilizes terrain data (DEM) and black rock series coverage data (BlkShlVec or BlkShale) to generate black rock series impact weight spatial data S5Path or S5logPath.

[0064] The present invention uses S5Path or S5logPath as an auxiliary variable for spatial prediction or interpolation to improve interpolation accuracy.

[0065] Example 1

[0066] Input data: 1) Digital terrain model (DEM) of the study area, with a resolution of 8 meters, obtained through collection; 2) Black rock series polygonal vector data (BlkShlVec), extracted from the 1:50,000 geological map; 3) Regional soil geochemical survey data (pnts) of the study area, obtained by soil sampling and testing, originated from the 1:50,000 land quality geochemical evaluation project, and organized into point vector data based on GIS, where attributes include the content of elements such as Cd.

[0067] 1-element spatial prediction (interpolation), Figure 1 It is a flow chart of element space prediction (interpolation) technology;

[0068] 1.1 The black rock distribution polygon vector file BlkShlVec is rasterized according to the DEM resolution and coordinate information. The grid value corresponding to the black rock system is set to 1, and the others are set to 0 or background (NoData) to form the raster data BlkShale, which is displayed as follows: Figure 7 .

[0069] 1.2 Based on the GIS hydrological analysis function, the DEM is used to fill the depressions, remove the sunken areas, and generate the file S1DEMfill.

[0070] 1.3 Calculate the water flow direction of each grid in S1DEMfill. It is recommended to use the D8 algorithm and use the code to give the water flow direction of each grid value to form the file S2FlowDir.

[0071] 1.4 Take BlkShale as the weight factor and use S2FlowDir to extract the water flow, that is, calculate the total number of corresponding water grids for each grid to form the grid data S3FlowAcc.

[0072] 1.5 Reclassify S3FlowAcc and binarize it, set 0 as background (NoData) and other values ​​greater than 0 as 1, resulting in the formation of raster data S4BlkRgn.

[0073] 1.6 Take S4BlkRgn as the "source" data, use the original DEM to provide surface information, calculate the "surface" path distance from each raster to the nearest "source", and form a raster file S5Path0. This step can be completed by PathDistance in ArcGIS or other similar software modules. During the calculation, the input DEM elevation can be multiplied by a coefficient to highlight or suppress the influence of terrain. When calculating the path, it is necessary to pay attention to the unification of the horizontal coordinate unit and the elevation unit of the input data, such as using meters. If the input data coordinates do not meet the above requirements, projection and coordinate conversion are required. Then perform a raster operation S5Path=S5Path0+1 on S5Path0 to obtain the image of the raster file S5Path(log(S5Path)) after removing the 0 value. Figure 8 shown).

[0074] 1.7 Overlay the soil geochemical survey point data pnts with S5Path and obtain the path distance attribute Path from it.

[0075] 1.8 Based on the element content attribute data in pnts, such as Cd content, and Path attribute, logarithmic transformation is performed (it can be natural logarithm or logarithm with base 10, but it must be unified) to obtain logCd and logPath. Establish the regression relationship between logCd and logPath: logCd = a × logPath + b, and calculate the regression coefficients a and b.

[0076] 1.9 Use the raster data S5Path to calculate the logarithm to obtain the raster data S5logPath, which is used as the logPath variable in the above regression model. Use the above regression coefficients a and b to perform raster calculations, and the predicted results of the obtained logCd values ​​form the raster data S6logpCd.

[0077] 1.10 Overlay the point data pnts with S6logpCd, and obtain the grid data value as the attribute of pnts, with the attribute name S6logpCd; subtract the logCd of pnts from S6logpCd to obtain the attribute reslogCd, that is, reslogCd = logCd - S6logpCd.

[0078] 1.11 Perform simple kriging interpolation on the reslogCd attribute in the pnts point data, and the interpolated raster data is recorded as S7reslogCd.

[0079] 1.12 Calculate the raster data S6logpCd and S7reslogCd:

[0080] S8logCd = S6logpCd + S7reslogCd, and the grid data S8logCd is obtained. The antilogarithmic transformation is performed to obtain the grid data S8Cd. S8Cd is the spatial prediction result of soil Cd content ( Figure 3 ).

[0081] 2 Division of black rock system influence area (division of high background area), Figure 2 It is divided into the black rock system influence area (high background area);

[0082] 2.1 Cut the grid data S8Cd and keep only the data in the study area; perform basic statistical analysis such as histogram on the S8Cd grid values, and determine the background threshold based on its probability distribution and reference to relevant soil environmental quality standards, such as GB15618. For example, 0.3 mg / kg can generally be selected as the background threshold of Cd element.

[0083] 2.2 Using a determined threshold, the area in S8Cd that is larger than the threshold is converted into polygonal vector data, recorded as S9Rgn.

[0084] 2.3 Select the polygon that intersects with BlkShlVecs in S9Rgn and save it as S10Bk as the black rock system influence range or the black rock system high background area ( Figure 6 ).

[0085] Note:

[0086] 1) In the above, Cd represents the cadmium content; this method can also be used for other trace elements such as selenium, which can be replaced by the corresponding element symbols during operation.

[0087] 2) The above method flow can be implemented using GIS software such as ArcGIS with the assistance of statistical software such as EXCEL, SPSS or R language; it can also be implemented using a dedicated software module developed based on this patent.

[0088] Figure 3 This is the soil Cd content prediction result (S8Cd) in Example 1. It is used to illustrate that the soil Cd prediction result obtained by the corresponding method of the present invention can show more element spatial distribution details.

[0089] Figure 4 The spatial distribution of soil Cd predicted by the ordinary kriging method is not based on the technology described in the present invention. The results use the same soil geochemical survey point data pnts as in Example 1 as input data. Figure 3 , this figure does not provide more details on the spatial distribution of elements.

[0090] Figure 5 The left figure is a scatter plot of the measured soil Cd value (S8logCd) and the predicted value (logCd) obtained using the technology described in the present invention; the right figure is a scatter plot of the measured soil Cd value (logCdordK) and the predicted value (logCd) obtained using ordinary kriging. The scatter points in the former are more concentrated and closer to the 1:1 line, which shows that the technology described in the present invention can obtain more accurate prediction results than ordinary methods.

[0091] Figure 6 The black rock system distribution range and the black rock system influence range corresponding to Example 1. The present invention provides a method for delineating the black rock system influence range. This figure shows the corresponding results of this method in Example 1.

[0092] Table 1 Comparison of evaluation parameters obtained by regression analysis of measured values ​​and predicted values ​​of the two prediction methods

[0093]

[0094] Table 1 is a comparison of the evaluation parameters obtained by regression analysis of the measured values ​​and predicted values ​​of the two prediction methods. Among them, the F value and adjusted R 2The values ​​are greater than those obtained by ordinary kriging; at the same time, the AIC and BIC of the corresponding results of this patent are smaller than those obtained by ordinary kriging. The above indicates that the method provided by this patent has better prediction accuracy.

[0095] It should be noted that the above embodiments can be freely combined as needed; the above introduction is only the preferred embodiment of the present invention and is not intended to limit the present invention. For those skilled in the art, the present invention may have various changes and variations; any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for spatial prediction of element content in black rock series high background soil and division of impact areas, characterized in that: The method includes the following steps: performing element space prediction based on the digital terrain model of the study area, the polygonal vector data of the black rock series and the regional soil geochemical survey data of the study area; retaining the data in the study area in the element space prediction results, determining the background threshold based on the probability distribution and the reference environmental soil standard, and then determining the impact range of the black rock series; The black rock series polygonal vector data is extracted from a geological map or obtained through field survey; Regional soil geochemical survey data of the study area, obtained from soil sampling and testing, are organized into point vector data; The element space prediction specifically includes the following steps: (1) The black rock distribution polygon vector file BlkShlVec is rasterized according to the DEM resolution and coordinate information. The grid value corresponding to the black rock is set to 1, and the others are set to 0 or background to form the raster data BlkShale; (2) Based on the GIS hydrological analysis function, the DEM is used to fill the depressions and remove the sunken areas, thus forming the file S1DEMfill; (3) Calculate the water flow direction of each grid in S1DEMfill. It is recommended to use the D8 algorithm and use the code to give the water flow direction of each grid value to form the file S2FlowDir; (4) Taking BlkShale as the weight factor, use S2FlowDir to extract the water flow, that is, calculate the total number of corresponding water grids for each grid to form the grid data S3FlowAcc; (5) S3FlowAcc is reclassified and binarized, with 0 set as the background and other values ​​greater than 0 set as 1, resulting in the formation of raster data S4BlkRgn; (6) Using S4BlkRgn as the "source" data and the original DEM to provide surface information, the "surface" path distance from each grid to the nearest "source" is calculated to form a grid file S5Path0; during the calculation, the input DEM elevation can be multiplied by a coefficient to highlight or suppress the influence of the terrain; when calculating the path, attention should be paid to the unity of the horizontal coordinate unit and the elevation unit of the input data, such as using meters. If the input data coordinates do not meet the requirements, projection and coordinate conversion should be performed; Then, S5Path0 is subjected to a raster operation S5Path=S5Path0+n, where n>0, to obtain a raster file S5Path after removing the 0 value; (7) Overlay the soil geochemical survey point data pnts with S5Path to obtain the path distance attribute Path; (8) Based on the element content attribute data in pnts, the X element content and the Path attribute are logarithmically transformed to obtain logX and logPath; the regression relationship between logX and logPath is established: logX = a × logPath + b, and the regression coefficients a and b are calculated; (9) Using the raster data S5Path to calculate the logarithm, the raster data S5logPath is obtained as the logPath variable in the regression model. The above regression coefficients a and b are used to perform raster calculations, and the predicted results of the obtained logX values ​​form the raster data S6logpX; (10) Superimpose the point data pnts and S6logpX, and obtain the grid data value as the attribute of pnts, the attribute name is S6logpX; subtract the logX of pnts from S6logpX to obtain the attribute reslogX, that is, reslogX = logX - S6logpX; (11) Perform simple Kriging interpolation on the reslogX attribute in the pnts point data, and the interpolated raster data is recorded as S7reslogX; (12) The raster data S6logpX and S7reslogX are operated: S8logX=S6logpX+S7reslogX to obtain the raster data S8logX, which is then subjected to an antilogarithmic transformation to obtain the raster data S8X; S8X is the spatial prediction result of the soil X content.

2. The method for spatial prediction of element content in black rock series high background soil and division of influence areas according to claim 1 is characterized in that: Determining the background threshold based on the probability distribution and reference environmental soil standards, and then determining the impact range of the black rock series, specifically includes the following steps: The grid data S8X was cropped to retain only the data in the study area; basic statistical analysis such as histogram was performed on the S8X grid values, and the background threshold was determined based on its probability distribution and in reference to relevant soil environmental quality standards; Using a determined threshold, the area corresponding to the value greater than the threshold in S8X is converted into polygonal vector data, which is recorded as S9Rgn; Select the polygon that intersects with BlkShlVecs in S9Rgn and save it as S10Bk as the influence range of the black rock system.

3. The method for spatial prediction of element content in black rock series high background soil and division of influence areas according to claim 2 is characterized in that: Element X is cadmium, selenium or other elements whose content in black rock series is higher than Clarke value.

4. The method for spatial prediction of element content in black rock series high background soil and division of influence zones according to claim 3 is characterized in that: In the step of performing element space prediction, after step (6), S5Path or S5logPath can be used as an auxiliary variable to perform regression kriging analysis on the soil sampling test value X or logX.

Images