Hyperspectral-based soil heavy metal content inversion modeling method and system

Abstract

The present application relates to heavy metal content inversion technical field, specifically relates to a kind of soil heavy metal content inversion modeling method and system based on hyperspectral.The spectral reflectance and elevation information are extracted under the unified grid unit by fusing hyperspectral image and digital elevation model;The source emission intensity index is constructed based on the spectral feature determination vegetation cover coefficient, clay adsorption coefficient and mineral abundance coefficient, and combined with topographic height index;Further, the migration resistance coefficient is calculated using topographic slope, vegetation and clay parameters, and the source emission intensity index and migration resistance coefficient are superimposed and attenuated according to the hydrological flow direction along the topological sequence, to generate the deposition accumulation index representing the secondary deposition accumulation trend of heavy metals;Finally, the dynamic accumulation index and static surface medium parameters are input into the regression algorithm to construct a nonlinear inversion model.The method effectively couples the geological background and hydrological migration process, and improves the accuracy of heavy metal inversion in complex terrain area.

Description

Technical Field

[0001] This invention relates to the field of heavy metal content inversion technology, specifically to a method and system for modeling and inverting soil heavy metal content based on hyperspectral imaging. Background Technology

[0002] The spatial distribution of heavy metal content in soil is controlled by the interaction of multiple environmental factors. In existing hyperspectral remote sensing inversion studies, establishing a statistical regression model between surface spectral reflectance and heavy metal content is the mainstream method. This type of method is usually based on an implicit assumption that the surface heavy metal content mainly depends on the in-situ soil mineral composition and organic matter content, and that its spectral response is spatially homogeneous.

[0003] However, in hilly areas with significant topographic relief and mineral-grain composite zones, this assumption often deviates from reality. The distribution of heavy metals on the Earth's surface is influenced not only by the primary geological background (such as parent material) but also significantly reshaped by hydrodynamic migration and sedimentation processes driven by surface runoff. Under the influence of rainfall runoff, weathering products or pollutants located in high-altitude areas migrate downstream with water flow and tend to undergo secondary enrichment in low-lying areas with dense vegetation or high clay content. This process leads to a spatially non-stationary mapping relationship between "surface spectral characteristics" and "actual heavy metal content." For example, in erosion and depositional zones, even if similar spectral mineral characteristics are observed, the corresponding actual accumulation of heavy metals may differ significantly. This makes existing methods for inverting and modeling soil heavy metal content unreliable. Summary of the Invention

[0004] To address the issue of poor reliability in soil heavy metal content inversion modeling methods in related technologies, this invention provides a method and system for soil heavy metal content inversion modeling based on hyperspectral imaging. The specific technical solution adopted is as follows: This invention proposes a method for inverting and modeling soil heavy metal content based on hyperspectral imaging. The method includes: Within different raster cells of the region to be analyzed, acquire spatially resolution registered hyperspectral images and digital elevation models, and determine the elevation value and spectral reflectance map of each raster cell. Based on the spectral reflectance distribution of each grid cell in different bands of the spectral reflectance map, the vegetation cover coefficient, clay adsorption coefficient, and mineral abundance coefficient are determined; based on the numerical difference in elevation between each grid cell and its adjacent grid cells, the topographic height index of the grid cell is determined; and combined with the topographic height index and mineral abundance coefficient, the source emission intensity index is determined. Based on the topographic slope, vegetation cover coefficient, and clay adsorption coefficient of the grid cells obtained from the analysis of the digital elevation model, the migration resistance coefficient of the grid cells is determined; according to the hydrological flow direction analysis, the migration resistance coefficient and the source emission intensity index of the grid cells are superimposed to obtain the deposition accumulation index of each grid cell. By combining vegetation cover coefficient, clay adsorption coefficient, mineral abundance coefficient and sedimentation accumulation index into a preset regression algorithm for nonlinear mapping analysis, an inversion model is constructed to determine the heavy metal content in the soil.

[0005] Further, determining the vegetation cover coefficient, clay adsorption coefficient, and mineral abundance coefficient based on the spectral reflectance distribution of each grid cell in different bands of the spectral reflectance map includes: The vegetation cover coefficient is determined based on the difference in spectral reflectance of the preset vegetation-affected bands in the spectral reflectance diagram. The clay adsorption coefficient is determined based on the absorption depth of the pre-defined clay characteristic influence band in the spectral reflectance diagram. The mineral abundance coefficient is determined based on the integral area of ​​the absorption valleys of the pre-defined mineral characteristic influence bands in the spectral reflectance diagram.

[0006] Furthermore, the preset vegetation influence bands include the near-infrared band and the red band, with a center wavelength of 670 nm for the red band and 860 nm for the near-infrared band; the preset clay characteristic influence bands include the Al-OH characteristic band with a wavelength of 2200 nm; the preset mineral characteristic influence bands include the mineral absorption band in the 800 nm-1000 nm band; the methods for obtaining the vegetation cover coefficient, clay adsorption coefficient, and mineral abundance coefficient include: Calculate the reflectance difference between the center wavelength of the near-infrared band and the center wavelength of the red band; The reflectance difference was normalized and used as the vegetation cover coefficient. The absorption depth of the trough in the preset clay characteristic influence band is obtained and normalized as the clay adsorption coefficient. In a spectral reflectance graph with wavelength as the abscissa and reflectance as the ordinate, the spectral curves of the mineral absorption bands and the area formed by the spectral curves and the continuum are obtained based on definite integrals, and the areas are normalized to obtain the mineral abundance coefficients.

[0007] Furthermore, based on the numerical difference in elevation values ​​between each raster cell and its adjacent raster cells, the terrain height index of the raster cell is determined, including: In the digital elevation model, the mean of all elevation values ​​of each grid cell within a preset neighborhood window is calculated and used as the neighborhood mean. The difference between the elevation value of a grid cell and the mean of its neighborhood is used as the terrain height index.

[0008] Furthermore, by combining the topographic height index and mineral abundance coefficient, the source emission intensity index is determined, including: Anomaly detection was performed on the mineral abundance coefficients of all raster cells using the three-standard-deviation principle, with the sum of the mean and three standard deviations used as the anomaly detection threshold. If the mineral abundance coefficient of any grid cell is greater than the anomaly judgment threshold and the terrain height index is greater than 0, then the source emission intensity index of the grid cell is set to the mineral abundance coefficient. Otherwise, the source emission intensity index of the grid cell is set to a preset background constant.

[0009] Furthermore, based on the terrain slope, vegetation cover coefficient, and clay adsorption coefficient of the raster cells obtained from the digital elevation model analysis, the migration resistance coefficient of the raster cells is determined, including: The environmental inhibition index is obtained by weighted summation of the vegetation cover coefficient and the clay adsorption coefficient. Based on the Horn algorithm, slope analysis is performed on the elevation values ​​of grid cells within a preset neighborhood window in the digital elevation model to determine the terrain slope; after converting the terrain slope to radians, the sum of the sine value and the preset safety constant term is calculated as the slope influence index. The ratio of the environmental inhibition index to the slope influence index is calculated and used as the migration resistance coefficient.

[0010] Furthermore, based on hydrological flow direction analysis, the migration resistance coefficient and source emission intensity index of the grid cells are superimposed to obtain the deposition accumulation index for each grid cell, including: Based on the hydrological flow direction analysis, the grid cells are sorted according to the flow direction to obtain the hydrological topology sequence along the hydrological flow direction; The migration resistance coefficient and source emission intensity index of each grid cell are superimposed according to the hydrological topological sequence to obtain the deposition accumulation index of each grid cell.

[0011] Furthermore, the migration resistance coefficient and source emission intensity index of the grid cells are superimposed according to the hydrological topological sequence to obtain the deposition accumulation index of each grid cell, including: The source emission intensity is used as the initial excitation index for the grid cell; Following the order of the hydrological topology sequence, the migration resistance coefficient of the first grid cell in the sequence is used as the exponential term, and the initial excitation index is used as the analytical excitation index of the first grid cell to calculate the exponential saturation decay, thus obtaining the deposition accumulation index of the first grid cell. The difference between the initial excitation index and the depositional accumulation index is calculated as the downstream transport index. The sum of the initial excitation index of the next grid cell in the sequence and the downstream transport index of all previous grid cells flowing to the next grid cell is used as the analytical excitation index of the next grid cell. This process is continued downstream to calculate the depositional accumulation index of each grid cell.

[0012] Furthermore, the preset regression algorithm is the random forest regression algorithm.

[0013] On the other hand, a soil heavy metal content inversion modeling system based on hyperspectral imaging is also provided. The system includes a memory, a processor, and a computer program stored in the memory and running on the processor. When the processor executes the computer program, it implements the steps of the method described in any of the foregoing descriptions.

[0014] The present invention has the following beneficial effects: This invention effectively overcomes the non-stationarity of the "spectrum-content" mapping caused by neglecting the hydraulic migration process in existing inversion methods by introducing the sedimentary accumulation index as a key process constraint feature. The sedimentary accumulation index constructs hydrological flow direction based on a digital elevation model and integrates surface migration resistance characterized by vegetation cover coefficient and clay adsorption coefficient, as well as source emission intensity determined by topographic height index and mineral abundance coefficient. Thus, it quantifies the secondary enrichment trend of heavy metals driven by runoff without the need for complex hydrological simulation. The dynamic accumulation index and static surface medium parameters are jointly input into the regression algorithm, enabling the model to distinguish areas with similar spectral characteristics but different hydrological locations (such as erosion areas and sedimentary areas), significantly improving the reliability and spatial accuracy of soil heavy metal content inversion in areas with large topographic relief. Attached Figure Description

[0015] To more clearly illustrate the technical solutions and advantages in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0016] Figure 1 The flowchart illustrates a method for inverting and modeling soil heavy metal content based on hyperspectral imaging, as provided in one embodiment of the present invention. Detailed Implementation

[0017] To further illustrate the technical means and effects adopted by the present invention to achieve its intended purpose, the following, in conjunction with the accompanying drawings and preferred embodiments, details the specific implementation, structure, features, and effects of a hyperspectral-based soil heavy metal content inversion modeling method and system proposed according to the present invention. In the following description, different "one embodiment" or "another embodiment" do not necessarily refer to the same embodiment. Furthermore, specific features, structures, or characteristics in one or more embodiments can be combined in any suitable form.

[0018] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0019] The following description, in conjunction with the accompanying drawings, details a specific scheme for a soil heavy metal content inversion modeling method based on hyperspectral imaging provided by this invention.

[0020] Please see Figure 1 The diagram illustrates a flowchart of a method for inverting and modeling soil heavy metal content based on hyperspectral imaging, according to an embodiment of the present invention. The method includes: S101: Within different raster cells of the area to be analyzed, acquire the spatially resolution-registered hyperspectral image and digital elevation model, and determine the elevation value and spectral reflectance map of each raster cell.

[0021] In this embodiment of the invention, a hyperspectral sensor can be used to acquire hyperspectral images and perform topographic mapping to obtain a digital elevation model.

[0022] Because hyperspectral sensors and topographic mapping sensors differ in imaging mechanisms and spatial resolution, this embodiment of the invention acquires a radiometrically calibrated and atmospherically corrected hyperspectral image of the area to be analyzed, defines its raster coordinate system as the global computational space, and obtains a digital elevation model (DEM) for the same area. Then, bicubic interpolation is used to resample the elevation values ​​in the DEM to match the spectral reflectance. Figure 1 To achieve the same spatial resolution (e.g., 30 meters), subpixel-level geometric registration is performed based on ground control points to ensure precise spatial correspondence between the two. The acquisition of hyperspectral imagery and digital elevation models, as well as their registration at the same spatial resolution, are well-known facts and will not be elaborated upon further.

[0023] In applications of inverting soil heavy metal content based on hyperspectral remote sensing, if the hyperspectral image and topographic data are not strictly aligned in space, it will lead to a misalignment between the surface material information and the analysis unit of the topographic-driven process, resulting in systematic errors in subsequent coupled modeling due to the "spectral-topographic" mismatch.

[0024] Therefore, this embodiment of the invention focuses on constructing spatial consistency of multi-source remote sensing data. Within the area to be analyzed, hyperspectral images and digital elevation models that have completed spatial resolution registration are acquired, and the elevation value and corresponding spectral reflectance map of each raster cell are simultaneously determined accordingly. This achieves precise alignment of hyperspectral land cover information and topographic morphology information at the pixel level, providing a unified, reliable, and geographically consistent basic data framework for the entire inversion model.

[0025] S102: Based on the spectral reflectance distribution of each grid cell in different bands of the spectral reflectance map, determine the vegetation cover coefficient, clay adsorption coefficient, and mineral abundance coefficient; based on the numerical difference in elevation between each grid cell and its adjacent grid cells, determine the topographic height index of the grid cell; and combine the topographic height index and the mineral abundance coefficient to determine the source emission intensity index.

[0026] In regions with complex topography, the primary sources and secondary migration processes of heavy metals in soil are intertwined. Relying solely on a single spectral or topographical indicator makes it difficult to effectively identify potential sources of material release, leading to misjudgments of heavy metal enrichment mechanisms. The surface medium plays a dual regulatory role in heavy metal migration: vegetation increases surface roughness, thereby hindering runoff; clay minerals in the soil provide chemical sites for adsorbing heavy metal ions; and the abundance of specific minerals indicates the primary geological background.

[0027] In this embodiment of the invention, the source emission potential of each grid cell is quantitatively characterized by analyzing the response characteristics of different bands in the hyperspectral reflectance map to extract the vegetation cover coefficient, clay adsorption coefficient, and mineral abundance coefficient, respectively, so as to reflect the composition and adsorption capacity of the surface medium. At the same time, the topographic height index is calculated based on the difference between the grid cell and its neighboring elevation values, and it is fused with the mineral abundance coefficient to construct the source emission intensity index.

[0028] Furthermore, in some embodiments of the present invention, the vegetation cover coefficient, clay adsorption coefficient, and mineral abundance coefficient are determined based on the spectral reflectance distribution of each grid cell in different bands of the spectral reflectance map. This includes: determining the vegetation cover coefficient based on the spectral reflectance differences of preset vegetation-affected bands in the spectral reflectance map; determining the clay adsorption coefficient based on the absorption depth of preset clay-affected bands in the spectral reflectance map; and determining the mineral abundance coefficient based on the integral area of ​​the absorption valleys of preset mineral-affected bands in the spectral reflectance map.

[0029] The first step is to calculate the vegetation cover coefficient by calculating the reflectance difference between the center wavelength of the near-infrared band and the center wavelength of the red band; then, the reflectance difference is normalized and used as the vegetation cover coefficient.

[0030] The red light band corresponds to the strong absorption peak of chlorophyll, while the near-infrared band corresponds to the strong reflection of mesophyll cell structure. Therefore, the normalized difference calculation between the red and near-infrared bands is based on the unique spectral response characteristics of green vegetation. In the red light band (approximately 670 nm), chlorophyll has a strong absorption effect on incident light, resulting in a significant reduction in reflectivity; while in the near-infrared band (approximately 860 nm), the spongy tissue structure of mesophyll cells causes strong multiple scattering and reflection of incident light, maintaining a high level of reflectivity.

[0031] In contrast, the reflectance difference between bare soil, rocks, and other non-vegetated backgrounds is relatively small in the two bands mentioned above. By calculating the difference between the two and normalizing it, the contrast between vegetation signals and background noise can be effectively enhanced. On the other hand, the normalization process can effectively offset the shadow effect caused by topographic relief and the difference in light intensity caused by changes in solar altitude angle, thereby extracting a purer spectral feature that reflects the surface vegetation coverage and growth status.

[0032] Specifically, in some embodiments of the present invention, the difference between the reflectance at wavelength 860nm and the reflectance at wavelength 670nm is directly calculated using the red light band (center wavelength approximately 670nm) and the near-infrared band (center wavelength approximately 860nm) in the spectral reflectance diagram as the reflectance difference. Then, the vegetation cover coefficient is obtained by performing maximum and minimum value normalization processing.

[0033] Of course, in other embodiments of the present invention, the red light band can also be determined based on the center wavelength. For example, the red light band at 670nm can specifically be 620nm-720nm. The average reflectance value of the corresponding band is taken as the red light reflectance, and the average reflectance value of the corresponding near-infrared band (810nm-910nm) is taken as the near-infrared light reflectance. The difference between the near-infrared light reflectance and the red light reflectance is normalized to the maximum and minimum values ​​to obtain the vegetation cover coefficient.

[0034] It should be noted that the normalization process in the embodiments of the present invention can specifically be, for example, maximum and minimum value normalization. Furthermore, the normalization in subsequent steps can all employ maximum and minimum value normalization. In other embodiments of the present invention, other normalization methods can be selected based on the specific range of values, which will not be elaborated further. The maximum and minimum value normalization in the embodiments of the present invention aims to perform standardization processing. The maximum and minimum values ​​of the normalization process can be set according to the specific scenario and data numerical characteristics. Adjusting, calibrating, or optimizing the maximum and minimum values ​​does not constitute a limitation of the present invention.

[0035] In this embodiment of the invention, when calculating the vegetation cover coefficient, the maximum and minimum values ​​can be selected based on the overall characteristics of all grid cells and normalized. The vegetation cover coefficient is used to characterize the physical blocking ability of surface vegetation on runoff; a larger value indicates a stronger blocking effect. During the model application phase, the maximum and minimum values ​​used for normalization should be fixed as the global extreme values ​​obtained from the training dataset.

[0036] In this embodiment of the invention, the preset clay characteristic influence band includes the Al-OH characteristic band with a wavelength of 2200nm. The absorption depth of the valley in the preset clay characteristic influence band is obtained and normalized as the clay adsorption coefficient.

[0037] The absorption depth at the Al-OH characteristic absorption valley near 2200 nm (the reflectance valley value at the wavelength closest to 2200 nm, commonly found in clay minerals such as kaolinite and montmorillonite) is extracted. Specifically, the wavelength positions of the left and right shoulders of the absorption valley are determined, a local continuum is constructed, and the ratio of the reflectance at the bottom wavelength of the absorption valley to the corresponding value on the continuum is calculated. A larger absorption depth indicates a more significant absorption characteristic (i.e., higher clay content). Obtaining the absorption depth is a well-known technique in this field and will not be elaborated upon further.

[0038] Among them, the wavelength positions of the left and right shoulders of the absorption valley are used to construct a local continuum, which represents the two local maxima with the highest reflectivity within the preset clay characteristic wavelength range (such as 2100nm-2300nm), and are respectively used as the wavelength positions of the left and right shoulders of the absorption valley.

[0039] In this embodiment of the invention, the crystal field transition absorption characteristics of the localized spectrum near the 900nm band (common in iron oxides such as hematite and goethite, often as associated minerals of heavy metals) are used as the preset mineral characteristic influence band. This preset mineral characteristic influence band includes the mineral absorption band in the 800nm-1000nm band. In a spectral reflectance diagram with wavelength as the abscissa and reflectance as the ordinate, the area formed by the spectral curve of the mineral absorption band and the continuum is obtained based on definite integrals. This area is then normalized to obtain the mineral abundance coefficient. The mineral abundance coefficient characterizes the relative abundance of primary surface minerals and is used to assist in identifying potential sources of heavy metals (high background areas).

[0040] In this context, the continuum in the mineral abundance coefficient calculation process is defined as the line connecting the local maxima of spectral reflectance on both sides of the 800nm-1000nm band.

[0041] It should be noted that the three parameters, vegetation cover coefficient, clay adsorption coefficient, and mineral abundance coefficient, are all normalized values ​​used to achieve data standardization. The normalization method is specifically maximum and minimum value normalization, and the minimum and maximum values ​​of the normalization are obtained from the values ​​corresponding to all raster cells.

[0042] Through the above steps, each grid cell in the entire domain is assigned three dimensionless static medium properties, namely vegetation cover coefficient, clay adsorption coefficient, and mineral abundance coefficient, providing standardized input parameters for subsequent spatial process modeling.

[0043] Heavy metal pollution is often accompanied by the enrichment of specific minerals, and is more prone to weathering and erosion in higher-altitude areas, spreading downstream with runoff. Therefore, in this embodiment of the invention, the topographic height index of a grid cell is determined by the numerical distribution of the elevation values ​​of each grid cell and its adjacent grid cells.

[0044] Furthermore, in some embodiments of the present invention, determining the terrain height index of a grid cell based on the numerical difference in elevation values ​​between each grid cell and its neighboring grid cells includes: in the digital elevation model, calculating the average of all elevation values ​​of each grid cell within a preset neighborhood window as the neighborhood average; and using the difference between the elevation value of the grid cell and the neighborhood average as the terrain height index.

[0045] The preset neighborhood window can be, for example, a 3×3 neighborhood window. The average elevation value of all grid cells within the preset neighborhood window is calculated as the neighborhood mean, which represents the overall neighborhood elevation level.

[0046] Then, the difference between the elevation value of the raster cell and the mean of its neighborhood is calculated as the terrain height index. A terrain height index greater than 0 indicates that the raster cell is higher than the overall level of its neighborhood, belonging to local highlands (such as ridges or slope tops), and possesses the topographical conditions to become an erosion source. Subsequently, combined with the mineral abundance coefficient, a statistical threshold method is used to identify potential high-potential material sources.

[0047] Furthermore, in some embodiments of the present invention, the source emission intensity index is determined by combining the terrain height index and the mineral abundance coefficient, including: performing anomaly judgment on the mineral abundance coefficient of all grid cells according to the three-standard-deviation principle, and using the sum of the mean and the three-standard-deviation as the anomaly judgment threshold; if the mineral abundance coefficient of any grid cell is greater than the anomaly judgment threshold and the terrain height index is greater than 0, then the source emission intensity index of the grid cell is set to the mineral abundance coefficient; otherwise, the source emission intensity index of the grid cell is set to a preset background constant.

[0048] In this embodiment of the invention, the mean value of the mineral abundance coefficients of all raster cells is calculated. with standard deviation Then, based on the three-standard-deviation principle, the default threshold for anomaly detection was set. For any grid cell, if both of the following conditions are met simultaneously: the terrain height index is greater than 0 and the mineral abundance coefficient is greater than the anomaly judgment threshold, then the grid cell is determined to be a potential source of high-potential energy matter, and the mineral abundance coefficient is directly used as the source emission intensity index; otherwise, the source emission intensity index is used as a preset background constant.

[0049] The preset background constant is a constant value under the background characteristics. Optionally, it can be less than the source emission intensity index corresponding to the high potential energy material source. Specifically, in this embodiment of the invention, the minimum value can be selected from the source emission intensity indices of all high potential energy material sources, and half of the minimum value can be used as the preset background constant.

[0050] Through the embodiments of the present invention, a comprehensive identification of potential heavy metal release source areas is achieved, organically combining geological background (minerals) and topographic location (relative elevation), providing a physically meaningful initial driving force input for subsequent migration process modeling.

[0051] S103: Based on the topographic slope, vegetation cover coefficient and clay adsorption coefficient of the grid cells obtained by the digital elevation model, determine the migration resistance coefficient of the grid cells; according to the hydrological flow direction analysis, superimpose the migration resistance coefficient and source emission intensity index of the grid cells to obtain the deposition accumulation index of each grid cell.

[0052] In areas with undulating terrain, heavy metal migration is influenced by a combination of physical resistance and chemical adsorption by the surface medium. If material transport is simulated solely based on topographic slope or flow direction, while ignoring the inhibitory effects of vegetation and clay on the migration process, the estimation of depositional accumulation trends will be distorted. To address this issue, this invention constructs a resistance and accumulation model that conforms to the actual physical mechanisms of migration. Based on the topographic slope extracted from the digital elevation model, and integrating the vegetation cover coefficient and clay adsorption coefficient obtained from hyperspectral inversion, a migration resistance coefficient characterizing the degree of migration obstruction is determined.

[0053] Furthermore, by combining hydrological flow direction analysis, the resistance coefficient and the source emission intensity index are superimposed and accumulated along the flow path to generate a sedimentary accumulation index reflecting the secondary enrichment potential of heavy metals. Through this step, a quantitative characterization of the migration-deposition process of heavy metals on the surface is achieved, enabling the calculation of sedimentary accumulation trends to incorporate both topographic driving logic and surface medium regulation effects, providing key process constraint characteristics for subsequent high-precision inversion.

[0054] First, based on the terrain slope, vegetation cover coefficient, and clay adsorption coefficient of the raster cells obtained from the analysis of the digital elevation model, the migration resistance coefficient of the raster cells is determined. This includes: weighting and summing the vegetation cover coefficient and clay adsorption coefficient to obtain an environmental inhibition index; performing slope analysis on the elevation values ​​of the raster cells within a preset neighborhood window in the digital elevation model based on the Horn algorithm to determine the terrain slope; converting the terrain slope to radians and calculating the sum of the sine value and the preset safety constant term as the slope influence index; and calculating the ratio of the environmental inhibition index to the slope influence index as the migration resistance coefficient.

[0055] It should be noted that the weight of the vegetation cover coefficient... Weight of clay adsorption coefficient You can directly use default values, such as the weight values ​​for vegetation cover coefficient and clay adsorption coefficient. and They are all the same and all equal to 1, or they can be adjusted according to the actual situation without any restrictions.

[0056] The preset neighborhood window is a 3×3 window. Based on the Horn algorithm, the elevation values ​​of the grid cells within the preset neighborhood window in the digital elevation model are analyzed for slope to determine the terrain slope. The calculation of terrain slope using the Horn algorithm is a well-known technique in this field and will not be elaborated further.

[0057] It should be noted that, based on the actual terrain analysis, the terrain slope radii in the scene are within the range of [0, π / 2], therefore, the sine value can be calculated.

[0058] Among them, the preset safety constant term is a safety factor to prevent problems in subsequent ratio calculations caused by the sine value of the terrain slope arc value being 0. Optionally, the preset safety constant term is 0.01. Thus, the sum of the sine value of the terrain slope arc value and 0.01 is directly used as the slope influence index, and the ratio of the environmental inhibition index to the slope influence index is calculated as the migration resistance coefficient.

[0059] When the environmental inhibition index is 0, the numerator can also be directly set to 0.01 to facilitate the calculation of the migration resistance coefficient. It should also be noted that the migration resistance coefficient is truncated to its maximum value, such as 100. If the calculated migration resistance coefficient is greater than 100, it is directly set to 100 to avoid numerical explosion.

[0060] The slope influence index characterizes the migration driving force provided by the gravity component. Using elevation values ​​to calculate slope accurately reflects surface undulations. Even in flat areas after depression filling, this value is small due to the gentleness of the original terrain, leading to an increased migration resistance coefficient. This aligns with the physical law that flat depressions are prone to deposition. The environmental inhibition index characterizes the inhibitory effect of vegetation and clay on the migration process. Therefore, the calculated migration resistance coefficient characterizes the environmental and topographical resistance to the migration of surface materials, which will be analyzed in subsequent sections.

[0061] Driven by gravity, the migration of surface materials exhibits an irreversible unidirectional nature. To avoid circular dependencies or computational deadlocks when calculating the cumulative index across the entire region, a linear calculation logic that strictly follows the hydrological flow direction from upstream to downstream must be established.

[0062] Furthermore, in some embodiments of the present invention, based on hydrological flow direction analysis, the migration resistance coefficient and source emission intensity index of the grid cells are superimposed to obtain the deposition accumulation index of each grid cell, including: sorting the grid cells according to the flow direction based on hydrological flow direction analysis to obtain a hydrological topological sequence along the hydrological flow direction; and superimposing the migration resistance coefficient and source emission intensity index of the grid cells according to the order of the hydrological topological sequence to obtain the deposition accumulation index of each grid cell.

[0063] In this embodiment of the invention, sink filling can be performed on the digital elevation model to eliminate non-realistic depressions and ensure the connectivity of water flow paths. Then, the D8 unidirectional flow algorithm (Deterministic 8-node) is used to calculate the flow direction of each grid cell in the region to be analyzed. For any grid cell... Find the area with an elevation lower than its 8 neighborhoods (up, down, left, right, and diagonal directions). And the grid cell with the largest slope drop ,definition for The only receiving node, denoted as If there is no downstream node that meets the condition (i.e., all neighboring nodes are not lower than the current node), then mark it. It is the outflow boundary unit (Sink) of the watershed.

[0064] Based on the established point-to-point flow relationships, a global directed acyclic graph (DAG) is constructed. In this graph, each grid cell is a node, and the flow relationships are represented by directed edges. Subsequently, a topological sorting algorithm (such as Kahn's algorithm or depth-first search) is applied to sort the DAG, generating a unidirectional hydrological topological sequence.

[0065] Hydrological topological sequences possess the following key topological properties: for any element raster cell in the sequence All its upstream confluence units (i.e., all flows ultimately point to) The position index of the raster cell in the sequence is strictly less than 1. This sorting ensures that when subsequent steps traverse and calculate in this sequence, the current grid cell has fully received and integrated all input information from its entire upstream catchment area before being processed, thus guaranteeing the logical correctness of the cumulative calculation.

[0066] Further, in some embodiments of the present invention, the migration resistance coefficient and source emission intensity index of the grid cells are superimposed according to the hydrological topological sequence to obtain the depositional accumulation index of each grid cell, including: using the source emission intensity value as the initial excitation index of the grid cell; according to the hydrological topological sequence, using the migration resistance coefficient of the first grid cell in the sequence as the index term, and using the initial excitation index as the analytical excitation index of the first grid cell, performing exponential saturation decay calculation to obtain the depositional accumulation index of the first grid cell; calculating the difference between the initial excitation index and the depositional accumulation index as the downstream transport index; using the sum of the initial excitation index of the next grid cell in the sequence and the downstream transport index of all previous grid cells flowing to the next grid cell as the analytical excitation index of the next grid cell; and continuing to pass this downstream to calculate the depositional accumulation index of each grid cell.

[0067] The initial excitation index is the index used for initial hydrological deposition in each cell. Then, the index is saturated and decayed according to the order of the hydrological topology sequence. The specific calculation formula includes: In the formula, This is a preset deposition rate adjustment coefficient (default value is 1.0, but it can also be optimized using the aforementioned grid search method). This represents the deposition accumulation index of the i-th grid cell. This represents the migration resistance coefficient of the i-th grid cell. This represents the analytical excitation index of the i-th grid cell (or the initial excitation index if i=1). This formula simulates the saturation effect of feature rejection: as the migration resistance coefficient... As input increases, the retention rate gradually approaches 100%, but the total amount retained will never exceed the total amount input.

[0068] Then, the downstream transmission index that was not intercepted and continued to be transmitted downstream was calculated. : .

[0069] Retrieve the current raster cell The only downstream receiving node If node Existence (i.e.) If the boundary is not an outflow boundary, then determine all unidirectional grid cells flowing towards the (i+1)th grid cell, and calculate the flow direction. All previous grid cells of the grid cell Sum value, as ,Will With the The sum of the initial excitation indices of the nth grid cell is used as the value of the nth... The analysis excitation index of the first grid cell. And further analysis based on the first... The analysis excitation index of each grid cell is used to calculate the exponential saturation decay and subsequent calculations until the outflow boundary is reached, where there is no next receiving node, thus obtaining the depositional accumulation index for each grid cell. The depositional accumulation index aggregates the spatial context information of grid cell i and all its upstream catchments, and can effectively indicate the secondary enrichment trend of heavy metals.

[0070] By constructing a rigorous hydrological topological sequence and implementing an exponential saturation decay transfer mechanism based on the migration resistance coefficient on this sequence, the physical process of heavy metals being released from the source, decaying along the path, and accumulating in the downstream under the drive of surface runoff was effectively simulated. This method overcomes the misjudgment problem of "different quantities in the same spectrum" caused by the neglect of medium resistance and material conservation in the migration path in traditional inversion models. It makes the sedimentation accumulation index not only reflect the topographic confluence trend, but also integrate the dynamic regulation of heavy metal transport by vegetation blockage and clay adsorption, thus providing process constraint characteristics with clear hydrogeochemical significance for hyperspectral inversion.

[0071] S104: Combine vegetation cover coefficient, clay adsorption coefficient, mineral abundance coefficient and sedimentation accumulation index into the preset regression algorithm for nonlinear mapping analysis, construct an inversion model and determine the soil heavy metal content.

[0072] It should be noted that, in this embodiment of the invention, due to the cumulative nature of the confluence process, the values ​​in the downstream channel or confluence center region may be much larger than those in the source region (potentially differing by several orders of magnitude). To eliminate the negative impact of such extreme magnitude differences on the weights of the subsequent machine learning model, a logarithmic transformation can be performed on the deposition accumulation index of each grid cell for standardization before outputting the final features: The final output is the standardized sedimentation accumulation index.

[0073] In this embodiment of the invention, vegetation cover coefficient, clay adsorption coefficient, mineral abundance coefficient, and standardized sedimentary accumulation index can be used as four dimensional components to construct a four-dimensional feature vector. In this feature system, It played a crucial "process instruction" role: when A higher value indicates that the region is on a strong confluence path, and the heavy metal content should be higher than the value inferred solely from the spectral background; when When the value approaches 0, it indicates that the area is a source erosion zone or transport channel, and the heavy metal content is mainly controlled by the primary geological background.

[0074] The default regression algorithm is the Random Forest Regression algorithm. Because the contribution of surface medium characteristics and migration / accumulation processes to heavy metal content involves complex nonlinear interactions, and the controlling factors differ across different geomorphic units (e.g., ridges are primarily controlled by clay adsorption coefficients, while valleys are significantly affected by depositional accumulation indices), this embodiment uses the Random Forest Regression algorithm, which has good adaptability to nonlinear relationships, for modeling. Of course, those skilled in the art can also use other machine learning algorithms such as Support Vector Regression (SVR), XGBoost, or deep neural networks.

[0075] For the specific random forest regression algorithm, a training dataset was first constructed. Several ground sampling points were set up in the study area to collect surface soil samples and conduct laboratory chemical analysis to obtain measured values ​​of heavy metal (such as lead, cadmium, arsenic, etc.) concentrations. As ground truth labels, extract the joint feature vector at the coordinates of each sampling point. . Pair all data from sampling point m Divide into training set and validation set (e.g., 7:3 ratio).

[0076] Define the hyperparameters of the random forest, including the number of decision trees (e.g., 500), maximum depth, and splitting criterion. Supervised training of the model is performed using the training set to build a joint feature vector. Measured values ​​of heavy metal concentration The nonlinear mapping function.

[0077] Next, model training and feature importance evaluation are performed. During training, the importance of each input feature is evaluated using out-of-bag error (OOBError) or the Gini index. If the results show... The presence of significant importance weights further validates the effectiveness of introducing spatial process features.

[0078] Then, a global inversion and deduction is performed, and the joint feature vectors of all raster cells in the area to be analyzed are obtained. Input the trained regression model to calculate the predicted value for each unit in the entire domain. The final result is a spatially continuous map of heavy metal concentration distribution. This outcome combines the ability of hyperspectral imaging to precisely identify surface materials with the ability of flux models to macroscopically constrain hydrological processes, achieving accurate spatial quantification of soil heavy metal content.

[0079] This invention effectively overcomes the non-stationarity of the "spectrum-content" mapping caused by neglecting the hydraulic migration process in existing inversion methods by introducing a sedimentary accumulation index as a key process constraint feature. This index constructs hydrological flow direction based on a digital elevation model and integrates surface migration resistance characterized by vegetation cover coefficient and clay adsorption coefficient with source emission intensity determined by topographic height index and mineral abundance coefficient. Thus, it quantifies the secondary enrichment trend of heavy metals driven by runoff without the need for complex hydrological simulation. By inputting this dynamic accumulation index and static surface medium parameters into the regression algorithm, the model can distinguish areas with similar spectral characteristics but different hydrological locations (such as erosion zones and sedimentary zones), significantly improving the reliability and spatial accuracy of soil heavy metal content inversion in areas with large topographic relief.

[0080] On the other hand, the present invention also provides a soil heavy metal content inversion modeling system based on hyperspectral imaging. The system includes a memory, a processor, and a computer program stored in the memory and running on the processor. When the processor executes the computer program, it implements the steps of the method described in any of the foregoing descriptions.

[0081] It should be noted that the order of the above embodiments of the present invention is merely for descriptive purposes and does not represent the superiority or inferiority of the embodiments. The processes depicted in the accompanying drawings do not necessarily require a specific or sequential order to achieve the desired result. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0082] The various embodiments in this specification are described in a progressive manner. The same or similar parts between the various embodiments can be referred to each other. Each embodiment focuses on describing the differences from other embodiments.

Claims

1. A method for inverting and modeling soil heavy metal content based on hyperspectral imaging, characterized in that, The method includes: Within different raster cells of the region to be analyzed, acquire spatially resolution registered hyperspectral images and digital elevation models, and determine the elevation value and spectral reflectance map of each raster cell. Based on the spectral reflectance distribution of each grid cell in different bands of the spectral reflectance map, the vegetation cover coefficient, clay adsorption coefficient, and mineral abundance coefficient are determined; based on the numerical difference in elevation between each grid cell and its adjacent grid cells, the topographic height index of the grid cell is determined; and combined with the topographic height index and mineral abundance coefficient, the source emission intensity index is determined. Based on the topographic slope, vegetation cover coefficient, and clay adsorption coefficient of the grid cells obtained from the analysis of the digital elevation model, the migration resistance coefficient of the grid cells is determined; according to the hydrological flow direction analysis, the migration resistance coefficient and the source emission intensity index of the grid cells are superimposed to obtain the deposition accumulation index of each grid cell. By combining vegetation cover coefficient, clay adsorption coefficient, mineral abundance coefficient and sedimentation accumulation index into a preset regression algorithm for nonlinear mapping analysis, an inversion model is constructed to determine the heavy metal content in the soil.

2. The method for inverting and modeling soil heavy metal content based on hyperspectral imaging as described in claim 1, characterized in that, The determination of vegetation cover coefficient, clay adsorption coefficient, and mineral abundance coefficient based on the spectral reflectance distribution of each grid cell in different bands of the spectral reflectance map includes: The vegetation cover coefficient is determined based on the difference in spectral reflectance of the preset vegetation-affected bands in the spectral reflectance diagram. The clay adsorption coefficient is determined based on the absorption depth of the pre-defined clay characteristic influence band in the spectral reflectance diagram. The mineral abundance coefficient is determined based on the integral area of ​​the absorption valleys of the pre-defined mineral characteristic influence bands in the spectral reflectance diagram.

3. The method for inverting and modeling soil heavy metal content based on hyperspectral imaging as described in claim 2, characterized in that, The preset vegetation influence bands include the near-infrared band and the red band, with a center wavelength of 670 nm in the red band and 860 nm in the near-infrared band; the preset clay characteristic influence bands include the Al-OH characteristic band with a wavelength of 2200 nm; the preset mineral characteristic influence bands include the mineral absorption band in the 800 nm-1000 nm band; the methods for obtaining the vegetation cover coefficient, clay adsorption coefficient, and mineral abundance coefficient include: Calculate the reflectance difference between the center wavelength of the near-infrared band and the center wavelength of the red band; The reflectance difference was normalized and used as the vegetation cover coefficient. The absorption depth of the trough in the preset clay characteristic influence band is obtained and normalized as the clay adsorption coefficient. In a spectral reflectance graph with wavelength as the abscissa and reflectance as the ordinate, the spectral curves of the mineral absorption bands and the area formed by the spectral curves and the continuum are obtained based on definite integrals, and the areas are normalized to obtain the mineral abundance coefficients.

4. The method for inverting and modeling soil heavy metal content based on hyperspectral imaging as described in claim 1, characterized in that, Based on the numerical difference in elevation values ​​between each raster cell and its adjacent raster cells, the terrain height index of the raster cell is determined, including: In the digital elevation model, the mean of all elevation values ​​of each grid cell within a preset neighborhood window is calculated and used as the neighborhood mean. The difference between the elevation value of a grid cell and the mean of its neighborhood is used as the terrain height index.

5. The method for inverting and modeling soil heavy metal content based on hyperspectral imaging as described in claim 1, characterized in that, By combining the topographic height index and mineral abundance coefficient, the source emission intensity index is determined, including: Anomaly detection was performed on the mineral abundance coefficients of all raster cells using the three-standard-deviation principle, with the sum of the mean and three standard deviations used as the anomaly detection threshold. If the mineral abundance coefficient of any grid cell is greater than the anomaly judgment threshold and the terrain height index is greater than 0, then the source emission intensity index of the grid cell is set to the mineral abundance coefficient. Otherwise, the source emission intensity index of the grid cell is set to a preset background constant.

6. The method for inverting and modeling soil heavy metal content based on hyperspectral imaging as described in claim 1, characterized in that, Based on the terrain slope, vegetation cover coefficient, and clay adsorption coefficient of the raster cells obtained from the analysis of the digital elevation model, the migration resistance coefficient of the raster cells is determined, including: The environmental inhibition index is obtained by weighted summation of the vegetation cover coefficient and the clay adsorption coefficient. Based on the Horn algorithm, slope analysis is performed on the elevation values ​​of grid cells within a preset neighborhood window in the digital elevation model to determine the terrain slope; after converting the terrain slope to radians, the sum of the sine value and the preset safety constant term is calculated as the slope influence index. The ratio of the environmental inhibition index to the slope influence index is calculated and used as the migration resistance coefficient.

7. The method for inverting and modeling soil heavy metal content based on hyperspectral imaging as described in claim 1, characterized in that, Based on hydrological flow direction analysis, the migration resistance coefficient and source emission intensity index of each grid cell are superimposed to obtain the deposition accumulation index for each grid cell, including: Based on the hydrological flow direction analysis, the grid cells are sorted according to the flow direction to obtain the hydrological topology sequence along the hydrological flow direction; The migration resistance coefficient and source emission intensity index of each grid cell are superimposed according to the hydrological topological sequence to obtain the deposition accumulation index of each grid cell.

8. The method for inverting and modeling soil heavy metal content based on hyperspectral imaging as described in claim 7, characterized in that, The migration resistance coefficient and source emission intensity index of each grid cell are superimposed according to the hydrological topological sequence to obtain the deposition accumulation index of each grid cell, including: The source emission intensity is used as the initial excitation index for the grid cell; Following the order of the hydrological topology sequence, the migration resistance coefficient of the first grid cell in the sequence is used as the exponential term, and the initial excitation index is used as the analytical excitation index of the first grid cell to calculate the exponential saturation decay, thus obtaining the deposition accumulation index of the first grid cell. The difference between the initial excitation index and the depositional accumulation index is calculated as the downstream transport index. The sum of the initial excitation index of the next grid cell in the sequence and the downstream transport index of all previous grid cells flowing to the next grid cell is used as the analytical excitation index of the next grid cell. This process is continued downstream to calculate the depositional accumulation index of each grid cell.

9. The method for inverting and modeling soil heavy metal content based on hyperspectral imaging as described in claim 1, characterized in that, The preset regression algorithm is the random forest regression algorithm.

10. A soil heavy metal content inversion modeling system based on hyperspectral imaging, the system comprising a memory, a processor, and a computer program stored in the memory and running on the processor, characterized in that, When the processor executes the computer program, it implements the steps of the method as described in any one of claims 1 to 9.

Images