A method for soil risk assessment for upland rice fields

Abstract

The application discloses a kind of soil risk assessment methods for red paddy field, specifically related to soil risk assessment field, for solving the problem that existing method cannot accurately identify the anoxic risk below tillage layer and its influence on root system;By constructing the combined monitoring mechanism of heat flux and oxygen concentration, the metabolic heat release power and oxygen consumption trend are extracted, the reduction dominant area below the tillage layer is identified, and the spatial boundary is fitted, and the three-dimensional distribution map of the anoxic dominant area is established;Combined with the distribution of crop root system, the spatial boundary of the overlap of root zone and anoxic area is extracted, and the profile expansion rate and root zone change characteristics are analyzed in the continuous monitoring period, and the soil degradation trend identification parameter set is constructed.Finally, according to the preset degradation judgment weight, the parameter sequence is scored and mapped, and the soil degradation risk grade atlas covering the entire red paddy field area is generated, with high spatial resolution and high discrimination accuracy of soil risk expression ability.

Description

Technical Field

[0001] This invention relates to the field of soil risk assessment technology, and more specifically, to a soil risk assessment method for red soil paddy fields. Background Technology

[0002] Due to the poor parent material, high acidity, and rapid decomposition of organic matter in red soil paddy fields, long-term cultivation often results in a stratified reduction zone within 20cm to 60cm below the topsoil. This zone is characterized by decreased oxygen permeability, altered microbial community structure, and enhanced activity and transformation of elements such as iron and manganese. Under continuous flooding or irrigation conditions, this reduction zone easily evolves into a stable deep anoxic zone, where a strong reduction process occurs, with iron and manganese acting as electron acceptors, producing large amounts of Fe. 2+ Mn 2 + This process, accompanied by the decomposition of organic reducing substances producing gases such as CH4 and H2S, hinders root growth in the depth direction, shrinks the effective root zone of crops, and induces nitrogen denitrification loss and heavy metal activation and migration, severely affecting the nutrient retention capacity and ecological stability of red soil paddy fields. However, these hypoxic areas often exhibit concealed, discontinuous, and highly heterogeneous spatial distribution. Their presence is often not accompanied by significant changes in conventional soil indicators such as surface pH, electrical conductivity, and organic matter content, making it difficult to identify their dynamic evolution process using conventional sampling and surface monitoring methods. Especially under the background of high temperature and humidity and long-term high nitrogen application, deep soil is more prone to oxygen depletion and strong reducing processes, inducing reducing stress and rhizosphere microenvironment disorder, exacerbating root distribution chaos and microecological imbalance, and forming a degraded soil structure state with a continuous evolutionary trend.

[0003] Therefore, it is necessary to construct a method to accurately reveal the interaction characteristics between deep hypoxia zones and crop root zones, establish a stable and continuous risk perception framework, and construct a quantifiable soil degradation trend identification parameter system in order to achieve a scientific assessment of the risk of red soil paddy fields. Summary of the Invention

[0004] In order to overcome the above-mentioned defects of the prior art, embodiments of the present invention provide a soil risk assessment method for red soil paddy fields to solve the problems mentioned in the background art.

[0005] To achieve the above objectives, the present invention provides the following technical solution:

[0006] A soil risk assessment method for red soil paddy fields includes the following steps:

[0007] S1. Collect the heat flux and oxygen concentration of the soil at pre-set sampling points in the red soil paddy field below the topsoil layer;

[0008] S2. Calculate the unit volume metabolic heat release power curve at different depths of the sampling point, and identify the reduction-dominant region by combining the time change trend of oxygen concentration.

[0009] S3. Spatial fitting is performed on the boundary of the reduction-dominant region to construct a three-dimensional distribution map of the deep hypoxia-dominant region.

[0010] S4. Construct a spatial distribution map of crop roots in red soil paddy fields, spatially register it with the three-dimensional distribution map of the hypoxia-dominant area, and extract the boundary contours of the overlapping areas.

[0011] S5. Calculate the boundary contour expansion rate and root density change amplitude of the overlapping area to construct a set of parameters for soil degradation trend identification.

[0012] S6. Based on the index sequence output by the degradation trend identification parameter set, calculate the degradation score of red soil paddy field and generate a soil risk level map.

[0013] In a preferred embodiment, step S1, specifically collecting the heat flux and oxygen concentration of the soil at predetermined sampling points in the red soil paddy field below the topsoil layer, includes:

[0014] Multiple evenly distributed sampling points were pre-set in red soil paddy fields to establish a planar coordinate system;

[0015] Below the topsoil layer, multiple depth layers were set in the longitudinal profile of the sampling points in the red soil paddy field. Temperature and oxygen concentration were collected at a fixed sampling frequency, and data from different depth layers were recorded.

[0016] All depth layer data are arranged according to acquisition time to form a parameter set with a unified time base and depth label;

[0017] Vertical heat flux is derived from temperature sequences at different depths at the same time point using spatial difference, and a heat flux distribution array is constructed with time as the horizontal axis and depth as the vertical axis.

[0018] In a preferred embodiment, the derivation of vertical heat flux from temperature sequences at different depths at the same time point using spatial difference specifically includes:

[0019] Perform differential calculations on the temperature sample values ​​of any two consecutive depth points in a set profile to obtain a vertical temperature gradient sequence;

[0020] Obtain the thermal diffusivity and volumetric heat capacity parameters corresponding to the red soil type, and convert the temperature gradient into heat flux value through the heat transfer calculation formula;

[0021] All heat flux values ​​are aggregated to construct a heat flux profile. The vertical structure of the heat flux is extracted at each time step, and a standardized heat flux distribution array is constructed.

[0022] In a preferred embodiment, step S2 involves calculating the unit volume metabolic heat release power curves at different depths of the sampling point, and, in conjunction with the time-varying trend of oxygen concentration, identifying the reduction-dominant region. This specifically includes:

[0023] Perform central differential calculations on the constructed heat flux distribution array in chronological order to obtain the rate of change of heat flux in sub-segments between a set number of consecutive adjacent depth layers;

[0024] The rate of change is the vertical heat flux gradient of each sub-section.

[0025] The rate of change was converted into the metabolic heat release power corresponding to the sub-segment volume, and the time evolution curve of the metabolic heat release power of the sub-segment was established.

[0026] In the parameter set, the time period corresponding to the heat flux change of each sub-segment is obtained, and the concentration change rate is extracted from the oxygen concentration sequence within the same time period to generate the oxygen consumption trend curve of the sub-segment.

[0027] Align the two types of curves with a uniform time resolution and calculate the covariance coefficient to identify sub-segments with increased metabolic heat release power and enhanced oxygen consumption, and mark them as reduction-dominant regions.

[0028] The covariance coefficients are numerical expressions after covariance normalization.

[0029] In a preferred embodiment, when identifying sub-segments with increased metabolic heat release power and enhanced oxygen consumption, if multiple sub-segments meet the criteria, the sub-segment with the largest covariance coefficient is selected as the reduction-dominant region.

[0030] In a preferred embodiment, step S3, which involves spatially fitting the boundary of the reduction-dominant region to construct a three-dimensional distribution map of the deep hypoxia-dominant region, specifically includes:

[0031] Based on the identified dominant restoration region in each sampling point, the upper boundary depth value and lower boundary depth value of the dominant restoration region are extracted and used as the upper boundary point set and lower boundary point set, respectively.

[0032] All boundary points are aligned according to the planar coordinate system of their sampling points, and spatial fitting is performed on each point to obtain the upper and lower boundary spatial fitting surfaces covering the dominant region of the entire sampling point restoration.

[0033] The stable spatial envelope region between the two fitted surfaces is used as a three-dimensional distribution map dominated by hypoxia.

[0034] In a preferred embodiment, step S4 involves constructing a spatial distribution map of crop roots in red soil paddy fields, spatially registering it with a three-dimensional distribution map of the hypoxia-dominant area, and extracting the boundary contours of overlapping areas. Specifically, this includes:

[0035] The root depth of crops at each sampling point is collected, and the root zone distribution at different depths is constructed with the root location of the crop at each sampling point as a reference.

[0036] The three-dimensional distribution map of the hypoxia-dominant area was retrieved, and the coordinates were aligned with the root zone distribution at different depths using an affine transformation algorithm. The spatial intersection of the overlapping areas was calculated and the contours were extracted.

[0037] In a preferred embodiment, step S5, calculating the boundary contour expansion rate and root density change magnitude of the overlapping area to construct a soil degradation trend identification parameter set specifically includes:

[0038] Within a preset continuous time period, differential analysis is performed on the boundary contour of the overlapping area in the time dimension to obtain the contour expansion rate.

[0039] Simultaneously analyze the variation of root density over time in the root zone distribution of all sampling points in the overlapping region;

[0040] The two types of indicators, namely the rate of outward expansion of the contour and the magnitude of change in root density, were organized into a set of parameters for identifying soil degradation trends according to the location of the sampling point and the time period label.

[0041] In a preferred embodiment, step S6, calculating the degradation score of red soil paddy field soil and generating a soil risk level map based on the index sequence output by the degradation trend identification parameter set, specifically includes:

[0042] The time series of each indicator in the degradation trend identification parameter set is standardized, and a multi-indicator weighted score is performed based on the preset degradation judgment weight.

[0043] A preset soil degradation risk level mapping range is used to map the scoring results to different soil degradation risk levels. In the overlapping area, the sampling points corresponding to all risk levels are mapped to the actual spatial location of the red soil paddy field, generating a soil degradation risk level map covering the entire red soil paddy field area.

[0044] The technical effects and advantages of the soil risk assessment method for red soil paddy fields of this invention are as follows:

[0045] This invention overcomes the limitations of traditional soil testing methods, which rely on surface physicochemical parameters, suffer from poor spatial continuity, and lack dynamic identification capabilities. Based on heat flux and oxygen concentration data below the topsoil layer, it achieves for the first time accurate identification and metabolic intensity assessment of deep reduction-dominant regions, improving the depth and accuracy of identifying hypoxic stress areas. By constructing a spatial registration layer between hypoxic-dominant regions and root distribution, an interactive expression mechanism between soil physical state and biologically active zones is established, effectively solving the problem of difficulty in quantifying root zone interference areas. Furthermore, by introducing contour change rate and root density evolution indicators in the time dimension, continuous monitoring and dynamic modeling of soil degradation trends are achieved. Finally, a soil risk level map is generated through multi-indicator sequence modeling, making the spatial expression of soil quality status more objective and operable. This method has the advantages of being quantifiable, visual, and traceable, providing high-precision, multi-dimensional technical support for monitoring and regional regulation of red soil degradation in paddy fields. Attached Figure Description

[0046] Figure 1 This is a schematic diagram of a soil risk assessment method for red soil paddy fields according to the present invention. Detailed Implementation

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

[0048] Example 1

[0049] Figure 1 This invention provides a soil risk assessment method for red soil paddy fields, which includes the following steps:

[0050] S1. Collect the heat flux and oxygen concentration of the soil at pre-set sampling points in the red soil paddy field below the topsoil layer;

[0051] S2. Calculate the unit volume metabolic heat release power curve at different depths of the sampling point, and identify the reduction-dominant region by combining the time change trend of oxygen concentration.

[0052] S3. Spatial fitting is performed on the boundary of the reduction-dominant region to construct a three-dimensional distribution map of the deep hypoxia-dominant region.

[0053] S4. Construct a spatial distribution map of crop roots in red soil paddy fields, spatially register it with the three-dimensional distribution map of the hypoxia-dominant area, and extract the boundary contours of the overlapping areas.

[0054] S5. Calculate the boundary contour expansion rate and root density change amplitude of the overlapping area to construct a set of parameters for soil degradation trend identification.

[0055] S6. Based on the index sequence output by the degradation trend identification parameter set, calculate the degradation score of red soil paddy field and generate a soil risk level map.

[0056] In step S1, the heat flux and oxygen concentration of the soil at preset sampling points in the red soil paddy field are collected below the cultivated layer.

[0057] Using the actual area of ​​the red soil paddy field as the spatial framework, and considering factors such as field area, shape, irrigation and drainage patterns, and historical agronomic practices, sampling points were laid out using an equidistant spacing principle. Specifically, a standard sampling grid was established within the entire field area using a 10m × 10m grid, with one sampling point placed at the center of each grid unit. For red soil paddy fields larger than 1 hectare, the grid density could be appropriately increased to a 20m interval to ensure a balance between representativeness and work efficiency. When laying out the sampling points, non-cultivated areas such as field ridges and drainage ditches were avoided, and sufficient lateral comparability between sampling points was ensured. To further establish a planar coordinate system, spatial positioning was performed on all sampling points using a positioning instrument for high-precision coordinate acquisition, with an accuracy controlled within ±2 cm. A two-dimensional coordinate system for the field was constructed using a corner point as the origin, providing a comparable spatial basis for subsequent sampling data. During this process, the layout of the sampling points was appropriately adjusted based on information such as topographic relief, irrigation impact range, and historical fertilization intensity within the field to ensure balanced representativeness and regional coverage of the acquired data. After the sampling points are set up, each sampling point is numbered and coded, and a sampling point distribution map is generated to provide a spatial index reference for subsequent sampling operations and data analysis.

[0058] Based on the established lateral sampling point coordinate system, vertical soil profile drilling was conducted at each sampling point. The profile depth was primarily set from below the topsoil layer to 60 cm. Considering that the reducing zone in red soil paddy fields typically extends to a depth of 20 cm to 60 cm, the default measurement depth levels were set to five fixed levels: 20 cm, 30 cm, 40 cm, 50 cm, and 60 cm. High-precision soil temperature sensors (with a sensing accuracy of at least 0.1℃) and oxygen concentration sensors were installed at each depth layer to ensure that the sensor probes could directly contact the representative soil of that layer, avoiding pore interference. The sensor sampling frequency was uniformly set to once every 5 minutes to ensure the capture of sufficiently detailed temporal dynamic characteristics, and all sensors were started synchronously to avoid time drift. The sampling operation was recommended to continue for at least 10 days to cover the diurnal variation characteristics during the crop root activity period and the response process to key meteorological events such as high temperatures and rainfall. All sensor data was recorded in real time using timestamps as indexes and stored in a unified format time-series file, with a data index established according to the sampling point number and depth layer label. During the data integration phase, data from different sampling points and depth layers are vertically merged in chronological order to ensure that a dataset of all sampling points across the entire profile can be obtained at each time point. Simultaneously, a parametric dataset structure based on a three-dimensional coordinate system of "time-depth-sampling point" is established for subsequent computational processing, forming a high spatiotemporal resolution monitoring dataset with a unified time reference and depth layer labels. This dataset not only provides fundamental support for heat flux and oxygen consumption trend analysis but can also be directly used in various computational tasks such as spatial interpolation, heat conduction modeling, and dynamic simulation of oxygen diffusion.

[0059] After acquiring temperature sampling data, temperature gradient calculations are performed to derive the heat flux distribution structure. A longitudinal profile depth level is selected in the red soil paddy field, assuming a 10cm interval between each sampling point, starting from 20cm below the tillage layer and ending at 60cm, forming five consecutive depth layers. For any two consecutive depth layers, such as the 20cm and 30cm layers, their temperature values ​​at the same time point are extracted, and a direct difference operation is performed to calculate the corresponding vertical temperature gradient. This difference operation is performed completely once in each time slice to form a temperature gradient sequence for different sub-segments. The temperature gradient value of each sub-segment represents the degree of temperature change per unit depth within that segment. After completing the temperature gradient sequence, to ensure the physical accuracy of the heat flux derivation, thermophysical parameters corresponding to the characteristics of red soil are introduced. Due to its low organic matter content and loose structure, the thermal diffusivity of red soil is typically in the range of (0.6-1.1)×10⁻⁶. - 6 m 2 The value is between / s, and the specific value can be determined based on experimental measurements. For red soil samples, use 0.85×10-6 m 2 / s is the set value in this embodiment. Meanwhile, to ensure the rationality of the heat capacity calculation, a volumetric heat capacity parameter needs to be introduced. The volumetric heat capacity of red soil often varies depending on the moisture content; in this implementation scenario, its value is set to 2.5 × 10⁻⁶. 6 J / (m 3 (·K). After clarifying the above two parameters, the temperature gradient is combined with the thermal diffusivity and volumetric heat capacity, and the heat flux per unit time per unit area for each sub-segment is derived using the one-dimensional heat conduction calculation formula. This type of calculation does not involve multidimensional modeling, but only focuses on the vertical profile dimension. The calculation results reflect the distribution of heat flux among sub-segments at different depths.

[0060] After the heat flux data is calculated, the heat flux of all cross-sectional structures needs to be summarized according to time points to construct a heat flux profile map under a continuous time series. The heat flux profile of each time slice should cover all depth sub-segments from 20cm to 60cm, forming a complete vertical heat flux structure representation. The horizontal axis of the profile map represents the time series, and the vertical axis represents the vertical depth position. The color or grayscale value corresponding to each data point is the magnitude of the heat flux. In this way, a heat flux distribution array can be constructed. This distribution array, as a two-dimensional spatial-temporal matrix, intuitively reflects the vertical distribution of heat flux at different time points.

[0061] In step S2, the metabolic heat release power curve per unit volume at different depths of the sampling point is calculated, and the reduction-dominant region is identified by combining the time change trend of oxygen concentration.

[0062] The established heat flux distribution array is processed chronologically to extract vertical heat flux profiles corresponding to different times. At each time slice, the heat flux values ​​are reorganized into a one-dimensional vector according to depth layers, and a central difference calculation is performed on the heat flux values ​​between adjacent depth layers. Specifically, a continuous depth layer with a preset depth difference (default not less than 10cm) is selected, for example, the heat flux value between 20cm and 30cm. The difference is divided by the preset length to form the heat flux gradient of that sub-segment. To ensure the stability and physical rationality of the calculation, all difference operations adopt the central difference mode, meaning the gradient at the current depth point is jointly determined by the heat flux values ​​of the layers above and below it. Using this method, multiple continuous sub-segment heat flux gradient sequences are obtained at each time point, forming a complete set of vertical change rate data structures. Since the thermal conductivity of the red soil profile may differ at different depths, it is recommended that the depth difference of the sub-segments be kept consistent during actual execution to improve the effectiveness of the comparison. All heat flux change rates are stored as intermediate variables for subsequent metabolic heat power conversion, and are managed by bidirectional indexing using depth labels and time labels to ensure that each data point can be uniquely located and retrieved.

[0063] The heat flux change rate sequence obtained above is converted into the metabolic heat release power per unit volume of each sub-section to more accurately reflect the local microbial metabolic intensity of the soil. First, the actual physical volume of each sub-section is determined. The volume calculation is based on a cylindrical model constructed using the sub-section thickness (e.g., 10 cm) and the lateral sampling influence radius of the sampling point (e.g., 0.5 m). The unit volume of each sub-section is calculated using standard geometry. Then, the heat flux change rate per unit area is divided by the volume of that sub-section to obtain the heat release rate per unit volume. This heat release rate is considered the metabolic heat release power, with units of W / m². 3 To ensure the accuracy of the conversion, known soil thermophysical parameters of the red soil paddy field profile, such as volumetric heat capacity and thermal diffusivity (based on the values ​​set in step S1), need to be considered. After the conversion, the metabolic heat release power of each sub-segment is arranged as a continuous function according to the time series, forming its time evolution curve. The curve is plotted with time on the horizontal axis and power value on the vertical axis, representing the change in heat release intensity of the sub-segment at each time point. Throughout the process, each curve maintains a time resolution consistent with the sampling frequency to facilitate subsequent co-analysis with the oxygen concentration change sequence. This processing clearly reveals the dynamic changes of the thermal metabolic process in the soil profile, providing quantitative indicators to support the determination of reduction behavior.

[0064] Oxygen concentration samples for the corresponding time period are retrieved from the existing parameter set. First, the heat flux change time period for each sub-segment is identified (the specific range can be controlled by setting the minimum temperature sensing accuracy). This time period is then aligned with the oxygen concentration sampling sequence using time labels. Differential calculations are performed on the oxygen concentration sequences within the same time interval to obtain the rate of decrease in oxygen concentration per unit time, which is taken as the oxygen consumption rate. To ensure data accuracy, the original oxygen concentration data should be filtered and smoothed before differential calculation. A sliding window averaging method is used to eliminate abrupt noise; the window size can be set to 3 to 5 data points depending on the sampling frequency. The differential method is backward differential to ensure that the rate of change reflects the actual consumption trend. The calculated rate of change is expressed as a negative value to indicate the decreasing trend of oxygen concentration, while the absolute value reflects the strength of the consumption rate. The oxygen change rate of each sub-segment over its time evolution is used as the vertical axis, and time as the horizontal axis to plot an oxygen consumption trend curve, forming a second type of trend expression corresponding to the metabolic heat release power curve. The aforementioned oxygen consumption trend curves must maintain the same time reference and resolution as the metabolic heat curves. All curves are organized and stored according to sub-segment labels to form a complete spatial-temporal oxygen consumption trend dataset.

[0065] For each sub-segment, two curves are used as a sample group. The value pairs at each unified time node are calculated to form a time series sample pair set. For each sample pair, the covariance coefficient of that sub-segment is calculated based on the normalized formula for covariance. The covariance coefficient is calculated as follows: first, the covariance of the metabolic heat release power and the rate of change of oxygen consumption are calculated; then, the standard deviations of the two series are calculated separately, and the covariance is divided by the product of the standard deviations to obtain a normalized correlation index. This coefficient is a real number between -1 and 1, where a value close to 1 indicates a strong consistency in the trends of the two types of curves, and is particularly suitable for identifying thermo-oxygen coupling segments where increased metabolic activity and simultaneously increased oxygen consumption occur. After calculating the covariance coefficient for all sub-segments, an identification threshold is set to screen segments that satisfy the thermo-oxygen synergistic trend of "increased metabolic heat release power and increased oxygen consumption". The identification threshold should be set based on experimental data and risk identification requirements. For example, when the covariance coefficient is greater than 0.5, and the metabolic heat release power and oxygen consumption both reach or exceed the average value of all depth layers at the sampling point, it can be considered to have significant thermo-oxygen synchronicity. Among the selected candidate segments, the covariance coefficient values ​​corresponding to each sub-segment are further compared, and the sub-segment with the largest covariance coefficient is selected as the reduction-dominant region for that sampling point. This selection method ensures that among all sub-segments with obvious thermo-oxygen coupling trends, the target region with the most significant coupling strength and the clearest indication is prioritized, effectively improving the representativeness and accuracy of reduction region identification.

[0066] In step S3, spatial fitting is performed on the boundary of the reduction-dominant region to construct a three-dimensional distribution map of the deep hypoxia-dominant region.

[0067] After identifying the dominant reduction region, the vertical profile data of each sampling point in the red soil paddy field were processed to extract the depths of the upper and lower boundaries of the corresponding dominant reduction region point by point. During this extraction process, the covariant condition of increased metabolic heat release power and enhanced oxygen consumption was used to determine whether each sub-segment was a dominant reduction region. Once a sub-segment was confirmed as a dominant reduction region, the top layer depth of that sub-segment was used as the upper boundary point, and the bottom layer depth as the lower boundary point. This allows for high-precision positioning of the boundary range of the dominant reduction region. In the specific recording process, the upper boundary depth value was added as one element to the upper boundary point set, and the lower boundary depth value was added as another element to the lower boundary point set. Both sets are structurally indexed by the lateral coordinates of the sampling points, and each record corresponds to a sampling point number, lateral position, and its corresponding boundary depth value.

[0068] After constructing the upper and lower boundary point sets, spatial registration is performed on all boundary points according to their corresponding sampling points in the red soil paddy field. Specifically, using a unified reference coordinate system within the paddy field area as a benchmark, the horizontal and vertical positions of all sampling points are mapped to the x and y axes in two-dimensional space, with depth values ​​used as the z-axis coordinates, forming data points in the form of (x, y, z) pairs. Spatial fitting operations are performed on both the upper and lower boundary point sets. The preferred fitting methods are spline surface fitting or bivariate smoothing interpolation to ensure high spatial continuity and fitting accuracy of the generated fitting surface. During spline fitting, the surface shape is controlled by setting node density and smoothing coefficients to avoid drastic fluctuations or fitting offsets caused by local outliers. Sampling points in non-dominant restoration regions are excluded during fitting. If abrupt changes or fitting errors exceed the preset tolerance range (e.g., fitting residual greater than 2 cm) appear in the spline surface fitting results, the boundary point set must be returned to perform local outlier removal and refitting to ensure the engineering applicability of the fitting results. In the final output, the fitting result corresponding to the upper boundary point set is defined as the upper boundary space surface, and the fitting result of the lower boundary point set is defined as the lower boundary space surface. The two surfaces are perfectly aligned in the coordinate dimension, which has the basic conditions for constructing a spatial closed envelope.

[0069] After obtaining the fitted surfaces of the upper and lower boundary spaces, an envelope construction operation is further performed. Specifically, the vertical depth difference between the two surfaces is extracted at a fixed resolution on the coordinate axes to form a continuous spatial vertical profile. At each location, the area defined by the upper and lower boundary surfaces is considered the effective space of the deep anoxic predominant zone on that cross-section. This process is then advanced point by point along the transverse axis, summarizing all profile envelope areas to form a spatially continuous three-dimensional data block, constituting a three-dimensional distribution map of the deep anoxic predominant zone in red soil paddy fields.

[0070] In step S4, a spatial distribution map of crop roots in red soil paddy fields is constructed, and the map is spatially registered with a three-dimensional distribution map of the hypoxia-dominant area to extract the boundary contours of overlapping areas.

[0071] Using a profile excavation and root stripping observation method, the soil was excavated from the surface down to the deepest point where the root system ceases to distribute. The vertical depth from the surface to the root base was measured and taken as the maximum root distribution depth at that sampling point. To ensure data representativeness and stability, a 50cm × 50cm observation area was extended along the edge of each sampling point, and at least three rice plants were selected within this area for averaging. The average value was taken as the final root depth at that sampling point. Within each depth layer, factors such as whether the roots penetrated the layer, the trend of the main root distribution density within the layer, and the distribution characteristics of the root type (taproot or fibrous root) were considered to determine whether the layer constituted an effective root zone. Depth layers with continuous root distribution were designated as root zone distribution segments, and each layer was coded and recorded according to its actual depth range.

[0072] Since the root zone layer and the restoration-dominant region are built on independent sampling data, their coordinate origins, direction vectors, and scale coefficients differ. Therefore, an affine transformation algorithm is required to unify their coordinate systems. Affine transformations include basic operations such as translation, rotation, scaling, and shearing, which can map the coordinate system of the root zone layer to the three-dimensional coordinate system of the anoxic region. Figure 1 A spatial reference frame is established. In practice, three sets of control points with known spatial correspondences (e.g., typical sampling point locations and corresponding depths) need to be selected for transformation parameter calculation, a transformation matrix is ​​constructed, and it is applied to the coordinates of all layers in the entire root band map group to ensure that the spatial position of the root band region at each depth layer is completely aligned with the dominant restoration region.

[0073] After coordinate alignment, the intersection of the root zone boundary and the anoxic dominant zone at each depth layer is calculated. The intersection region represents the overlapping portion of the root zone distribution and the anoxic region at that depth layer, indicating the spatial segment where root growth encounters anoxic stress. In the intersection calculation, spatial Boolean operations are used to find the intersection of polygons, resulting in closed overlapping contour lines. The boundary point sets and geometric attributes are then extracted to establish the corresponding spatial contour structure. To maintain the depth correspondence between contours, contour information is stored separately according to layer number and organized into a three-dimensional overlapping region contour set in a unified format.

[0074] In step S5, the expansion rate of the boundary contour of the overlapping area and the change in root density are calculated to construct a set of parameters for identifying soil degradation trends.

[0075] An observation period is set, for example, 60 consecutive days, with time steps divided into time periods of 10 days each. The boundary contour maps of overlapping areas generated within each time step are subjected to temporal differencing. Volume differences are calculated for the boundary contours of any two adjacent time periods, and an extrapolation vector field is constructed to determine the expansion rate of the contours in different spatial orientations. This expansion rate is measured as the average displacement of the contour boundary per unit time in the direction perpendicular to its local normal, in millimeters per hour or centimeters per day, reflecting the trend of hypoxia-dominant processes advancing towards the root zone over time.

[0076] Based on previously recorded root distribution depth information at each sampling point, root strip segments intersecting with overlapping areas are extracted according to predetermined depth levels, and root strip distribution maps for each time period are constructed, uniformly formatted as a three-dimensional grid. Subsequently, the difference in root strip occupancy ratio between two adjacent time periods for each grid cell is calculated to obtain the variation amplitude of root strip distribution within a local area. The variation amplitude is taken as the absolute percentage change in the root strip proportion within a unit space, effectively reflecting the shrinkage or expansion trend of root strips under stress. During the calculation process, if the root strip coverage in a certain area increases from 40% to 70%, the corresponding variation amplitude for that grid cell is 30%. The variation amplitudes of all cells are then spatially weighted and integrated (specific settings are made under special planting conditions, such as fertilization; by default, the weights are equal) to form root strip variation amplitude indicators for the corresponding sampling points, arranged sequentially by time labels, ultimately outputting a time-series parameter table of root strip changes covering all sampling points and monitoring periods. The outward expansion rate of the contour and the variation amplitude of root density are organized into a soil degradation trend identification parameter set according to the sampling point location and time label.

[0077] In step S6, the process of calculating the degradation score of red soil paddy fields and generating a soil risk level map based on the index sequence output by the degradation trend identification parameter set specifically includes:

[0078] After constructing the parameter set for identifying degradation trends, the data of each included indicator were first standardized. This parameter set consists of two main indicators: the outward expansion rate of the overlapping area boundary and the variation in the number of sampling points for root zone distribution within the overlapping area. Each indicator has different numerical scales and distribution characteristics at different sampling points and time periods, therefore, it needs to be converted to a unified dimension. Specifically, the zero-mean unit variance method was used, with the mean and standard deviation of each indicator across all sampling points and time periods as a benchmark to calculate standardized values, ensuring that each indicator has an equal starting point for contribution in subsequent weighted scoring. Subsequently, a multi-indicator weighted scoring mechanism was introduced, pre-setting the relative weight of each indicator to the soil degradation process. For example, in the high-temperature, high-humidity, and deep-reduction-prone red soil paddy field environment, the outward expansion rate of the overlapping area outline is considered to have a significant impact on the degree of root zone disturbance, therefore it is assigned a weight of 0.6; while the variation in root zone density at sampling points reflects the degree of shrinkage in the actual usable space of the root zone, and is assigned a weight of 0.4. The sum of all weights is 1. Subsequently, the standardized index values ​​are multiplied by their corresponding weights and summed to obtain the soil degradation score for each sampling point in each time period, forming a continuous degradation score sequence, which constitutes the basic data for risk level mapping.

[0079] After obtaining the degradation score sequences of all sampling points at various time periods, these need to be converted into practically readable risk levels to assist agricultural managers in intuitive identification and differentiated management. To this end, a mapping range for soil degradation risk levels is first established. This mapping range is formulated based on long-term regional monitoring data and degradation threshold research results, dividing the degradation score into five risk level ranges: 0.0–0.2 corresponds to “extremely low risk,” 0.2–0.4 corresponds to “low risk,” 0.4–0.6 corresponds to “medium risk,” 0.6–0.8 corresponds to “high risk,” and 0.8–1.0 corresponds to “extremely high risk.” The grading standard needs to meet interpretability and field guidance significance. The degradation score value of each sampling point in each time period is mapped to the above range, outputting the corresponding risk level label. Since the degradation score has time series characteristics, the risk level of each sampling point throughout the entire assessment period is recorded to ensure that the map reflects the most severe degradation situation. Finally, by combining the horizontal and vertical spatial coordinates of the sampling points in the red soil paddy field, the points with risk level labels are mapped to the actual spatial distribution area of ​​the paddy field, forming a visualization layer.

[0080] The above formulas are all dimensionless calculations. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters and thresholds in the formulas are set by those skilled in the art according to the actual situation.

[0081] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented using software, the above embodiments can be implemented, in whole or in part, as a computer program product. The computer program product includes one or more computer instructions or computer programs. When the computer instructions or computer programs are loaded or executed on a computer, all or part of the processes or functions described in the embodiments of this application are generated. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. The computer instructions can be stored in a computer-readable storage medium or transmitted from one computer-readable storage medium to another. For example, the computer instructions can be transmitted from one website, computer, server, or data center to another website, computer, server, or data center via wired (e.g., infrared, wireless, microwave, etc.) means. The computer-readable storage medium can be any available medium that a computer can access or a data storage device such as a server or data center that includes one or more sets of available media. The available medium can be a magnetic medium (e.g., floppy disk, hard disk, magnetic tape), an optical medium (e.g., DVD), or a semiconductor medium. A semiconductor medium can be a solid-state drive.

[0082] Those skilled in the art will recognize that the modules and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0083] Those skilled in the art will understand that, for the sake of convenience and brevity, the specific working processes of the systems, devices, and modules described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0084] In the several embodiments provided in this application, it should be understood that the disclosed systems, apparatuses, and methods can be implemented in other ways. For example, the apparatus embodiments described above are merely illustrative; for instance, the division of modules is only a logical functional division, and in actual implementation, there may be other division methods. For example, multiple modules or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the coupling or direct coupling or communication connection shown or discussed may be through some interfaces; the indirect coupling or communication connection between apparatuses or modules may be electrical, mechanical, or other forms.

[0085] The modules described as separate components may or may not be physically separate. The components shown as modules may or may not be physical modules; they may be located in one place or distributed across multiple network modules. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs.

[0086] In addition, the functional modules in the various embodiments of this application can be integrated into one processing module, or each module can exist physically separately, or two or more modules can be integrated into one module.

[0087] If the aforementioned functions are implemented as software functional modules and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, or a portion of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of this application. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0088] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

[0089] In conclusion, the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A soil risk assessment method for red soil paddy fields, characterized in that, Includes the following steps: S1. Collect the heat flux and oxygen concentration of the soil at pre-set sampling points in the red soil paddy field below the topsoil layer; S2. Calculate the unit volume metabolic heat release power curve at different depths of the sampling point, and identify the reduction-dominant region by combining the time change trend of oxygen concentration. S3. Spatial fitting is performed on the boundary of the reduction-dominant region to construct a three-dimensional distribution map of the deep hypoxia-dominant region. S4. Construct a spatial distribution map of crop roots in red soil paddy fields, spatially register it with the three-dimensional distribution map of the hypoxia-dominant area, and extract the boundary contours of the overlapping areas. S5. Calculate the boundary contour expansion rate and root density change amplitude of the overlapping area to construct a set of parameters for soil degradation trend identification. S6. Based on the index sequence output by the degradation trend identification parameter set, calculate the degradation score of red soil paddy field soil and generate a soil risk level map. In step S2, the metabolic heat release power curve per unit volume at different depths of the sampling point is calculated, and combined with the time-varying trend of oxygen concentration, the reduction-dominant region is identified, specifically including: Vertical heat flux is derived from temperature sequences at different depths at the same time point by spatial difference, and a heat flux distribution array is constructed with time as the horizontal axis and depth as the vertical axis. Perform central differential calculations on the constructed heat flux distribution array in chronological order to obtain the rate of change of heat flux in sub-segments between a set number of consecutive adjacent depth layers; The rate of change is the vertical heat flux gradient of each sub-section. The rate of change was converted into the metabolic heat release power corresponding to the sub-segment volume, and the time evolution curve of the metabolic heat release power of the sub-segment was established. In the parameter set, the time period corresponding to the heat flux change of each sub-segment is obtained, and the concentration change rate is extracted from the oxygen concentration sequence within the same time period to generate the oxygen consumption trend curve of the sub-segment. Align the two types of curves with a uniform time resolution and calculate the covariance coefficient to identify sub-segments with increased metabolic heat release power and enhanced oxygen consumption, and mark them as reduction-dominant regions. The covariance coefficients are numerical expressions after covariance normalization.

2. The soil risk assessment method for red soil paddy fields according to claim 1, characterized in that, In step S1, the collection of heat flux and oxygen concentration of soil at pre-set sampling points in red soil paddy fields below the topsoil layer specifically includes: Multiple evenly distributed sampling points were pre-set in red soil paddy fields to establish a planar coordinate system; Below the topsoil layer, multiple depth layers were set in the longitudinal profile of the sampling points in the red soil paddy field. Temperature and oxygen concentration were collected at a fixed sampling frequency, and data from different depth layers were recorded. All depth layer data are arranged according to acquisition time to form a parameter set with a unified time base and depth label.

3. The soil risk assessment method for red soil paddy fields according to claim 1, characterized in that, The derivation of vertical heat flux from temperature sequences at different depths at the same time point using spatial difference specifically includes: Perform differential calculations on the temperature sample values ​​of any two consecutive depth points in a set profile to obtain a vertical temperature gradient sequence; Obtain the thermal diffusivity and volumetric heat capacity parameters corresponding to the red soil type, and convert the temperature gradient into heat flux value through the heat transfer calculation formula; All heat flux values ​​are aggregated to construct a heat flux profile. The vertical structure of the heat flux is extracted at each time step, and a standardized heat flux distribution array is constructed.

4. The soil risk assessment method for red soil paddy fields according to claim 1, characterized in that, When identifying sub-segments with increased metabolic heat release power and enhanced oxygen consumption, if multiple sub-segments meet the criteria, the sub-segment with the largest covariance coefficient is selected as the reduction-dominant region.

5. The soil risk assessment method for red soil paddy fields according to claim 1, characterized in that, In step S3, spatial fitting of the boundary of the reduction-dominant region to construct a three-dimensional distribution map of the deep hypoxia-dominant region specifically includes: Based on the identified dominant restoration region in each sampling point, the upper boundary depth value and lower boundary depth value of the dominant restoration region are extracted and used as the upper boundary point set and lower boundary point set, respectively. All boundary points are aligned according to the planar coordinate system of their sampling points, and spatial fitting is performed on each point to obtain the upper and lower boundary spatial fitting surfaces covering the dominant region of the entire sampling point restoration. The stable spatial envelope region between the two fitted surfaces is used as a three-dimensional distribution map dominated by hypoxia.

6. The soil risk assessment method for red soil paddy fields according to claim 1, characterized in that, In step S4, a spatial distribution map of crop roots in red soil paddy fields is constructed, and this map is spatially registered with a three-dimensional distribution map of the hypoxia-dominant area. Extracting the boundary contours of overlapping areas specifically includes: The root depth of crops at each sampling point is collected, and the root zone distribution at different depths is constructed with the root location of the crop at each sampling point as a reference. The three-dimensional distribution map of the hypoxia-dominant area was retrieved, and the coordinates were aligned with the root zone distribution at different depths using an affine transformation algorithm. The spatial intersection of the overlapping areas was calculated and the contours were extracted.

7. A soil risk assessment method for red soil paddy fields according to claim 1, characterized in that, In step S5, calculating the boundary contour expansion rate and root density change magnitude of the overlapping area to construct a soil degradation trend identification parameter set specifically includes: Within a preset continuous time period, differential analysis is performed on the boundary contour of the overlapping area in the time dimension to obtain the contour expansion rate. Simultaneously analyze the variation of root density over time in the root zone distribution of all sampling points in the overlapping region; The two types of indicators, namely the rate of outward expansion of the contour and the magnitude of change in root density, were organized into a set of parameters for identifying soil degradation trends according to the location of the sampling point and the time period label.

8. A soil risk assessment method for red soil paddy fields according to claim 1, characterized in that, In step S6, the process of calculating the degradation score of red soil paddy fields and generating a soil risk level map based on the index sequence output by the degradation trend identification parameter set specifically includes: The time series of each indicator in the degradation trend identification parameter set is standardized, and a multi-indicator weighted score is performed based on the preset degradation judgment weight. A preset soil degradation risk level mapping range is used to map the scoring results to different soil degradation risk levels. In the overlapping area, the sampling points corresponding to all risk levels are mapped to the actual spatial location of the red soil paddy field, generating a soil degradation risk level map covering the entire red soil paddy field area.

Images