Citrus plantation soil fertility characterization method and system

Abstract

The invention provides a characterization method and system for soil fertility of a citrus plantation, and relates to the technical field of characterization of soil fertility, and the characterization method comprises the following steps: collecting deep soil samples according to uniform grids in the citrus plantation, arranging three-parameter sensors, and obtaining a multispectral remote sensing image at the same time; secondly, respectively constructing a space continuous function, a time continuous function and a spectrum continuous function by using common Kriging interpolation, Gaussian process regression and spectral index inversion; then selecting high-yield area samples with the first 10% of yield in continuous three years, generating high fertility reference distribution based on kernel density estimation, and fusing the three types of functions by using an entropy evaluation method to obtain a comprehensive fertility function; finally, principal components are extracted through function principal component analysis, the Tukey depth of the principal components relative to reference distribution is calculated, a soil fertility index is generated, and five-level fertility areas are divided through an improved natural fracture method.

Description

Technical Field

[0001] This invention relates to the field of soil fertility characterization technology, specifically to a method and system for characterizing soil fertility in citrus orchards. Background Technology

[0002] Citrus, as one of the main economic fruit trees cultivated in southern my country, is heavily influenced by soil fertility in terms of yield and quality. Due to long-term improper crop rotation, extensive farming practices, and uneven fertilization, the spatial distribution of soil nutrients within citrus orchards exhibits significant heterogeneity. This heterogeneity directly leads to uneven growth of citrus trees, unstable fruiting capacity, and declining fruit quality, severely hindering the improvement of the citrus industry's quality and efficiency. Therefore, accurately assessing the spatial pattern of soil fertility within citrus orchards and constructing an indicator system reflecting the fertility levels of different regions is fundamental to achieving precision fertilization and integrated water and fertilizer management. Especially against the backdrop of the accelerated advancement of smart agriculture and agricultural informatization, relying on multi-source data, such as soil sensors, remote sensing imagery, and laboratory analysis data, to conduct spatial continuous modeling and dynamic monitoring of soil fertility has become a core element in realizing differentiated orchard management and precision agricultural cloud decision-making.

[0003] Current research on citrus orchard soil fertility primarily employs laboratory physicochemical analysis methods based on field sampling, combined with geostatistics, such as ordinary kriging interpolation, to spatially complete soil properties. For continuous temporal monitoring, three-parameter soil sensors are often used, combined with Gaussian process regression to predict future conditions. In remote sensing, normalized difference vegetation index (NDVI) or soil diagnostic index is frequently used as proxy indicators to indirectly infer soil fertility levels by reflecting vegetation cover. However, these methods are typically conducted independently, lacking an effective data fusion mechanism, particularly in the fusion of multi-source information, where a unified function expression and dynamic weight adjustment mechanism are lacking. Furthermore, the evaluation process often relies on subjective scoring or fixed empirical values ​​to assign weights to each dimension, making it difficult to address the challenge of varying information quality in real-world scenarios.

[0004] The existing technology has the following key shortcomings: First, the fertility expression method based on a single information dimension cannot comprehensively reflect the true fertility status of the target area, resulting in large spatial prediction errors and insufficient basis for management zoning; Second, there is a lack of a method to dynamically determine the contribution of each dimension based on the information content of the data itself, which makes the final fertility function construction process subject to human bias and affects the objectivity of the evaluation; Third, existing studies rarely establish quantifiable reference distribution models from the perspective of high-yield reference samples, which makes it impossible to accurately assess the similarity between any region and typical high-yield soils, and also makes it impossible to scientifically guide the classification of fertility levels and the optimization of regional management.

[0005] The information disclosed in the background section is only intended to enhance the understanding of the background of this disclosure, and therefore may include information that does not constitute prior art known to those skilled in the art. Summary of the Invention

[0006] The purpose of this invention is to provide a method and system for characterizing soil fertility in citrus orchards, in order to solve the problems mentioned in the background art.

[0007] To achieve the above objectives, the present invention provides the following technical solution: A method for characterizing soil fertility in citrus orchards, comprising the following steps: Step 1: Divide the target citrus orchard into a uniform sampling grid, set up a sampling point at the center of each grid to collect soil samples at a depth of W, and conduct laboratory tests on the soil samples to obtain the soil physicochemical index values ​​at each point. Simultaneously, deploy a three-parameter soil sensor at the center of each grid to monitor soil moisture, electrical conductivity and temperature data, and acquire multispectral remote sensing images covering the entire target citrus orchard. Step 2: Use ordinary kriging interpolation to construct spatially continuous functions for soil physicochemical indices at each point, apply Gaussian process regression to sensor data to generate temporally continuous functions, perform atmospheric correction and geometric registration on remote sensing images to calculate normalized vegetation index and soil diagnostic index, and then construct spectral continuous functions. Step 3: Select soil samples from the top 10% of citrus production in each of the past three years as reference samples for high-yield areas. Based on these reference samples, construct a high fertility reference distribution with kernel density estimation. Extract the spatial function values, temporal function values, and spectral function values ​​of the soil samples and combine them into a comprehensive fertility function using the entropy method. Step 4: Perform principal component decomposition on the comprehensive fertility function, extract the first three principal components to obtain the dimension-reduced vector, and calculate the depth value of the dimension-reduced vector relative to the high fertility reference distribution using the Tukey depth function to generate the soil fertility index. Use the Jenks natural fracture classification method to divide each grid into five levels according to the soil fertility index, where level 4 and level 5 represent high fertility areas and level 1 represents low fertility areas.

[0008] Furthermore, the logic for dividing the target citrus orchard into a uniform sampling grid and setting a sampling point at the center of each grid to collect soil samples at a depth of W is as follows: Within the target citrus orchard, a uniform square sampling grid was established using a dynamic grid partitioning algorithm: the grid side length was defined. Total area of ​​the target citrus plantation The functional relationship is as follows: ,in, Indicates rounding up; Set origin At the northwest corner of the target citrus orchard, the center coordinates of each grid are set as follows: ; in, Let i be the coordinates of the center of the i-th grid. Here, represents the row and column number of the i-th grid, respectively, where i is the grid index. , This represents the total number of sampling grids. At each center point At this location, samples were collected at a depth of [depth value missing] using vertical drilling. Furthermore, soil column samples covering the main root system of citrus were collected, and the following five soil physicochemical indicators were measured in the laboratory: soil pH, organic matter content, total nitrogen content, available phosphorus content, and available potassium content. Three-parameter wireless sensor nodes were deployed simultaneously at each grid center to monitor soil moisture, soil conductivity, and soil temperature. A drone equipped with a 5-band multispectral camera was used to retrieve vegetation and soil indices from the multispectral images.

[0009] Furthermore, the logic for performing ordinary kriging interpolation on the soil physicochemical index data is as follows: Based on soil physicochemical index data from the grid center, a spatial autocorrelation model is constructed to generate spatially continuous functions of the physicochemical indices. Leave-one-out cross-validation was used to evaluate the interpolation results; The humidity, conductivity, and temperature data monitored by the sensor are normalized, a time covariance function is constructed, and a time continuous function is generated. , For continuous time variables; Reflectance correction and geometric registration were performed on UAV multispectral imagery, and vegetation indices and soil diagnostic indices were fused to generate a spectral continuity function. ; The three types of functions mentioned above are unified to the same coordinate system and resolution, and data integration is achieved through spatiotemporal alignment and outlier removal.

[0010] Furthermore, the spatial function values, temporal function values, and spectral function values ​​of each grid are extracted and combined into a comprehensive fertility function using the entropy method. The target period for fruit enlargement in the citrus orchard is set as follows: Extract spatially continuous function values ​​for each grid. The values ​​of the time-continuous function and the spectral continuous function are in Integral average and Furthermore, it is dynamically weighted according to the entropy method and integrated into a comprehensive fertility function. : ; in, The weights are dynamically fused to reflect the differences in information content across spatial, temporal, and spectral dimensions. make This represents the output of the i-th grid in year b. For year b, calculate the 90th percentile of the output of all grids in that year. , The high-yield area sample set is defined as areas where the yield in each of the three consecutive years is no less than the 90th percentile of that year. Grid: And satisfy the high-yield area sample set Total number of samples , Grid index for high-yield areas; For each high-yield grid The dimensionality-reduced vector of the high-yield area is obtained by using the PCA function. , The three-dimensional Epanechnikov kernel was selected, and the bandwidth was chosen according to the Silverman rule of thumb. Then refer to the distribution density for: ; in, Given an input vector with an arbitrary grid, Corresponding to 3D volume scaling, For indicator functions, It is the Euclidean norm.

[0011] Furthermore, the logic for performing principal component decomposition on the comprehensive fertility function and extracting the first three principal components to obtain the dimension-reduced vector is as follows: For each grid center During the full growing season of the target citrus orchard The observed comprehensive fertility function Notation: ; Perform principal component analysis on the function and extract the first three principal components to obtain the dimensionality-reduced vector. Define a dimension reduction vector Based on high fertility reference distribution Tukey depth; Measured by Tukey depth Similarity to the sample set of high-fertility, high-yield areas The formula used is: ; in, For inclusion The set of all closed half-spaces, that is, the half-spaces in the principal component space divided by the hyperplane. The reference is distributed in a closed half-space. The probability mass in; Introducing Mahalanobis distance as a fast approximation: ; in, , ; These are the mean and covariance of the reference distribution in high-yield areas, representing the spatial distribution center and morphology of the principal components; Mahaldo depth Mapped to soil fertility index : ; in, When the fertility at that location is most similar to the reference distribution of high-yield areas, the fertility is optimal; when... When the fertility at that location is significantly lower than the benchmark for high-yield areas, it indicates that the fertility at that location is significantly lower than that of high-yield areas.

[0012] Furthermore, the logic for classifying it into five levels using Jenks' natural fracture classification method is as follows: SFI values ​​for all grids Sort the SFI values ​​from smallest to largest; divide them into 5 categories, with the category boundaries as follows: Let o be the category index, and let the o-th category be: ; in, This refers to the index of the grid's position in the sorted sequence; The improved Jenks natural fracture method is used to classify fertility levels, and the objective function is defined. The weighted sum of intra-class tightness and inter-class separation: ; in, It is the mean of class o. The weighting is determined by the proportion of samples in class o to the total number of samples, making the intra-class tightness more focused on the intermediate stable samples. Let be the number of samples in class o. It serves as a balancing factor for inter-class separation. The simulated annealing algorithm is used to search for the optimal classification boundary, and the boundary is finally obtained. Then, the optimal five-level classification result was obtained. arrive ,satisfy arrive It is a highly fertile area and spatially continuous; The soil fertility classification rules for the target planting area are as follows: grade SFI range is This is a priority improvement zone; [Level] SFI range is This area requires additional fertilization; grade SFI range is This is a routine management area; level SFI range is To optimize the maintenance zone; level SFI range is It is a high-yield protected area.

[0013] The present invention also provides a soil fertility characterization system for citrus orchards, the system being used to perform the above-described method for characterizing soil fertility in citrus orchards, comprising: The data acquisition module is used to divide the target citrus orchard into a uniform sampling grid, set a sampling point at the center of each grid to collect soil samples at a depth of W, and conduct laboratory tests on the soil samples to obtain the soil physicochemical index values ​​at each point. Simultaneously, a three-parameter soil sensor is deployed at the center of each grid to monitor soil moisture, electrical conductivity and temperature data, and acquire multispectral remote sensing images covering the entire target citrus orchard. The data processing module is used to construct spatially continuous functions of soil physicochemical indicators at each point using ordinary Kriging interpolation, generate time-continuous functions by applying Gaussian process regression to sensor data, calculate normalized vegetation index and soil diagnostic index after atmospheric correction and geometric registration of remote sensing images, and then construct spectral continuous functions. The function calculation module is used to select soil samples from the top 10% of citrus production in each of the past three years as reference samples for high-yield areas. Based on these reference samples, a high fertility reference distribution with kernel density estimation is constructed. Spatial function values, temporal function values, and spectral function values ​​of the soil samples are extracted and combined into a comprehensive fertility function using the entropy method. The fertility assessment module is used to perform principal component decomposition on the comprehensive fertility function, extract the first three principal components to obtain a dimensionality-reduced vector, and calculate the depth value of the dimensionality-reduced vector relative to the high fertility reference distribution using the Tukey depth function to generate a soil fertility index. The Jenks natural fracture classification method is used to divide each grid into five levels according to the soil fertility index, where level 4 and level 5 represent high fertility areas and level 1 represents low fertility areas.

[0014] Compared with the prior art, the beneficial effects of the present invention are: This invention achieves the acquisition and correlation of spatial, temporal, and spectral data in the same coordinate system by simultaneously deploying a three-parameter soil sensor at the center of a sampling grid and acquiring laboratory soil samples and multispectral remote sensing data. This ensures the consistency of the basic data for subsequent fusion analysis. Three-dimensional continuous functions are constructed using ordinary Kriging interpolation, Gaussian process regression, and remote sensing spectral index inversion, respectively. This effectively solves the problems of spatial sparsity of soil physicochemical indicators, temporal discreteness of sensor data, and inconsistency in spatial accuracy of remote sensing images, laying a functional foundation for comprehensive evaluation.

[0015] This invention constructs a kernel density estimation reference distribution using high-yield grids that rank in the top 10% for three consecutive years, representing ideal fertility characteristics. It also constructs a multidimensional comprehensive fertility function using the entropy method, which automatically adjusts the weights of each dimension according to actual differences, effectively avoiding bias from human weighting and improving the scientific nature of the fusion results. This invention extracts multidimensional feature vectors for each region using principal component analysis and Tukey depth calculation, and generates a soil fertility index based on the depth value of the region in the high-yield reference distribution. This achieves a classification method based on similarity rather than absolute value. Combined with Jenks' natural fracture classification method, the region is divided into five levels. This improves the clarity of classification boundaries while achieving spatial aggregation, providing a guiding basis for subsequent fertilization optimization and agronomic management. Attached Figure Description

[0016] Figure 1 This is a schematic diagram of the overall method flow of the present invention; Figure 2 This is the fitting curve of the soil fertility index-comprehensive fertility function of the present invention; Figure 3 This is a histogram of the comprehensive fertility function of the present invention; Figure 4 This is the reference distribution density-relative distance-influence intensity curve for this invention; Figure 5 This is the fitting curve of the comprehensive fertility function-score coefficient of the present invention; Figure 6 This is the soil fertility index-score coefficient fitting curve of the present invention; Figure 7 This is a flowchart of the overall system modules of the present invention. Detailed Implementation

[0017] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments.

[0018] It should be noted that, unless otherwise defined, the technical or scientific terms used in this invention should have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in this invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed following the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0019] Example: Please see Figures 1-6 The present invention provides a technical solution: A method for characterizing soil fertility in citrus orchards, comprising the following steps: Step 1: Divide the target citrus orchard into a uniform sampling grid, set up a sampling point at the center of each grid to collect soil samples at a depth of W, and conduct laboratory tests on the soil samples to obtain the soil physicochemical index values ​​at each point. Simultaneously, deploy a three-parameter soil sensor at the center of each grid to monitor soil moisture, electrical conductivity and temperature data, and acquire multispectral remote sensing images covering the entire target citrus orchard. The logic for dividing the target citrus orchard into a uniform sampling grid and setting a sampling point at the center of each grid to collect soil samples at a depth of W is as follows: Within the target citrus orchard, a uniform square sampling grid was established using a dynamic grid partitioning algorithm: the grid side length was defined. Total area of ​​the target citrus plantation The functional relationship is as follows: ,in, This indicates rounding up. Small-area park, adopting High-density mesh, for The large-scale park adopts Low-density grid; Based on the total area of ​​the target citrus orchard Automatic grid density adjustment: For small-area planting areas, i.e. More detailed measurements are needed; the minimum grid side length should be 10m. This is suitable for medium-sized planting areas. The grid is divided into segments based on area, with a grid side length of [missing information]. For large-scale planting areas, that is A coarser grid can be used, with the side length set to the maximum value of 30m; Grid side length Determines the center spacing of the grid. The larger the size, the wider the grid and the sparser the center of the grid, making it suitable for large-area scenarios with low sampling cost requirements; The smaller the size, the denser the mesh, with a denser mesh center, suitable for small-area scenarios requiring high precision; when When smaller, , ,when Larger makes hour, The calculation will Cut off to 30m, for the intermediate case, along with It increases in size in a step-like manner; Set origin At the northwest corner of the target citrus plantation, numbered by row and column. The center coordinates of each grid cell are calculated sequentially, and the center coordinates of each grid cell are set as follows: ; in, This represents the coordinates of the center of the i-th grid. Here, represents the row and column number of the i-th grid, respectively, where i is the grid index. , This represents the total number of sampling grids. The x-direction increases from west to east, and the y-direction decreases from north to south; each step moves one unit. Each unit ensures that its center point falls precisely within its respective grid cell. At each center point At this location, samples were collected at a depth of [depth value missing] using vertical drilling. Soil columnar samples covering 20–60 cm of the main root system of citrus trees were collected. The following five soil physicochemical indicators were measured in the laboratory: soil pH, organic matter content, total nitrogen content, available phosphorus content, and available potassium content. Soil pH was calculated using the potentiometry method, organic matter content was calculated using the Walkley-Black wet oxidation method, total nitrogen content was calculated using the Kjeldahl method, available phosphorus content was calculated using the Olsen method, and available potassium content was calculated using the flame photometry method. The central value of 0.6m corresponds to 60cm, which can cover the 20-60cm range of citrus roots. This indicates a 10cm allowance for flexibility to adapt to actual terrain conditions and depth. The larger the depth, the better the assessment of deeper soil properties, but this also increases cost and disturbance. The smaller the sample size, the more likely it is to miss the impact of the root system by only sampling the shallow layer. Soil pH reflects the soil's acid-base balance; organic matter content reflects the humus and fertility base in the soil, with a higher content indicating better soil fertility; total nitrogen content reflects the nitrogen nutrient supply capacity, with a higher content indicating a more abundant nitrogen supply; available phosphorus content reflects the proportion of phosphorus that can be absorbed by plants, with a higher content indicating more abundant phosphorus; and available potassium content reflects the potassium element that plants can utilize, with a higher content indicating a more abundant potassium fertilizer supply. Each dependent variable (index) is a measured value, while the independent variable is the physical and chemical state of the sample soil itself. Three-parameter wireless sensor nodes are deployed synchronously at each grid center to monitor soil moisture, soil conductivity and soil temperature. The monitoring frequency is set to once per hour, and the data is collected synchronously with the timestamp aligned to the hour. The higher the soil moisture content, electrical conductivity, and dissolved salt concentration in the monitoring area, the more likely it is to affect root activity and microbial metabolism. Using a drone equipped with a 5-band multispectral camera, flying at an altitude of 120 meters, with a ground resolution of less than or equal to 0.5 meters, based on the reflectance data of the ground calibration plate, eliminating the influence of atmospheric scattering, the true surface reflectance was calculated, and the normalized vegetation index and soil diagnostic index were calculated to reflect vegetation coverage and organic matter content, respectively. The higher the normalized vegetation index (NDI) of the target planting area, the stronger the chlorophyll blocks red light and reflects near-infrared light, indicating more lush vegetation. The lower the NDI, the less vegetation or water body there is. The higher the soil diagnostic index, the higher the proportion of bare soil reflectance, indicating that the organic matter may be lower. The lower the soil diagnostic index, the richer the organic matter and the better the soil cover.

[0020] Step 2: Use ordinary kriging interpolation to construct spatially continuous functions for soil physicochemical indices at each point, apply Gaussian process regression to sensor data to generate temporally continuous functions, perform atmospheric correction and geometric registration on remote sensing images to calculate normalized vegetation index and soil diagnostic index, and then construct spectral continuous functions. The logic for performing ordinary kriging interpolation on soil physicochemical index data is as follows: The steps for spatial interpolation analysis of soil physicochemical indicators measured in the laboratory include: calculating the experimental semivariance value at different spatial distances based on the physicochemical indicator data of each grid center, quantifying the spatial variation characteristics, selecting a spherical model to describe the spatial variation law, and determining the model parameters through nonlinear optimization, including nugget value, sill value and range, wherein the nugget value represents random variation, the sill value represents total spatial variation, and the range represents the spatial autocorrelation range; Calculate the experimental semivariance for each physicochemical index: ; in, This is the experimental semivariance value, used to measure the distance between grid centers. The degree of numerical difference, It is the spatial distance between the center pairs of the grid. For all distances The grid center pair set, The measured values ​​of the physicochemical properties at the center of the i-th grid are... Let be the measured value of the physicochemical properties of the j-th grid center point, and let be the distance between it and point i. ; The larger the value, the stronger the spatial variability of soil properties at that scale, and the weaker the correlation. hour, This indicates that there is no difference at the same location; The theoretical semivariance is fitted using a spherical model: ; in, It is the nugget value, representing measurement error or microscale variation. It is the sill value, representing the degree of spatial variability in the population. It is the range, representing the maximum distance at which spatial correlation disappears; when hour, The increase with increasing distance indicates enhanced mutation. When, it indicates that the maximum mutation has been reached; when At that time, space is no longer relevant; Based on the spatial autocorrelation model, a system of linear equations is constructed to solve for the weight coefficients at each location, generating a spatial distribution map of physicochemical indicators with a resolution of 1 meter, i.e., a spatially continuous function. The interpolation results were evaluated using leave-one-out cross-validation to ensure that the root mean square error did not exceed 1.5 times the laboratory measurement error. For any position The formula for ordinary Kriging interpolation prediction is: ; in, For predicted values, It is the weight coefficient of the i-th grid center, reflecting the strength of spatial correlation. It is the measured value at the center of the i-th grid. This represents the total number of grid centers; The weights are solved using the following linear system: ; in, Here, q represents the Lagrange multiplier used to keep the weight sum to 1, and q is the equation number in the system of equations. , The larger the value, the greater the influence of the i-th point on the prediction. Vector symbol; The leave-one-out method is used to verify the interpolation accuracy, and the root mean square error is calculated: ; in, For predicted values, These are experimental values; requirements for each indicator are as follows: Less than or equal to 1.5 times the laboratory measurement error; The steps for time-series modeling of soil moisture, conductivity and temperature data monitored by sensors include: performing mean removal and variance normalization on the raw data to eliminate dimensional differences, constructing a time covariance function, quantifying the correlation strength and noise impact at different time points, including signal variance, time scale parameter and noise variance, wherein the signal variance reflects the temporal fluctuation amplitude, the time scale parameter controls the correlation decay rate, and the noise variance characterizes the measurement error; Define a Gaussian process where the time covariance function uses a squared exponential kernel and a white noise term: ; in, The interval between two time points. For signal variance, Indicates the range of change in soil parameters. Using time as the time scale, the rate of correlation decay is controlled; a larger value indicates a more stable change. The noise variance represents the magnitude of the measurement error. For a unit impulse function, when The value is 1 for time and 0 for everything else; By maximizing the marginal likelihood function, the model parameters are automatically adjusted to ensure that the confidence interval coverage of the time prediction is greater than or equal to 95%, generating a continuous prediction curve for soil parameters with a time resolution of 10 minutes, i.e., a time-continuous function. This includes the mean function and the range of uncertainty. For continuous time variables; The optimal parameters are obtained by maximizing the marginal likelihood function. : ; in, Let be the vector of sensor measured values, and be the covariance matrix, which depends on the kernel function and the time difference. Given the set of hyperparameters to be optimized, the optimal model parameters can be obtained by maximizing this function. The normalization constant is denoted by log; in the formula, log represents the natural logarithm, i.e. The L-BFGS algorithm is used for iterative solution, with a convergence threshold set to... ; Regarding time-based forecasts and uncertainties: ; in, For the future The predicted value, Let covariance be the vector between the predicted point and the observed point. The variance of the predicted point itself; the mean and confidence interval are provided for uncertainty analysis; The steps for spectral analysis of multispectral images acquired by UAVs include: using reflectance data from a ground calibration plate to eliminate atmospheric scattering and absorption effects, calculating the true reflectance of the ground surface, and using a polynomial transformation model to correct geometric distortion of the image based on the coordinate information of ground control points, with the registration residual controlled within 0.5 pixels. The normalized vegetation index (NVI) is calculated by inverting vegetation cover using the near-infrared to red band reflectance ratio, and the soil diagnostic index is calculated by assessing soil organic matter content using the short-wave infrared to red band reflectance ratio. The NVI and soil diagnostic index are then fused to generate a continuous spectral function characterizing the spectral features. ; Atmospheric correction uses the 6S radiative transfer model to eliminate the effects of atmospheric scattering and absorption. Geometric registration is based on ground control points, and a quadratic polynomial transformation model is established. Calculation of Normalized Difference Vegetation Index The formula used is: ; in, The surface reflectance is in the near-infrared band. The surface reflectance is in the red light band. The larger the area, the higher the vegetation coverage; Calculate soil diagnostic index The formula used is: ; in, For shortwave infrared reflectivity, The larger the value, the higher the proportion of bare soil and the lower the organic matter content. The fusion spectrum continuity function is: ; in, The weighting parameters are dynamically optimized using the entropy method, adjusting the contribution ratio of vegetation and soil according to the growth stage. All data were converted to the UTM Zone 50N coordinate system with a grid resolution of 1 m × 1 m. Based on the sensor data timestamps, the Kriging interpolation results and spectral data were linearly interpolated and aligned. Data points that deviated from the spatial mean by more than 3 standard deviations were removed, as were data points that exceeded the Gaussian process prediction confidence interval by more than 2.5 standard deviations.

[0021] Step 3: Select soil samples from the top 10% of citrus production in each of the past three years as reference samples for high-yield areas. Based on these reference samples, construct a high fertility reference distribution with kernel density estimation. Extract the spatial function values, temporal function values, and spectral function values ​​of the soil samples and combine them into a comprehensive fertility function using the entropy method. The logic for extracting the spatial function values, temporal function values, and spectral function values ​​of each grid and combining them into a comprehensive fertility function using the entropy method is as follows: The target period for fruit enlargement in the citrus orchard is set as follows: Extract spatially continuous function values ​​for each grid. The values ​​of the time-continuous function and the spectral continuous function are in Integral average and Furthermore, it is dynamically weighted according to the entropy method and integrated into a comprehensive fertility function. : ; in, To dynamically fuse weights, reflecting the differences in information content across spatial, temporal, and spectral dimensions, and ; ; This indicates the spatial distribution values ​​of soil physicochemical indicators at this location. The larger the value, the stronger the fertility signal in that dimension, and the better it is for... The greater the positive contribution.

[0022] Table 1 shows some sample data and specific data of the comprehensive fertility function.

[0023]

[0024] From sample 1 to sample 15, although the spatial continuity function values ​​fluctuated (e.g., rising from 1.062 to 1.466, then falling to 0.879 before rising again), overall, higher spatial continuity function values ​​tended to correspond to higher comprehensive fertility functions. For example, in sample 2... hour, It is also at a relatively high level; sample 15 hour, The value was significantly higher than most other points, indicating that spatial fertility characteristics made a significant contribution to overall fertility. Although time function The contribution weight to overall fertility and Same, but overall change The effects are relatively mild; for example, in sample 4... Higher, while the sample's It significantly outperformed neighboring samples (such as sample 3). , Partly due to the higher spectral values, it can be seen that... The fluctuation energy, under conditions similar to spatial and spectral characteristics, is... It serves as a fine-tuning tool; When the spectral function exhibits extreme values ​​in some samples, it can significantly influence overall fertility. For example, in sample 3... Significantly higher than surrounding samples, The fact that the spectral index exceeds that of samples 2 and 4 indicates that when the spectral index is high (possibly corresponding to better vegetation cover or soil organic matter), it will increase non-linearly. ; In the section from sample 9 to sample 11, we can see that all three indicators in sample 9 have declined. ),lead to It decreased to 0.999; while any indicator in samples 10 and 11 increased significantly (e.g., sample 11's...). ), making The rapid rebound to 1.375 further confirms the synergistic effect of three-dimensional indicators on comprehensive fertility. Although each dimension is related to The contribution is fixed, but due to the slight differences in the adaptive weights of the entropy method, the above samples... The significant difference between the maximum (1.375) and minimum (0.992) values ​​reflects that, under different environments, the fusion of spatial-temporal-spectral three-dimensional information can generate more distinctive fertility levels, providing strong support for refined fertilization and management.

[0025] Construct the original index matrix as follows ,in , , The exponent is normalized to: Entropy calculation: ; Weight confirmation: , , thus calculating , , ; in, It is the value of the u-th sample under the v-th index. It is the normalized weight, which measures the proportion of sample u in the v-type indicator. It is the information entropy of the v-th index. This refers to the information redundancy of the v-th type of indicator; the larger the value, the more significant the differences between the indicators. It is the entropy weight of the v-th type index; make This represents the output of the i-th grid in year b. For year b, calculate the 90th percentile of the output of all grids in that year. , The high-yield area sample set is defined as areas where the yield in each of the three consecutive years is no less than the 90th percentile of that year. Grid: And satisfy the high-yield area sample set Total number of samples , Grid index for high-yield areas; For each high-yield grid The dimensionality-reduced vector of the high-yield area is obtained by using the PCA function. , The three-dimensional Epanechnikov kernel was selected, and the bandwidth was chosen according to the Silverman rule of thumb. Then refer to the distribution density for: ; in, The bandwidth for kernel density estimation is optimized according to the Silverman criterion. Corresponding to 3D volume scaling, It is the three-dimensional principal component score of the nth high-yield area sample. For indicator functions; The input vector is an arbitrary grid. Points to be evaluated The kernel density estimate on the distribution of high-yield samples is larger, indicating that the point is more similar to the high-yield samples. The larger the value, the closer the fertility characteristics of the point are to the core characteristics of the high-yield area. When the value is close to 0, the fertility characteristics of the point deviate from the typical pattern of the high-yield area. middle, The score of the first principal component represents the direction of major fertility variation. The score of the second principal component represents the direction of minor fertility variation. The score of the third principal component represents the variation in residual fertility; middle, The projection value of the first principal component. The projection value of the second principal component. The projection value of the third principal component; The core terms of the kernel function are: ,in ,when ,Right now , The maximum impact intensity indicates that the point to be evaluated completely overlaps with the sample from the high-yield area; when , The effect is a 25% decrease in intensity; when , The effect disappears, indicator function It is truncated to 0; when the spatial distribution of samples in high-yield areas is dense, that is When the bandwidth is larger, the Silverman criterion automatically reduces the bandwidth. Decreasing the kernel size shrinks the effective range of the kernel function, making density estimation more focused on local features; conversely, increasing the kernel size increases the bandwidth and covers a wider area.

[0026] The specific data of some sample data and reference distribution density are shown in Table 2.

[0027]

[0028] The table shows relative distances. Epanechnikov nuclear influence intensity The corresponding three-dimensional principal component vectors and the simplified reference density As the sample numbers increase from 1 to 15, the parameters... As the value gradually increases from 0.000 to 0.483, the Epanechnikov kernel function value gradually decreases from a maximum of 0.7500 to 0.6341. This is due to the effect of the squared term. It is a parabola opening downwards, with a relatively slow initial decay (e.g., from...). Up to 0.483, (From 0.7491 to 0.7357) and then the decay accelerated (e.g. from) Up to 0.483, (From 0.6881 to 0.6341).

[0029] The reference distribution density is simplified to Therefore, with Increase It also shows the same non-linear decreasing trend, from 0.7500 to 0.6341, indicating that the farther the point to be evaluated is from the center of the high-yield area, the lower its probability of belonging to the high-yield area or its similarity. The vector of each sample is a distance... The scaling result of unit vectors of equal length but random orientation, for example, sample 5. Its vector is (0.118, 0.019, 0.069), and its Euclidean norm is approximately 0.138; Sample 15 The norm of the vector (0.373, 0.280, -0.239) is also equal to 0.483; Overall, and The curve shape is a smooth decay, which conforms to the "strong center, weak edge" characteristic of kernel function design, and can be used to describe the similarity decay process between any point and the high-yield reference area.

[0030] Step 4: Perform principal component decomposition on the comprehensive fertility function, extract the first 3 principal components to obtain the dimension reduction vector, and use the Tukey depth function to calculate the depth value of the dimension reduction vector relative to the high fertility reference distribution to generate the soil fertility index. Use the Jenks natural fracture classification method to divide each grid into five levels according to the soil fertility index, where level 4 and level 5 represent high fertility areas and level 1 represents low fertility areas. The logic for performing principal component decomposition on the comprehensive fertility function and extracting the first three principal components to obtain the dimensionality-reduced vector is as follows: For each grid center During the full growing season of the target citrus orchard The observed comprehensive fertility function Notation: ; By using principal component analysis, a set of orthogonal characteristic functions is constructed. It satisfies the orthogonal normalization condition: ; in, It is the Kronecker delta function, and Orthogonality is The time integral is 0, and the normalization property is... The time integral is 1; For each grid i, calculate the score of the k-th principal component: ; in, The score is given on the k-th principal component in the i-th grid. Extract the first three principal components to obtain the dimensionality-reduced vector. Define a dimension reduction vector Based on high fertility reference distribution Tukey depth; This may reflect seasonal changes, such as increased fertility in spring. This may correspond to a rainfall response pattern. It may indicate the decline trend after fertilization; orthogonality ensures that different patterns are independent of each other; for The larger the value, the better the fertility change at that point matches the model; the smaller the value, the more it deviates from the typical change model. It is the actual fertility. It's a rainfall pattern; Indicates strong seasonality. Indicates out of season; Measured by Tukey depth Similarity to the sample set of high-fertility, high-yield areas The formula used is: ; in, For inclusion The set of all closed half-spaces, that is, the half-spaces in the principal component space divided by the hyperplane. The reference is distributed in a closed half-space. The probability mass in; It measures the minimum distance between a point and the core of a high-yield area. The closer the value is to 1, the more likely the point is located in the core of the distribution, i.e., the better the fertility. The closer the value is to 0, the more likely the point is located on the edge of the distribution, i.e., the worse the fertility. Introducing Mahalanobis distance as a fast approximation: ; in, , ; These are the mean and covariance of the reference distribution in high-yield areas, representing the spatial distribution center and morphology of the principal components; When the location is The Mahalanobis distance is 0.0. The optimal fertility state is achieved when the location is at one standard deviation, with a Mahalanobis distance of approximately 1.0. The fertility status is good; when the location is at 2 standard deviations, the Mahalanobis distance is approximately 2.0. The fertility status is moderate. When the location is at 3 standard deviations, the Mahalanobis distance is greater than 3.0. The soil fertility condition is one that needs improvement. Mahaldo depth Mapped to soil fertility index : ; in, The closer the value is to 1, the more similar the fertility at that location is to the reference distribution in high-yield areas, indicating optimal fertility; when... The closer the value is to 0, the more significantly the fertility at that location is lower than the benchmark for high-yield areas; The range and meaning are as follows: For optimal fertility , For excellent fertility , Medium fertility , Fertilization needs to be improved , Fertilization needs to be improved first. .

[0031] Table 3 shows some sample data and specific data on soil fertility index.

[0032]

[0033] As the score coefficient increased from -0.200 (Sample 3) to approximately 1.385 (Sample 2), the comprehensive fertility function also gradually increased from a minimum of 0.941 to a maximum of 1.389 (Sample 2), which validates the model's... The linear relationship; for example, the score coefficient of sample 15 is 1.077, corresponding to It shows a significant improvement compared to adjacent samples with lower score coefficients; Mahalanobis depth is strictly negatively correlated with the score coefficient; as the score coefficient increases, the denominator increases, leading to a decrease in Mahalanobis depth. Because In essence It is also negatively correlated with the score coefficient; a high score coefficient (meaning a greater deviation from the high-yield center) will result in a lower soil fertility index. pass The value reveals the logic behind the fertility level classification: when (e.g., Sample 3, Sample 11), designated as G5 (high-yield protection zone); when (e.g., samples 5, 10, and 12) are classified as G4; while when If the value is between 0.4 and 0.7 (e.g., for sample 1, sample 4, sample 6, etc.), then it is G3; lower values ​​indicate a lower value. It will correspond to G2 or G1; this classification reflects both the numerical change of Mahalanobis depth and the combined effect of the scoring coefficient and the center and deviation of the comprehensive fertility function, ensuring that the fertility zoning takes into account both linear mapping and nonlinear reduction mechanisms.

[0034] The logic behind classifying fractures into five levels using Jenks' natural fracture classification method is as follows: SFI values ​​for all grids Sort the SFI values ​​from smallest to largest; divide them into 5 categories, with the category boundaries as follows: Let o be the category index, and let the o-th category be: ; in, This refers to the index of the grid's position in the sorted sequence; The improved Jenks natural fracture method is used to classify fertility levels, and the objective function is defined. The weighted sum of intra-class tightness and inter-class separation: ; in, For intra-class compactness, For inter-class separation degree, It is the mean of class o. The weighting factor is the proportion of the number of samples in class o to the total number of samples. An index of 1.5 is used to increase the penalty for outliers, making the intra-class tightness more focused on the intermediate stable samples. Let be the number of samples in class o. It serves as a balancing factor for inter-class separation. It is a cost function that combines "intra-class dispersion" and "between-class separation," and its value reflects the quality of the current partitioning—the larger the value, the "worse" (intra-class dispersion, inter-class mixing). The larger the independent variables (sample bias, mean difference), the larger the corresponding penalty or reward, thus affecting... Size; When a certain sample point The mean of its category The difference increases. The item increases, leading to Increase; when the difference between the means of two classes increases. Increase The term increases, but this part is minimized. This will prompt the algorithm to choose a larger class interval; Intra-class tightness reflects the number of samples within each class. With the mean of this class The degree of dispersion between them, when a certain point The further away from its class mean, The larger the value, the greater the corresponding intra-class penalty; therefore, the larger the independent variable (the deviation between the sample value and the class mean), the larger this term will be. The inter-class separation reflects the distance between the means of two adjacent classes; the greater the difference between the means of two adjacent classes, the larger this term becomes. However, because it is added to... In fact, it is the separation reward that "needs to be subtracted". When minimizing it, it will prompt us to increase the mean gap. The larger the independent variable (mean gap), the larger this term is. The simulated annealing algorithm is used to search for the optimal classification boundary, with the following parameters set: initial temperature. Temperature decay coefficient Termination temperature In each iteration, at the current solution Above, a random perturbation of a certain boundary. ,Keep Calculate the objective function under the new solution. , and In comparison, if If the solution is correct, then accept the new solution; otherwise, use probability. Accept, to escape local optima; press Gradually cool down until ; The objective function value corresponding to the current classification boundary. The objective function value under the new perturbation classification boundary; Finally, the boundary is obtained. Then, the optimal five-level classification result was obtained. arrive ,satisfy arrive It is a highly fertile area and spatially continuous; The soil fertility classification rules for the target planting area are as follows: grade SFI range is The color code is red. This is a priority improvement zone; [Level] SFI range is Color code: orange This area requires additional fertilization; grade SFI range is The color code is yellow. This is a routine management area; level SFI range is The color code is light green. To optimize the maintenance zone; level SFI range is The color code is green. It is a high-yield protected area.

[0035] Please see Figure 7 The present invention also provides a soil fertility characterization system for citrus orchards, the system being used to perform the above-described method for characterizing soil fertility in citrus orchards, comprising: The data acquisition module is used to divide the target citrus orchard into a uniform sampling grid, set a sampling point at the center of each grid to collect soil samples at a depth of W, and conduct laboratory tests on the soil samples to obtain the soil physicochemical index values ​​at each point. Simultaneously, a three-parameter soil sensor is deployed at the center of each grid to monitor soil moisture, electrical conductivity and temperature data, and acquire multispectral remote sensing images covering the entire target citrus orchard. The data processing module is used to construct spatially continuous functions of soil physicochemical indicators at each point using ordinary Kriging interpolation, generate time-continuous functions by applying Gaussian process regression to sensor data, calculate normalized vegetation index and soil diagnostic index after atmospheric correction and geometric registration of remote sensing images, and then construct spectral continuous functions. The function calculation module is used to select soil samples from the top 10% of citrus production in each of the past three years as reference samples for high-yield areas. Based on these reference samples, a high fertility reference distribution with kernel density estimation is constructed. Spatial function values, temporal function values, and spectral function values ​​of the soil samples are extracted and combined into a comprehensive fertility function using the entropy method. The fertility assessment module is used to perform principal component decomposition on the comprehensive fertility function, extract the first three principal components to obtain a dimensionality-reduced vector, and calculate the depth value of the dimensionality-reduced vector relative to the high fertility reference distribution using the Tukey depth function to generate a soil fertility index. The Jenks natural fracture classification method is used to divide each grid into five levels according to the soil fertility index, where level 4 and level 5 represent high fertility areas and level 1 represents low fertility areas.

[0036] 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 in the formulas are set by those skilled in the art according to the actual situation.

[0037] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented in software, the above embodiments can be implemented, in whole or in part, as a computer program product. Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented by 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.

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

[0039] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes 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.

Claims

1. A method for characterizing soil fertility in citrus orchards, characterized in that, The specific steps include: Step 1: Divide the target citrus orchard into a uniform sampling grid. Set up a sampling point at the center of each grid to collect soil samples at a depth of W. Conduct laboratory tests on the soil samples to obtain the soil physicochemical index values ​​at each point. Simultaneously, deploy a three-parameter soil sensor at the center of each grid to monitor soil moisture, electrical conductivity and temperature data, and acquire multispectral remote sensing images covering the entire target citrus orchard. Step 2: Use ordinary kriging interpolation to construct spatially continuous functions for soil physicochemical indices at each point, apply Gaussian process regression to sensor data to generate temporally continuous functions, perform atmospheric correction and geometric registration on remote sensing images to calculate normalized vegetation index and soil diagnostic index, and then construct spectral continuous functions. Step 3: Select soil samples from the top 10% of citrus production in each of the past three years as reference samples for high-yield areas. Based on these reference samples, construct a high fertility reference distribution with kernel density estimation. Extract the spatial function values, temporal function values, and spectral function values ​​of the soil samples and combine them into a comprehensive fertility function using the entropy method. Step 4: Perform principal component decomposition on the comprehensive fertility function, extract the first three principal components to obtain the dimension-reduced vector, and calculate the depth value of the dimension-reduced vector relative to the high fertility reference distribution using the Tukey depth function to generate the soil fertility index. Use the Jenks natural fracture classification method to divide each grid into five levels according to the soil fertility index, where level 4 and level 5 represent high fertility areas and level 1 represents low fertility areas.

2. The method for characterizing soil fertility in a citrus orchard according to claim 1, characterized in that: The logic for dividing the target citrus orchard into a uniform sampling grid and setting a sampling point at the center of each grid to collect soil samples at a depth of W is as follows: Within the target citrus orchard, a uniform square sampling grid was established using a dynamic grid partitioning algorithm: the grid side length was defined. Total area of ​​the target citrus plantation The functional relationship is as follows: ,in, Indicates rounding up; Set origin At the northwest corner of the target citrus orchard, the center coordinates of each grid are set as follows: ； in, Let i be the coordinates of the center of the i-th grid. Here, represents the row and column number of the i-th grid, respectively, where i is the grid index. , This represents the total number of sampling grids. At each center point At this location, samples were collected at a depth of [depth value missing] using vertical drilling. Furthermore, soil column samples covering the main root system of citrus were collected, and the following five soil physicochemical indicators were measured in the laboratory: soil pH, organic matter content, total nitrogen content, available phosphorus content, and available potassium content. Three-parameter wireless sensor nodes were deployed simultaneously at each grid center to monitor soil moisture, soil conductivity, and soil temperature. A drone equipped with a 5-band multispectral camera was used to retrieve vegetation and soil indices from the multispectral images.

3. The method for characterizing soil fertility in a citrus orchard according to claim 2, characterized in that: The logic for performing ordinary kriging interpolation on soil physicochemical index data is as follows: Based on soil physicochemical index data from the grid center, a spatial autocorrelation model is constructed to generate spatially continuous functions of the physicochemical indices. Leave-one-out cross-validation was used to evaluate the interpolation results; The humidity, conductivity, and temperature data monitored by the sensor are normalized, a time covariance function is constructed, and a time continuous function is generated. , For continuous time variables; Reflectance correction and geometric registration were performed on UAV multispectral imagery, and vegetation indices and soil diagnostic indices were fused to generate a spectral continuity function. ; The three types of functions mentioned above are unified to the same coordinate system and resolution, and data integration is achieved through spatiotemporal alignment and outlier removal.

4. The method for characterizing soil fertility in a citrus orchard according to claim 3, characterized in that: The logic for extracting the spatial function values, temporal function values, and spectral function values ​​of each grid and combining them into a comprehensive fertility function using the entropy method is as follows: The target period for fruit enlargement in the citrus orchard is set as follows: Extract spatially continuous function values ​​for each grid. The values ​​of the time-continuous function and the spectral continuous function are in Integral average and Furthermore, it is dynamically weighted according to the entropy method and integrated into a comprehensive fertility function. : ； in, The weights are dynamically fused to reflect the differences in information content across spatial, temporal, and spectral dimensions. make This represents the output of the i-th grid in year b. For year b, calculate the 90th percentile of the output of all grids in that year. , The high-yield area sample set is defined as areas where the yield in each of the three consecutive years is no less than the 90th percentile of that year. Grid: And satisfy the high-yield area sample set Total number of samples , Grid index for high-yield areas; For each high-yield grid The dimensionality-reduced vector of the high-yield area is obtained by using the PCA function. , The three-dimensional Epanechnikov kernel was selected, and the bandwidth was chosen according to the Silverman rule of thumb. Then refer to the distribution density for: ； in, Given an input vector with an arbitrary grid, Corresponding to 3D volume scaling, For indicator functions, It is the Euclidean norm.

5. The method for characterizing soil fertility in a citrus orchard according to claim 4, characterized in that: The logic for performing principal component decomposition on the comprehensive fertility function and extracting the first three principal components to obtain the dimensionality-reduced vector is as follows: For each grid center During the full growing season of the target citrus orchard The observed comprehensive fertility function Notation: ； Principal component analysis of the function is performed, and the first three principal components are extracted to obtain the dimensionality-reduced vector. Define a dimension reduction vector Based on high fertility reference distribution Tukey depth; Measured by Tukey depth Similarity to the sample set of high-fertility, high-yield areas The formula used is: ； in, For inclusion The set of all closed half-spaces, that is, the half-spaces in the principal component space divided by the hyperplane. The reference is distributed in a closed half-space. The probability mass in; Introducing Mahalanobis distance as a fast approximation: ； in, , ; These are the mean and covariance of the reference distribution in high-yield areas, representing the spatial distribution center and morphology of the principal components; Mahaldo depth Mapped to soil fertility index : ； in, When the fertility at that location is most similar to the reference distribution of high-yield areas, the fertility is optimal; when... When the fertility at that location is significantly lower than the benchmark for high-yield areas, it indicates that the fertility at that location is significantly lower than that of high-yield areas.

6. The method for characterizing soil fertility in a citrus orchard according to claim 5, characterized in that: The logic behind classifying fractures into five levels using Jenks' natural fracture classification method is as follows: SFI values ​​for all grids Sort the SFI values ​​from smallest to largest; divide them into 5 categories, with the category boundaries as follows: Let o be the category index, and let the o-th category be: ； in, This refers to the index of the grid's position in the sorted sequence; The improved Jenks natural fracture method is used to classify fertility levels, and the objective function is defined. The weighted sum of intra-class tightness and inter-class separation: ； in, It is the mean of class o. The weighting is determined by the proportion of samples in class o to the total number of samples, making the intra-class tightness more focused on the intermediate stable samples. Let be the number of samples in class o. It serves as a balancing factor for inter-class separation. The simulated annealing algorithm is used to search for the optimal classification boundary, and the boundary is finally obtained. Then, the optimal five-level classification result was obtained. arrive ,satisfy arrive It is a highly fertile area and spatially continuous; The soil fertility classification rules for the target planting area are as follows: grade SFI range is This is a priority improvement zone; level SFI range is This area requires additional fertilization; grade SFI range is This is a routine management area; level SFI range is To optimize the maintenance zone; level SFI range is It is a high-yield protected area.

7. A soil fertility characterization system for citrus orchards, characterized in that: The system is used to perform a method for characterizing soil fertility in a citrus orchard as described in any one of claims 1-6, comprising: The data acquisition module is used to divide the target citrus orchard into a uniform sampling grid, set a sampling point at the center of each grid to collect soil samples at a depth of W, and conduct laboratory tests on the soil samples to obtain the soil physicochemical index values ​​at each point. Simultaneously, a three-parameter soil sensor is deployed at the center of each grid to monitor soil moisture, electrical conductivity and temperature data, and acquire multispectral remote sensing images covering the entire target citrus orchard. The data processing module is used to construct spatially continuous functions of soil physicochemical indicators at each point using ordinary Kriging interpolation, generate time-continuous functions by applying Gaussian process regression to sensor data, calculate normalized vegetation index and soil diagnostic index after atmospheric correction and geometric registration of remote sensing images, and then construct spectral continuous functions. The function calculation module is used to select soil samples from the top 10% of citrus production in each of the past three years as reference samples for high-yield areas. Based on these reference samples, a high fertility reference distribution with kernel density estimation is constructed. Spatial function values, temporal function values, and spectral function values ​​of the soil samples are extracted and combined into a comprehensive fertility function using the entropy method. The fertility assessment module is used to perform principal component decomposition on the comprehensive fertility function, extract the first three principal components to obtain a dimensionality-reduced vector, and calculate the depth value of the dimensionality-reduced vector relative to the high fertility reference distribution using the Tukey depth function to generate a soil fertility index. The Jenks natural fracture classification method is used to divide each grid into five levels according to the soil fertility index, where level 4 and level 5 represent high fertility areas and level 1 represents low fertility areas.

Images