Supplemental cultivated land soil fertility rapid evaluation device and method based on unmanned aerial vehicle multispectral imaging

Abstract

The present application relates to a kind of based on unmanned aerial vehicle multispectral imaging supplementary farmland soil fertility rapid evaluation device and method, belong to resource environment remote sensing and precision agriculture technical field.The method is aimed at solving the problem of low model calibration accuracy caused by insufficient representativeness of preset sampling points in the supplementary farmland with large soil spatial variation.The technical scheme includes: configuring unmanned aerial vehicle system, and collecting multispectral image by ground control terminal control;Real-time calculation of the spatial variation coefficient of near-infrared band reflectivity of image, when exceeding threshold, additional sampling point coordinates are dynamically generated and guided to collect spectrum in the field;Integrate basic and additional sampling point data, and construct soil reflectance calibration dataset;Using the dataset to calibrate the localized soil reflectance model, the soil organic matter and alkali-hydrolyzed nitrogen content are inverted, and the soil fertility spatial distribution map is generated.The method is mainly used to realize the rapid and accurate evaluation of supplementary farmland soil fertility.

Description

Technical Field

[0001] This invention belongs to the field of resource and environmental remote sensing and precision agriculture technology, specifically relating to a device and method for rapid assessment of supplementary arable land soil fertility based on UAV multispectral imaging. Background Technology

[0002] Existing methods for soil fertility assessment using remote sensing technology typically rely on pre-deploying a fixed number of soil sampling points within the target area and collecting their spectral data to establish an inversion model between remote sensing image reflectance and soil fertility parameters. However, supplementary cultivated land is often created by backfilling or modifying soil from different sources, and its spatial heterogeneity of soil properties is usually much higher than that of conventional cultivated land that has undergone long-term cultivation and management. Under such conditions of high spatial variability, if a fixed number of basic sampling points are still pre-set, the spatial distribution of these sampling points may not be sufficient or effective in capturing the complex spatial variation pattern of soil fertility. This leads to insufficient representativeness of the dataset used for model calibration in terms of the overall regional conditions, and the established model parameters are easily affected by local extreme values ​​or the variability characteristics of unsampled areas, thus significantly reducing the overall accuracy and reliability of the model when inverting key fertility indicators such as soil organic matter and available nitrogen across the entire supplementary cultivated land area. Since soil spatial variability patterns are often unknown and difficult to predict before flight operations, how to construct a calibration dataset that can adaptively cover highly variable areas and is sufficiently representative with limited and efficient sampling costs has always been a technical challenge in improving the accuracy of rapid assessment of soil fertility in supplemented farmland. Summary of the Invention

[0003] One object of the present invention is to solve at least the above-mentioned problems and to provide at least the advantages that will be described later.

[0004] Another objective of this invention is to provide a rapid assessment method for supplementary arable land soil fertility based on UAV multispectral imaging. This method can dynamically guide ground supplementary sampling by analyzing the spectral spatial variation of flight-acquired images in real time, thereby solving the problem of low model calibration accuracy caused by insufficient representativeness of preset sampling points in supplementary arable land with large soil spatial variation, and improving the accuracy of soil fertility parameter inversion.

[0005] To achieve these objectives and other advantages of the present invention, a rapid assessment method for supplementary arable land soil fertility based on UAV multispectral imaging is provided, comprising the following steps: The configuration includes a drone platform, a multispectral camera mounted on the platform, a positioning module, a data transmission radio, and a ground control terminal. The flight path is planned through the ground control terminal, and the UAV platform is controlled to fly and collect multispectral images of the farmland. At the same time, the positioning module records the center position coordinates of each image. Multispectral images are transmitted in real time to the ground control terminal via a data transmission radio. The ground control terminal calculates the spatial variation coefficient of reflectance in the near-infrared band for each image in real time. The spatial variation coefficient is the ratio of the standard deviation to the mean. When the spatial variation coefficient is greater than a preset threshold, the ground control terminal generates additional sample point coordinates based on the high variation area in the current image and sends them to the mobile terminal of the ground operator. The preset threshold is 0.15-0.30. Ground operators use handheld spectrometers to collect on-site spectral data based on coordinates and simultaneously transmit the coordinates and data back. Before the drone platform takes off, 10-30 basic soil sample points are set up in the supplementary cultivated land, and the field spectral data and location coordinates of all basic sample points are collected. Based on the position coordinates of all basic and additional sample points, the ground control terminal extracts the corresponding pixel reflectance from the multispectral image to form a pixel reflectance set. The set of pixel reflectance data is combined with all the corresponding field spectral data to construct a soil reflectance calibration dataset; The parameters of the preset soil reflectance model were calibrated using the soil reflectance calibration dataset to obtain the localized calibrated model. The localized and calibrated model was used to perform pixel-by-pixel calculations on all multispectral images to retrieve the content data of soil organic matter and soil available nitrogen. Based on soil organic matter and available nitrogen content data, a spatial distribution map of soil fertility for supplementary cultivated land is generated.

[0006] In assessing soil fertility in supplementary cultivated land, the complex origins and uneven backfilling of soils lead to significant spatial heterogeneity. Traditional methods based on fixed, pre-set sampling points often struggle to obtain sufficiently spatially representative calibration data, limiting the predictive accuracy of the established remote sensing inversion model in areas of local variation. To address this challenge, this method analyzes the spatial variation of near-infrared reflectance in multispectral images acquired during UAV flight in real time. When the variation exceeds a preset threshold, it automatically identifies areas with prominent spectral variations and dynamically generates the coordinates of additional sampling points to guide ground-based supplementary sampling. By combining these dynamically triggered additional sampling points with pre-deployed basic sampling points, a soil reflectance calibration dataset is constructed, enabling the data used for model calibration to adaptively cover highly variable areas in cultivated land. This approach effectively improves the representativeness of the calibration data for the entire assessment area, thereby enhancing the robustness and accuracy of the soil reflectance model subsequently calibrated based on this dataset. Ultimately, it achieves more reliable spatial inversion and mapping of organic matter and available nitrogen content in supplementary cultivated land soil.

[0007] Preferably, in the step of planning the flight route through the ground control terminal, the flight altitude of the UAV platform is set to 50-150 meters based on the preset accuracy of generating additional soil sample point coordinates and the spatial resolution of the multispectral camera. The forward overlap is set to 65%-80% and the lateral overlap is set to 55%-70% according to the requirements of subsequent real-time calculation of spatial variation coefficient. This ensures that the acquired multispectral images can meet the geometric positioning accuracy required to extract additional sample point coordinates from the images, and can also provide image data with sufficient redundancy for real-time calculation of the spatial variation coefficient of near-infrared reflectance of each image.

[0008] When planning the UAV flight path, the flight altitude setting needs to comprehensively consider the geometric positioning accuracy required for precise extraction of additional sampling point coordinates from the imagery, as well as the inherent spatial resolution of the multispectral camera. Therefore, the flight altitude is limited to between 50 and 150 meters. Simultaneously, to ensure real-time and reliable calculation of the spectral spatial variation coefficient of each image during flight, sufficient redundancy and continuity of the image data are required. Therefore, the forward overlap is set to 65% to 80%, and the lateral overlap is set to 55% to 70%. This series of collaboratively set flight parameters ensures that the acquired multispectral imagery data not only meets the requirements of high-precision geolocation but also provides a stable and sufficient data foundation for real-time spectral variation analysis, thereby supporting the smooth implementation of the entire rapid evaluation process and the accuracy of the final results.

[0009] Preferably, the step of the ground control terminal calculating the spatial variation coefficient of reflectance in the near-infrared band for each image in real time specifically includes: After receiving the current multispectral image, the ground control terminal first performs soil pixel determination on each pixel in the image. The determination criteria are that the ratio of the reflectance value of the pixel in the red band to the reflectance value in the near-infrared band is greater than 0.7, the reflectance value of the pixel in the green band is less than 0.4, and the reflectance value of the pixel in the blue band is less than 0.3. Pixels that simultaneously meet the above three criteria are marked as soil pixels; Based on the set of all pixels labeled as soil pixels, the reflectance value of each pixel in the near-infrared band is extracted. Calculate the average of the near-infrared reflectance values ​​for all near-infrared bands in this soil pixel set; Calculate the standard deviation of the near-infrared reflectance values ​​for all near-infrared bands in this soil pixel set; The ratio obtained by dividing the standard deviation by the mean is the spatial variation coefficient of the image used to trigger additional sampling; The subsequent step of generating additional sample point coordinates is performed only if the spatial variation coefficient calculated based on the soil pixel set is greater than a preset threshold.

[0010] When processing multispectral images transmitted by UAVs in real time and calculating the spatial variation coefficient of their near-infrared reflectance, directly using all pixels in the image for calculation would incorporate spectral information from non-soil features such as vegetation, water bodies, and shadows. This would result in the coefficient failing to purely represent the spectral variation characteristics of the soil itself, potentially leading to misinterpretations of dynamic sampling commands. Therefore, before calculation, a multi-band joint judgment rule is established based on the differences in spectral responses between soil and typical non-soil features in the red, green, blue, and near-infrared bands. Only pixels that simultaneously satisfy the conditions of a red-to-near-infrared reflectance ratio greater than 0.7, a green reflectance less than 0.4, and a blue reflectance less than 0.3 are identified as high-confidence soil pixels. Subsequently, the mean and standard deviation of the reflectance are calculated only based on this subset of soil pixels, and the ratio of the standard deviation to the mean is used as the final spatial variation coefficient. This approach effectively eliminates interference from non-soil factors by pre-screening data sources, enabling the calculated coefficient of variation to more accurately and specifically reflect the spectral differences within the soil. This provides a more reliable and accurate basis for subsequent judgments on whether to trigger dynamic supplementary sampling, enhancing the robustness and scientific nature of the entire assessment system.

[0011] Preferably, the step of the ground control terminal generating additional sample point coordinates based on highly variable areas in the current image specifically includes: The ground control terminal performs image block processing on the current multispectral image with a spatial variation coefficient greater than a preset threshold, dividing it into several rectangular sub-regions; For each rectangular sub-region, the local anomaly factor value of the reflectance of all pixels within it in the near-infrared band is calculated. The calculation of the local anomaly factor value is based on the k-nearest neighbor distance of the pixel reflectance within the rectangular sub-region, which is used to characterize the degree of deviation of the reflectance of a single pixel relative to the pixels in its local neighborhood. Within each rectangular sub-region, the calculated local anomaly factor values ​​are sorted in descending order, and the top N pixels with the highest local anomaly factor values ​​are selected, where N is an integer between 1 and 5. The pixel row and column coordinates of each selected pixel in the current multispectral image are combined with the center position coordinates of the image and the spatial resolution and imaging geometry model of the multispectral camera to calculate the corresponding geographic coordinates. The ground control terminal uses the calculated geographic coordinates as the coordinates of additional sample points.

[0012] When a drone identifies an image whose overall spectral variation exceeds a preset threshold and triggers a supplementary sampling mechanism, accurately determining the specific sampling location within that image becomes crucial. Relying solely on the overall variation coefficient fails to distinguish whether the variation is widespread or concentrated in certain local micro-regions, potentially leading to inaccurate sampling point placement. To address this, the system further divides the entire image into multiple rectangular sub-regions. Within each sub-region, based on the reflectance value of pixels in the near-infrared band, it calculates an index—the local anomaly factor—that quantifies the spectral deviation of each pixel from its local neighborhood. By selecting the pixels with the highest factor values ​​within each sub-region, the system can precisely locate anomaly points whose spectral characteristics most differ from their surroundings. Subsequently, the coordinates of these selected pixels are combined with the image's georeferenced information to calculate the actual geographic coordinates. This process deepens the understanding from identifying regional variation to locating specific anomalies, enabling dispatched sampling tasks to directly target the micro-locations with the most significant spectral variations and the greatest need for on-site verification. This significantly improves the targeting of supplementary sampling operations and the spatial representativeness of the collected data, laying a solid foundation for building a high-quality calibration dataset.

[0013] Preferably, in the step of calculating the local anomaly factor value of the reflectance in the near-infrared band for all pixels within each rectangular sub-region, the k value on which the k-nearest neighbor distance for calculating the local anomaly factor value depends is dynamically determined for any rectangular sub-region as follows: The ground control terminal first calculates the average and standard deviation of the near-infrared reflectance of all pixels within the rectangular sub-region; Based on the mean and standard deviation, a first reflectance threshold and a second reflectance threshold are set. The first reflectance threshold is the mean plus one standard deviation, and the second reflectance threshold is the mean minus one standard deviation. The number of pixels within the rectangular sub-region whose reflectance values ​​fall between the first reflectance threshold and the second reflectance threshold is counted and recorded as the core pixel count of the sub-region. The number of core pixels in the rectangular sub-region is compared with the preset minimum nearest neighbor number K. min And the number of maximum nearest neighbors K max Compare; K min For 5, K max It is an integer corresponding to 20% of the total number of cells in the rectangular sub-region, and does not exceed 50; If the number of core pixels is less than or equal to K min Then, the k value used to calculate the local anomaly factor value for this rectangular sub-region is determined to be K. min ; If the number of core pixels is greater than or equal to K maxThen, the k value used to calculate the local anomaly factor value for this rectangular sub-region is determined to be K. max ; If the core pixel count is between K min With K max Between these values, the k value used to calculate the local anomaly factor value for the rectangular sub-region is determined to be equal to the number of core pixels; After the value of k is determined, for each pixel in the rectangular sub-region, the ground control terminal calculates the Euclidean distance from it to all other pixels in the sub-region in terms of near-infrared reflectivity. Based on the calculated Euclidean distance, find the k nearest neighbors of the pixel as its k-nearest neighbors, and calculate the local reachability density of the pixel based on the distances of the k nearest neighbors; The local anomaly factor value of the pixel is calculated based on the local reachability density and the local reachability density of the pixel's k nearest neighbors.

[0014] When using the local anomaly factor algorithm to perform detailed analysis of image sub-regions to locate spectral anomalies, using a uniform, fixed parameter to define the range of the local neighborhood is difficult to adapt to the varying soil reflectance distributions within different sub-regions. To ensure that the algorithm can reasonably define the local background and accurately identify true anomalous pixels under different conditions such as concentrated or dispersed reflectance distributions, the system independently calculates the average and standard deviation of the near-infrared reflectance for each rectangular sub-region to be analyzed. Based on this statistical characteristic, the system sets a core interval reflecting the central trend of the data and counts the number of pixels falling within this interval. This number is compared with a reasonable range preset based on experience, and finally, an optimal neighborhood size parameter is dynamically determined for the sub-region. This parameter determination mechanism allows the concept of "local" to flexibly adapt to the spectral variation scale within the current sub-region. Based on this dynamically determined parameter, the algorithm then calculates the deviation of each pixel from its adaptive neighborhood, thereby achieving accurate capture of local spectral anomaly patterns. This method improves the adaptability of anomaly identification, avoids misjudgment or omission due to parameter rigidity, and makes the final generated additional sampling point coordinates more representative of those micro-regions with real and significant spectral variations.

[0015] Preferably, the step of the ground control terminal extracting the corresponding pixel reflectance from the multispectral image based on the position coordinates of all basic sample points and supplementary sample points to form a pixel reflectance set specifically includes: For the location coordinates of each sample point, the ground control terminal performs the following operations: Based on the positioning accuracy parameters of the positioning module and the spatial resolution of the multispectral camera, the expected position uncertainty radius of the sample point is calculated. The expected position uncertainty radius is equal to the sum of the positioning accuracy error value and the three pixel sizes to cover the random errors that may exist in the positioning system and the residual errors of image geometric correction, ensuring that the real sampling point position falls within the search neighborhood with a high probability. Using the coordinates of the sample point as the center and the radius of the uncertainty of the expected location as the radius, a circular neighborhood is delineated on the corresponding multispectral image; Extract all pixels located within a circular neighborhood from the multispectral image; Based on the soil pixel determination criteria, each extracted pixel is determined and the pixels belonging to soil pixels are selected; the soil pixel determination criteria are based on the reflectance values ​​of the pixel in the red band, near-infrared band, green band and blue band. Determine the number of soil pixels selected: If the number is greater than or equal to 3, the median of the reflectance of these soil pixels in the near-infrared band is calculated, and this median is used as the pixel reflectance value representing the sample point. If the number is less than 3, keep the center coordinates of the circular neighborhood unchanged, expand the radius of the uncertain location to 1.5 times its original value, redefine the circular neighborhood and repeat the steps of pixel extraction, soil pixel determination and quantity judgment; this radius expansion step is repeated a maximum of 3 times. If the number of soil pixels selected is still less than 3 after the radius is expanded 3 times, then the sample point is marked as an invalid sample point. For all sample points that are not marked as invalid sample points, their corresponding pixel reflectance values ​​are collected to form a pixel reflectance set.

[0016] When extracting spectral reflectance from UAV imagery based on preset and dynamically supplemented sample point coordinates, directly obtaining the value of a single pixel may lack representativeness due to positioning deviations, residual errors from image geometric correction, or uneven soil mixing within the pixel. To address this challenge, the system constructs a dynamic search strategy for each sample point. First, based on the UAV's positioning accuracy parameters and camera resolution, an initial search radius that integrates positioning errors and multiple pixel scales is calculated, and a circular area is delineated on the image using this radius. After extracting all pixels from this area, pixels belonging to soil are selected according to specific spectral rules. When the number of selected soil pixels reaches a certain requirement, the median is taken as the representative spectral value of the sample point; if the number is insufficient, the search range is gradually expanded proportionally, and this process is repeated to provide multiple opportunities to obtain valid data. If, ultimately, not enough valid pixels can be obtained, the sample point is carefully excluded. This method significantly reduces random errors introduced by misalignment of a single coordinate point or impure pixels by constructing a buffer zone that tolerates spatial uncertainty and using robust statistics to represent the point spectrum, thereby improving the reliability and accuracy of registration between spectral data extracted from images and field sampling points.

[0017] Preferably, the step of combining the pixel reflectance set with all corresponding field spectral data to construct a soil reflectance calibration dataset specifically includes: For each pixel reflectance value in the pixel reflectance set and its paired field spectral data, the ground control terminal performs the following operations to construct a valid data pair: Obtain the acquisition timestamp of the multispectral image corresponding to the pixel reflectance value, as well as the timestamp generated when the field spectral data is acquired by the handheld spectrometer; Calculate the absolute difference between two timestamps. If the difference is less than or equal to 5 minutes, mark the data pair as a valid time-synchronized data pair. For each time-synchronized valid data pair, acquire the solar elevation angle data recorded by the UAV platform at the corresponding multispectral image acquisition time, and the solar elevation angle data recorded by the mobile terminal at the on-site spectral data acquisition time; Calculate the difference between two solar altitude angle data. If the difference is less than or equal to 5 degrees, then mark the data pair as a matching data pair under illumination conditions. For each pair of illumination condition matching data, the ground control terminal performs the following operations to generate calibration data: The average reflectance value of the field spectral data within the corresponding spectral range of the near-infrared band of the multispectral camera is extracted and used as the field reference reflectance. Read the corresponding pixel reflectance value of this data and use it as the image reflectance. The on-site reference reflectance and the image-extracted reflectance are input into a preset linear normalization model; the parameters of the linear normalization model are determined by simultaneously measuring and fitting the same standard diffuse reflectance reference plate using the UAV platform and a handheld spectrometer before each flight mission. A linear normalization model is applied to correct the image reflectance, resulting in a normalized image reflectance value. The normalized image reflectance value is paired with its corresponding ground reference reflectance value to form a tuple for calibration. All the pairs obtained from the above steps are summarized to form a soil reflectance calibration dataset.

[0018] Constructing a soil reflectance calibration dataset for model calibration presents challenges when directly merging reflectance data extracted from UAV imagery with spectral data measured by ground-based handheld devices. This is because the acquisition time, lighting conditions, and the radiometric response characteristics of the sensors may differ. To ensure physical comparability between paired data, the system employs a rigorous screening and calibration process. First, it requires a high degree of synchronization between image acquisition and ground measurement, with the time difference controlled within an extremely short range. Simultaneously, the solar illumination angle at the two acquisition times must be essentially consistent. For data pairs that pass the spatiotemporal consistency test, the system does not directly use the original image reflectance but instead inputs it along with the corresponding ground reference reflectance into a pre-calibrated transformation model. This model, established by synchronously measuring both sensors using the same standard reference before flight, corrects the reflectance values ​​acquired by the image sensors to a physical scale consistent with the ground-based handheld spectrometer reference. Through this rigorous spatiotemporal matching and spectral normalization between sensors, the final calibration dataset exhibits high internal consistency and physical uniformity for each data point, laying a solid foundation of data quality for subsequent reliable model parameter calibration.

[0019] Preferably, the parameters of a pre-defined soil reflectance model are calibrated using a soil reflectance calibration dataset to obtain a locally calibrated model, specifically including: The ground control terminal uses all normalized image reflectance values ​​in the soil reflectance calibration dataset as model input data and all corresponding field reference reflectance values ​​as expected model output data. The preset soil reflectance model contains a set of model parameters to be determined. The model parameters are optimized and calibrated using a weighted iterative re-estimation method, which involves performing the following steps through multiple iterative cycles: a) Based on the model parameters of the current iteration, calculate the deviation between the model-predicted reflectance value and its field reference reflectance value for each data point in the calibration dataset; b) Based on the deviations of all data points obtained in this iteration, calculate the median of the absolute values ​​of all deviations, which serves as a measure of the dispersion of the current data deviations; c) Based on the calculated deviation dispersion scale and applying the preset robust weight calculation rules, calculate a specific weight value for each data point in the calibration dataset; the weight calculation rules ensure that individual data points with deviation patterns significantly different from the majority of data points are assigned lower weights. d) Using the weights calculated for all data points in this iteration, construct a weighted optimization objective that aims to minimize the sum of squared weighted prediction biases for all data points; e) An iterative algorithm suitable for nonlinear parameter optimization is adopted to minimize the weighted optimization objective, and the model parameters are updated once to obtain a new set of model parameters; f) Use the updated model parameters as the starting point for the next iteration loop, and repeat steps a) to e). Set iteration termination conditions, including the change between model parameters obtained from two adjacent iterations being less than a preset small threshold, or the total number of iterations reaching a preset upper limit; When any iteration termination condition is met, the iteration stops, and the model parameters obtained from the last iteration are used as the final calibration result. Substitute the final calibration results into the preset soil reflectance model to obtain the localized soil reflectance model with completed parameter calibration.

[0020] When calibrating local parameters for a soil reflectance model, the dataset used for calibration may contain a few data points that deviate from the overall trend due to soil spatial heterogeneity or measurement randomness. If a conventional optimization method assigns equal importance to all sample points, these local outliers can excessively influence the final determination of model parameters, leading to a bias in the calibrated model's reflection of the overall regional patterns. Therefore, the system employs an iterative parameter optimization strategy. In each iteration, the system first calculates the prediction bias for each calibration data point based on the current temporary model parameters. Then, based on the bias distribution of all data points, the system dynamically calculates and assigns a weight value to each point. This weighting rule aims to reduce the influence of data points that significantly deviate from the mainstream bias pattern. Using this dynamic weighting, the system constructs a new weighted optimization objective and solves for an updated set of model parameters under this objective. This process is repeated iteratively, with model parameters and data point weights adjusting and optimizing each other during iteration. When the parameter updates stabilize, the iteration stops, and the final model parameters are output. This mechanism enables the parameter optimization process to automatically identify and suppress interference from potential outlier data, ensuring that the calibrated model parameters more robustly capture the inherent trends of the data subject, and improving the model's generalization ability and prediction reliability on new data.

[0021] Preferably, the localized and calibrated model is used to perform pixel-by-pixel calculations on all multispectral images to retrieve the content data of soil organic matter and soil available nitrogen, specifically including: The ground control terminal takes the reflectance value of each pixel in all multispectral images, which has been corrected in each band of the multispectral camera, as the input vector and inputs it into the locally calibrated soil reflectance model pixel by pixel. For each pixel, the model calculates and outputs the predicted values ​​of soil organic matter content and soil available nitrogen content. For each soil organic matter content prediction data layer and soil available nitrogen content prediction data layer obtained from the inversion of a multispectral image, the ground control terminal performs the following post-processing steps: First, identify and label spatially isolated outlier pixels in the prediction data layer. The criteria for identifying spatially isolated outlier pixels is: the absolute difference between the predicted value of the pixel itself and the median of the predicted values ​​of the other 8 pixels in a 3x3 pixel window centered on the pixel exceeds twice the standard deviation of the predicted values ​​of all pixels in the window. Secondly, for each spatially isolated anomalous cell that is marked, its predicted value is replaced by the median of the predicted values ​​of all unmarked anomalous cells in its 3x3 cell neighborhood, excluding itself. Next, edge-aware spatial smoothing is performed on the prediction data layer. The specific steps of this process include: dividing the prediction data layer into several non-overlapping 32x32 pixel square blocks; for each pixel in each block, calculating the spatial weighted average of its predicted value with that of all other pixels in the block, with the weights calculated according to the inverse proportional function of the spatial Euclidean distance between pixels, and the closer the distance, the higher the weight; replacing the original predicted value of the pixel with the calculated spatial weighted average, but retaining areas where the gradient change at the edges between blocks is greater than a preset threshold without smoothing; The soil organic matter content prediction data layer and the soil available nitrogen content prediction data layer, which have undergone the above post-processing, are used as the final soil organic matter content data and soil available nitrogen content data obtained from the inversion, respectively.

[0022] When optimizing soil reflectance model parameters using calibration datasets, a few data points may deviate from the main distribution due to local soil anomalies or measurement errors. If a traditional equal-weight fitting method is used, these deviations can excessively influence the determination of model parameters, leading to a decrease in the representativeness of the final model for the overall region. To address this, a dynamic weighting iterative optimization strategy is adopted. In each iteration, this strategy automatically calculates and assigns a weight to each data point based on the current model's predicted residual distribution for all data points, giving lower weights to points that significantly deviate from the mainstream trend. The model parameters are updated using this dynamic weighting, and the residuals are re-evaluated and the weights adjusted based on the new parameters. This process is repeated until the model parameters stabilize. This process allows parameter optimization to focus on data points reflecting the main trend while effectively reducing the interference of individual outliers, resulting in a more robust and generalized localized model, improving the accuracy and reliability of the model in practical inversion applications.

[0023] A rapid assessment device for supplementary arable land soil fertility based on UAV multispectral imaging, used to implement the method described in this invention, comprising: Unmanned aerial vehicle (UAV) platform; A multispectral camera, mounted on a drone platform, is used to collect multispectral images of supplementary farmland. The positioning module, located on the drone platform, is used to record the center coordinates of each multispectral image; Data radio, set up on the drone platform, is used to transmit multispectral images to the ground in real time; The ground control terminal, which communicates with the data radio, is configured as follows: Plan the flight path of the unmanned aerial vehicle (UAV) platform and control its flight; It receives multispectral images and calculates the spatial variation coefficient of reflectance in the near-infrared band for each image in real time. When the spatial variation coefficient is greater than a preset threshold, additional sample point coordinates are generated based on the high variation region in the current image. Based on the position coordinates of all basic and additional sample points, the corresponding pixel reflectance is extracted from the multispectral image to form a pixel reflectance set. The set of pixel reflectance data is combined with all the corresponding field spectral data to construct a soil reflectance calibration dataset; The parameters of the preset soil reflectance model were calibrated using the soil reflectance calibration dataset to obtain the localized calibrated model. The localized and calibrated model was used to perform pixel-by-pixel calculations on all multispectral images to retrieve the content data of soil organic matter and soil available nitrogen. A spatial distribution map of soil fertility in the supplementary cultivated land is generated based on the content data; A handheld spectrometer, operated by a ground operator, is used to collect field spectral data based on the coordinates of additional sample points and to simultaneously transmit the coordinates and data back. The mobile terminal communicates with the ground control terminal and the handheld spectrometer to receive additional sample point coordinates and guide the ground operator, and to transmit back the field spectral data and its acquisition location coordinates.

[0024] The present invention has at least the following beneficial effects: First, the core innovation of this invention lies in realizing a paradigm shift from preset fixed sampling to dynamic adaptive sampling. By analyzing the spectral variation of images in real time during UAV flight and automatically identifying highly variable areas to dynamically generate supplementary sampling points, the calibration dataset constructed by the system can proactively cover areas with drastic spatial variations in soil properties. This fundamentally solves the problem of low model calibration accuracy caused by insufficient representativeness of preset sampling points in highly heterogeneous environments such as supplementary cultivated land, thus significantly improving the overall accuracy and regional adaptability of subsequent soil organic matter and alkaline nitrogen inversion models.

[0025] Secondly, the method of this invention comprehensively ensures the quality and reliability of each stage from data acquisition and processing to model application through a series of meticulous data processing steps. Specifically, pre-flight parameter co-planning ensures the geometric and spectral usability of image data; rigorous soil pixel screening and precise localization based on local anomalies ensure the targeting and effectiveness of sampling commands; a highly consistent calibration dataset is constructed by extracting representative spectral values ​​from a buffer zone tolerant of localization uncertainties and implementing rigorous spatiotemporal synchronization and sensor normalization correction; finally, a spatially continuous and reliable soil fertility distribution map is obtained by employing a robust iterative weighted algorithm for model calibration and edge-preserving noise reduction of the inversion results. This entire process is interconnected and works together to greatly enhance the robustness of the entire technical system and the credibility of the output results.

[0026] Third, this method improves accuracy while also optimizing the efficiency and cost-effectiveness of the assessment work. The UAV platform enables rapid coverage and data collection over a large area, while the intelligent sampling mode with real-time analysis and dynamic guidance allows limited ground manpower and time resources to be precisely targeted at the key variation points that most need verification, avoiding the blindness and redundant work that may exist in traditional uniform grid sampling. This air-ground collaborative and intelligently guided operation mode achieves efficient allocation of manpower and time costs while ensuring assessment accuracy, providing a practical and feasible technical solution for the rapid and large-scale assessment of farmland soil fertility.

[0027] Other advantages, objectives and features of the present invention will become apparent in part from the following description, and in part from those skilled in the art through study and practice of the invention. Detailed Implementation

[0028] The present invention will now be described in further detail so that those skilled in the art can implement it based on the description.

[0029] It should be understood that terms such as “having,” “comprising,” and “including” as used herein do not exclude the presence or addition of one or more other elements or combinations thereof.

[0030] In existing technologies, a fixed number of soil sample points are typically pre-planned and deployed before the assessment. After the UAV acquires multispectral imagery, only the spectral data from these pre-determined locations are matched with the image reflectance to establish and calibrate the soil property inversion model. This method assumes that the sampling points can uniformly represent the entire area and does not consider new information acquired during flight.

[0031] However, due to the complex soil sources and uneven backfilling, the fertility properties of supplemented farmland exhibit extremely high spatial heterogeneity and unpredictable patterns. The aforementioned method of using fixed, pre-set sampling points often fails to effectively capture those unknown, locally volatile areas. This directly results in insufficient representativeness of the dataset used for model calibration regarding the overall regional conditions, making the established model parameters susceptible to the influence of unsampled extreme values. Consequently, this significantly reduces the overall accuracy and reliability of the model when retrieving key indicators such as soil organic matter and available nitrogen across the entire region.

[0032] This invention provides a rapid assessment method for supplementary arable land soil fertility based on UAV multispectral imaging. First, operators configure an integrated UAV system, including a UAV platform, an onboard multispectral camera, a positioning module, a data transmission radio, and a ground control terminal. Before the flight operation begins, operators, based on experience, pre-select 10-30 basic soil sample points evenly within the supplementary arable land area to be assessed, and use a handheld spectrometer to collect the on-site spectral data of these points, while simultaneously recording their precise geographical coordinates. After preparation, the ground control terminal plans the UAV's automatic flight path, sets appropriate flight altitude, heading, and lateral overlap, and then controls the UAV to take off and perform the image acquisition task. During flight, the positioning module synchronously records the center point coordinates of each acquired multispectral image. The image data is transmitted back to the ground control terminal in real time via the data transmission radio. Upon receiving each new image, the ground control terminal immediately processes it. It first judges each pixel in the image based on specific threshold relationships of reflectance in the red, green, blue, and near-infrared bands, filtering out pixels belonging to bare soil. Next, based solely on these soil pixels, the average and standard deviation of their reflectance in the near-infrared band are calculated, and the ratio of the standard deviation to the mean is used as the spatial variability coefficient of the image. The system presets a variability coefficient threshold. This threshold can be fine-tuned based on historical soil data for supplementary cultivated land or a small amount of pre-flight image data to achieve an optimal balance between sampling cost and model accuracy. When the calculated variability coefficient exceeds this threshold, the ground control terminal determines that there is significant spectral spatial variability in the current image coverage area. At this point, it automatically performs further analysis on the image, using block segmentation and local anomaly detection algorithms to calculate the locations of several specific pixels with the most anomalous spectral characteristics. The row and column coordinates of these pixels are then combined with image geographic information to calculate actual geographic coordinates, generating a set of additional sample point coordinates. These coordinates are transmitted in real-time to the mobile terminal held by the ground operator via a communication link. Following the navigation guidance of the mobile terminal, the ground operator quickly proceeds to the locations indicated by these coordinates, uses a handheld spectrometer to conduct on-site spectral measurements, and synchronously transmits the measurement data and confirmed location coordinates back. After the UAV completed its flight across the entire area and collected in-situ spectral data from all base points and dynamically triggered additional points, the ground control terminal began data integration and model building. Based on the location coordinates of all sample points, it robustly extracted the image reflectance values ​​representing each point from the corresponding multispectral images using a circular neighborhood search method that considers positioning uncertainty, forming a set of pixel reflectance values. Subsequently, the system rigorously paired these image-extracted reflectance values ​​with the in-situ measured spectral data, retaining only data pairs that highly matched the acquisition time and lighting conditions. It then used a sensor normalization model calibrated before flight to correct the image reflectance, ultimately constructing a high-quality soil reflectance calibration dataset.Using this dataset, the system optimizes and calibrates the parameters of the preset soil fertility inversion model to obtain a customized model suitable for the current assessment area. Finally, the calibrated model is applied to perform pixel-by-pixel calculations on all multispectral images of the entire area to invert the content of soil organic matter and available nitrogen, and to generate an intuitive spatial distribution map of soil fertility.

[0033] Through the above-described embodiments, this invention can perceive the spatial variability of soil spectra in real time during flight operations and dynamically guide ground sampling forces to precisely supplement the key locations with the most dramatic variability. This allows the final dataset used for model calibration to adaptively cover highly heterogeneous areas in cultivated land. This effectively overcomes the inherent deficiency of insufficient representativeness in fixed sampling, thereby significantly improving the accuracy of model calibration and the reliability and accuracy of the final spatial inversion results of soil fertility without significantly increasing sampling costs. The entire process achieves a closed loop from data acquisition and intelligent analysis to decision feedback, improving the overall efficiency and intelligence level of the assessment operation.

[0034] In current technologies, when planning UAV flight paths to acquire multispectral imagery, the flight altitude and image overlap are typically set based on general surveying standards or empirical values. For example, a high overlap rate is set to ensure map formation, or a low flight altitude is set to pursue efficiency. This type of setting is often an isolated step with the primary goal of obtaining complete image coverage, without fully considering the specific data quality requirements of subsequent analysis tasks (such as real-time spectral variation calculation and accurate point extraction).

[0035] However, in the dynamic sampling and rapid evaluation system constructed in this invention, the setting of flight parameters faces a dual challenge: on the one hand, it is necessary to accurately extract the coordinates of dynamically triggered additional sampling points from the images, which requires the images to have sufficient geometric positioning accuracy; on the other hand, it is necessary to calculate the spatial variation coefficient of each image in real time and reliably during flight, which requires the image data to have sufficient redundancy and continuity to resist local interference and ensure computational stability. If general parameters are simply adopted, the spatial resolution of the images may be insufficient to support centimeter-level point calculation, or insufficient overlap may result in an excessively small effective data area and distorted variation analysis during real-time calculation, thereby directly affecting the reliability of the entire dynamic sampling decision chain.

[0036] To address the aforementioned issues, furthermore, in the implementation of this invention, flight path planning is a meticulously calculated, end-to-end collaborative design process. During the planning phase, the ground control terminal comprehensively analyzes the inherent physical parameters of the multispectral camera and the data requirements of the entire method for subsequent steps. First, to meet the requirement of accurately calculating the geographic coordinates of additional sample points from the imagery, the system will, based on the required planar positioning accuracy and in conjunction with parameters such as the multispectral camera pixel size and lens focal length, reverse-calculate the allowable flight altitude range of the UAV platform. Specifically, the system will ensure that the ground size corresponding to a single pixel in the image acquired at this altitude meets the accuracy tolerance for calculating the coordinates of the additional sampling points, ultimately setting the flight altitude to a suitable value between 50 and 150 meters. Second, to ensure robust real-time calculation of the spatial variation coefficient for each newly acquired image, the system imposes higher requirements on the image overlap. Considering that the calculations require a sufficient number of pure soil pixels and that potentially distorted areas at image edges need to be excluded, the system sets high forward and lateral overlap, for example, forward overlap between 65% and 80% and lateral overlap between 55% and 70%. This ensures that any point within the flight area is captured by multiple images from different angles, providing data redundancy for real-time calculations. Simultaneously, the central area of ​​a single image is large enough to provide sufficient pixels for soil screening and statistical calculations, avoiding the problem of excessive randomness or failure in coefficient of variation calculations due to an excessively small effective area.

[0037] Through the aforementioned collaborative settings, the acquired multispectral imagery data possesses dual safeguards from the source: it retains sufficient geometric detail to support accurate subsequent geographic coordinate inversion, while also ensuring ample data redundancy and integrity to guarantee the reliability of real-time spectral analysis. This targeted optimization of flight parameters does not pursue flight efficiency or image coverage in isolation, but rather closely serves the core logic of subsequent real-time analysis and dynamic sampling. This lays a solid data foundation for the efficient and accurate operation of the entire evaluation method, enabling the raw data collected by the UAV to support a series of subsequent intelligent processing and decision-making steps to the greatest extent possible.

[0038] In existing technologies, when assessing soil fertility using UAV multispectral imagery, the statistical characteristics (such as mean and standard deviation) of near-infrared reflectance are typically calculated directly using all pixels in the entire image. Based on this, the spatial variation coefficient is then calculated as a basis for determining whether additional sampling is necessary. This method does not distinguish between different types of land features in the image and assumes that all pixels can represent the soil.

[0039] However, supplementary farmland often contains sparse vegetation, residual crops, puddles, temporary debris, or shaded areas, whose spectral characteristics differ significantly from those of the soil. Including these non-soil pixels in the calculation results in a coefficient of variation that not only reflects the soil's own spectral variation but also includes spectral interference from other features, making it impossible to accurately and purely represent the soil's internal variability. Triggering dynamic sampling based on such coefficients may lead the system to misjudge areas of homogeneous soil as having high variability and trigger redundant sampling, or to miss sampling in areas of significant actual soil variability due to non-soil pixels flattening statistical values. This reduces the accuracy of the entire dynamic sampling mechanism and the effectiveness of subsequent calibration datasets.

[0040] To address the aforementioned issues, this invention further proposes that after the UAV flies along a planned route and transmits multispectral imagery in real time, the ground control terminal does not immediately perform variation calculations on all pixels of the entire image. Instead, it first executes a rigorous soil pixel screening process. Specifically, for each newly received multispectral image, the ground control terminal checks the reflectance values ​​of each pixel in the red, green, blue, and near-infrared bands. The system is configured with the following judgment logic: a pixel is classified as a soil pixel only if it simultaneously meets three conditions. First, the ratio of the pixel's reflectance value in the red band to its reflectance value in the near-infrared band must be greater than 0.7; this ratio helps distinguish between vegetation and soil. Second, the pixel's reflectance value in the green band must be less than 0.4. Third, the pixel's reflectance value in the blue band must be less than 0.3. The latter two conditions together exclude interference from water bodies, highly reflective debris, and certain shadows. Only pixels that fully meet all three conditions are marked by the system and included in a dedicated set for subsequent calculations.

[0041] After soil pixel screening, the ground control terminal performs statistical calculations only based on this subset of soil pixels. It extracts the reflectance value of each pixel in the near-infrared band, calculates the average of these values, and then calculates their standard deviation. Finally, it divides the standard deviation by the average to obtain the spatial variation coefficient specifically used to trigger additional sampling for that image. The subsequent judgment logic remains consistent with the previous one: the system will only execute the subsequent step of generating additional sample point coordinates if the spatial variation coefficient calculated based on pure soil pixels is greater than a preset threshold.

[0042] Through the above implementation methods, this invention introduces a ground feature spectral filtering mechanism at the initial stage of calculating the spatial variation coefficient, ensuring that each data point involved in the statistics has high-confidence soil spectral properties. This allows the final spatial variation coefficient to more accurately and purely reflect the spectral heterogeneity of the supplemented cultivated land soil itself, thereby providing a reliable basis for dynamic sampling triggering decisions, effectively avoiding false triggering or missed triggering caused by non-soil ground feature contamination, and improving the robustness of the entire assessment process and the reliability of the results.

[0043] In existing technologies, when UAV image analysis indicates that a certain area has high overall spectral variation and requires supplementary sampling, the common practice is to randomly place sampling points within the image range or select several locations as additional sampling points based on empirical rules (such as grids). This approach only utilizes the macroscopic judgment of high overall regional variation without further exploring the specific distribution pattern of variation within the region, and cannot identify whether the variation is widespread and uniform or concentrated in certain specific local patches or anomalies.

[0044] This leads to a significant problem: additional sampling points may be placed in regions where overall spectral variability is high but local features are not prominent, while micro-hotspots that are truly spectrally abnormal and extremely valuable for model calibration may be missed. The sampling action lacks focus, and the obtained data may fail to effectively capture the key local information that leads to increased overall variability, thus weakening the actual role of supplementary sampling in improving the representativeness of the calibration dataset.

[0045] To address the aforementioned issues, this invention, after implementing the steps described above—namely, the ground control terminal calculating the spatial variation coefficient of a certain image based on soil pixels in real time and determining that it exceeds a preset threshold—does not simply assign one or a few general sampling locations to the entire image. Instead, the system immediately performs a more refined local analysis on the image that triggered the alarm. First, the ground control terminal automatically divides the multispectral image into multiple regularly sized rectangular sub-regions, each containing a certain number of pixels. Next, the system focuses on the near-infrared band and independently calculates a specific index—the local anomaly factor value—for each pixel within each rectangular sub-region. The calculation principle of this value is to quantify the degree of anomaly of a pixel relative to its immediate local background by measuring the average deviation in reflectance between a pixel and its nearest neighbors. After calculating the local anomaly factor values ​​for all pixels within each sub-region, the system sorts these values ​​from high to low and selects the top few pixels, typically one to five. These pixels are considered the local anomaly points within the sub-region whose spectral characteristics are most different from their surrounding environment.

[0046] Subsequently, the system needs to convert the locations of these outliers in the images into usable geographic coordinates. The ground control terminal, combining the known geographic coordinates of the image's center point, the spatial resolution parameters of the multispectral camera, and the camera's imaging geometry model, precisely calculates the row and column numbers of the selected outlier pixels in the image into their corresponding latitude and longitude coordinates. This series of coordinates, automatically generated by the system, becomes the additional sample point coordinates sent to the ground operator.

[0047] Through the above implementation methods, after identifying highly variable images requiring supplemental sampling, this invention does not remain at the overall regional level but delves into the image's interior, calculating local anomaly factors to locate the most anomalous microscopic locations within their respective small neighborhoods. This allows the assigned sampling tasks to directly target the focal points of spectral variation, significantly improving the targeting and accuracy of supplemental sampling. The collected field spectral data can more effectively represent the key local features leading to increased overall variation, thereby significantly enhancing the ability of the subsequently constructed calibration dataset to characterize complex spatial variation patterns and laying a solid foundation for obtaining a more accurate localized model.

[0048] When using the local anomaly factor algorithm to identify spectral anomalies within image sub-regions, existing techniques typically set a fixed neighborhood size parameter for the algorithm. This means that throughout the analysis, regardless of whether the reflectance values ​​within each sub-region are concentrated or dispersed, the same number of nearest neighbor pixels are used to define the local background of each pixel. This approach is based on the simplified assumption that the scale of local variation is similar across all regions.

[0049] However, the distribution characteristics of soil reflectance vary significantly across different plots within the supplementary cultivated land. Some sub-regions may exhibit highly concentrated reflectance values ​​with minimal variation, while others may show highly dispersed values. If a fixed neighborhood size parameter is used, in areas with concentrated reflectance, an excessively large neighborhood might misclassify normal, continuous, and gradual variations as isolated outliers; conversely, in areas with dispersed reflectance, an excessively small neighborhood might fail to capture genuine local anomaly patterns, leading to missed detections. This mismatch between the parameter and data characteristics directly reduces the accuracy of the local anomaly factor algorithm in locating spectral anomalies, thus affecting the quality of subsequently generated additional sampling point coordinates.

[0050] To address the aforementioned issues, this invention further introduces a dynamic adaptive mechanism for determining the neighborhood size when performing image segmentation and preparing to calculate the local anomaly factor values ​​of pixels within each rectangular sub-region. For any rectangular sub-region currently being processed, the ground control terminal first calculates the average and standard deviation of the near-infrared reflectance of all pixels within that region. Then, the system defines a core reflectance interval based on the average plus or minus one standard deviation, and counts the number of pixels whose reflectance values ​​fall within this interval; this number is referred to as the core pixel count of the sub-region. This core pixel count reflects the concentration of spectral values ​​of most pixels in the region.

[0051] Subsequently, the system compares this core pixel count with a preset reasonable range of neighborhood numbers. The lower limit of this range ensures that there are enough nearest neighbors for statistics, while the upper limit prevents the neighborhood from becoming too large and losing local significance. Specifically, if the core pixel count is lower than or equal to a preset minimum value, this minimum value is used as the neighborhood size parameter for calculating the local anomaly factor in this sub-region; if the core pixel count is higher than or equal to a preset maximum value, this maximum value is used; if the core pixel count is between the minimum and maximum values, this core pixel count is directly used as the parameter. After the parameter is determined, the system calculates the Euclidean distance in reflectance from each pixel to other pixels within the sub-region based on this dynamically determined neighborhood size, finds its k nearest neighbors, and finally completes the calculation of the local anomaly factor value for each pixel.

[0052] Through the above method, this invention enables the neighborhood size used to define local extents to be no longer a rigid, fixed value, but rather to be adaptively adjusted according to the actual density of reflectance data distribution within each image sub-region. Smaller neighborhoods are used in concentrated areas of the dataset to maintain sensitivity to subtle local anomalies, while larger neighborhoods are used in dispersed areas to avoid misclassifying widely distributed continuous variations as anomalies. This mechanism significantly improves the adaptability and recognition accuracy of the local anomaly factor algorithm under different soil spectral distribution characteristics, thereby ensuring that the located additional sampling points are indeed key locations where spectral features are truly distinctive, providing more reliable target information for constructing high-quality calibration datasets.

[0053] In existing technologies, when it is necessary to correlate the location coordinates of ground sampling points with UAV multispectral imagery, the common practice is to directly locate the corresponding individual pixel in the geometrically corrected imagery based on the coordinate values, and then use the reflectance value of that pixel as the representative spectral value of the sampling point. This method assumes that the positioning is absolutely accurate and that a single pixel can completely represent the soil conditions at the location.

[0054] However, this method faces two practical challenges. First, the drone's positioning system and the georegistration process for the image inevitably contain errors, which may lead to a deviation of several pixels between the theoretical coordinates and the actual soil sampling points on the image. Second, because the soil surface is not completely homogeneous and the spatial resolution of the image is limited, the ground covered by a single pixel may contain tiny stones, vegetation remnants, or moisture differences, causing its reflectance to not purely and stably represent the soil spectrum that the sampling point intends to measure. Directly using a single pixel value will directly introduce these positioning errors and intra-pixel noise into the subsequent calibration dataset, affecting the reliability of the model calibration.

[0055] To address the aforementioned issues, this invention further integrates the coordinates of the base and additional sampling points without directly extracting individual pixel values. The ground control terminal designs an intelligent spectral value extraction process for each sample point. First, the system calculates a comprehensive initial search radius based on the nominal accuracy parameters of the UAV positioning module and the ground dimension corresponding to a single pixel of the multispectral camera. This radius takes into account potential positioning errors and adds a buffer of several pixels to define a circular neighborhood centered on the sampling point coordinates. Subsequently, the system extracts all image pixels falling within this circular neighborhood.

[0056] Next, instead of directly using these pixels, the system applies the soil pixel determination rules described earlier to filter them, retaining only those identified as high-confidence soil pixels. After filtering, the system checks the number of remaining valid soil pixels. If the number is three or more, the median of the reflectance in the near-infrared band for these pixels is calculated, and this median is used as the final image reflectance value representing that sampling point. Using the median instead of the average is to further mitigate potential interference from individual anomalous pixels.

[0057] If fewer than three valid soil pixels are found through the screening process, the system will not easily give up. It will keep the center point unchanged, proportionally expand the circular search radius, and then repeat the extraction and screening process described above. This attempt to expand the search can be performed a limited number of times to balance the probability of obtaining data and maintaining spatial representativeness. If, after several expansion searches, a sufficient number of valid soil pixels still cannot be found, the system will mark the sample point as invalid and not use it. Finally, the data of all sample points that successfully obtained representative reflectance values ​​through this process are compiled to form a pixel reflectance set for subsequent modeling.

[0058] Through this implementation method, the present invention proactively accommodates unavoidable spatial uncertainties when extracting point spectral information from images, and ensures robust data acquisition through local area statistics and multiple rounds of trials. It effectively mitigates the risks posed by inaccurate single-point positioning and mixed pixel problems, ensuring that the image reflectance values ​​ultimately used to construct the calibration dataset more realistically and reliably correspond to the soil spectra sampled in the field, thus laying a solid data foundation for obtaining a high-precision localized model.

[0059] In existing technologies, the common practice for constructing soil reflectance calibration datasets is to directly combine reflectance data from specific locations extracted from UAV imagery with reflectance data measured at the same location on the ground using a handheld spectrometer. This method assumes that the two data sources can be directly compared in terms of physical quantities, focusing primarily on the matching of location points, while neglecting to adequately consider the differences in external environmental conditions during data acquisition and the characteristics of the sensors themselves.

[0060] This approach has significant drawbacks. The acquisition of drone imagery and the measurement by ground-based handheld devices are not perfectly synchronized in time. Changes in the solar altitude angle during this period alter lighting conditions, leading to discrepancies in the reflectance measurements of the same ground object. More importantly, the multispectral camera on the drone and the ground-based handheld spectrometer are two independent sensors with different optical systems, detector responses, and radiometric calibration standards. This results in potentially non-negligible systematic deviations in the reflectance values ​​output by both, even when measuring the same stable target. Directly using these data pairs with temporal calibration differences and sensor biases for model calibration introduces fundamental errors, severely limiting the upper limit of the final inversion model's accuracy.

[0061] To address the aforementioned issues, this invention further implements a rigorous data quality control and normalization process after obtaining the pixel reflectance set and corresponding field spectral data to construct a truly usable calibration dataset. The ground control terminal first performs a spatiotemporal consistency check on each potential image-ground data pair. The system compares the time of the UAV-captured image with the time of the ground spectrometer-recorded data, retaining only data pairs with a very short time difference to ensure that the observed target has not undergone a change in state within a short period. Next, the system further acquires the solar elevation angle information recorded at these two times. Only when the solar illumination angle is also substantially consistent is the data pair considered to have experienced similar illumination conditions.

[0062] For data pairs that pass the spatiotemporal consistency test, the system does not directly use the raw image reflectance values. It has pre-established a conversion relationship between the two sensors. Before each flight mission, the operator controls the drone and handheld spectrometer to simultaneously measure the same standard diffuse reflectance reference plate. Based on these two sets of synchronized measurement data, the system fits a simple linear transformation model. When processing each pair of field data that passes the test, the system inputs the reflectance value extracted from the image into this pre-calibrated model, converting it to the reflectance reference scale of the handheld spectrometer to generate a normalized image reflectance value. Finally, this normalized value is paired with the ground-measured reference reflectance value to form a physically consistent calibration data unit that can be used for high-precision model calibration. All calibration data units generated through this process are summarized to form a high-quality soil reflectance calibration dataset.

[0063] Through this implementation method, the present invention systematically eliminates spectral measurement deviations caused by differences in observation time, variations in lighting conditions, and sensor heterogeneity in the core stage of data fusion. It ensures that every data point in the calibration dataset has high internal consistency and physical comparability, providing a clean and reliable data foundation for subsequent model parameter calibration, thereby guaranteeing the accuracy of the final soil fertility inversion model from the source.

[0064] In existing technologies, conventional optimization algorithms such as least squares are typically used when calibrating soil reflectance models for localization. These methods treat all data points used for calibration as equally important, aiming to minimize the sum of squared prediction errors for all points. The implicit assumption is that all samples in the calibration dataset are accurate and identically distributed.

[0065] However, this assumption often fails to hold true in the highly heterogeneous environment of supplementary cultivated land. Even after the aforementioned rigorous screening, the constructed calibration dataset may still contain a few data points caused by local extreme soil variations that cannot be completely eliminated, unavoidable microscale measurement errors, or residual minor registration biases. These points deviate significantly from the main data trend. If all points are treated equally during calibration, these deviation points will have an excessive impact on the optimization process, "pulling" the model parameters towards them. This will cause the final calibrated model to be biased in reflecting the main soil spectral property relationships in the region, impairing its inversion accuracy and robustness in most normal areas.

[0066] To address the aforementioned issues, this invention further employs a weighted iterative re-estimation method for model parameter optimization after obtaining a high-quality soil reflectance calibration dataset. This algorithm does not solve the problem all at once, but rather iterates through multiple loops to gradually approach the robust optimal parameters. At the beginning of each iteration, the system calculates the deviation between the model prediction value and the field reference value for each data point in the dataset, based on the model parameters used in the current iteration. Then, the system analyzes the overall distribution of the deviations across all data points and calculates a statistic that reflects the central tendency of the deviations.

[0067] Based on this measure of the dispersion of current data bias, the system dynamically calculates and assigns a specific weight value to each data point in the dataset. The core rule for weight calculation is: good points whose prediction bias is consistent with the pattern of most points will receive higher weights, while suspicious points whose bias deviates significantly from the norm will be automatically assigned lower weights. Subsequently, the system uses this set of dynamic weights to construct a weighted optimization objective function, which aims to minimize the sum of squared prediction biases of all data points. Because outliers have very low weights, their influence on the results is significantly weakened in the new objective function.

[0068] The system then invokes an optimization algorithm suitable for nonlinear models to solve for the new weighted objective, obtaining a set of updated model parameters. These new parameters are more focused on fitting the data points that are given higher weights and represent the main trend. The algorithm then uses these new parameters as a starting point for the next iteration: recalculating the deviation of all points using the new parameters, re-evaluating and reassigning weights, and solving for the updated parameters again. This process is repeated iteratively, with model parameters and data point weights adjusting and optimizing each other during iteration. The loop stops when the parameter update becomes negligible or when the preset safe iteration limit is reached. The model parameters obtained in the last iteration are the final calibration result. Substituting this result into the preset soil reflectance model structure yields the final model for localized calibration of the newly added farmland.

[0069] Through the above implementation methods, this invention introduces an intelligent, data-driven anti-interference mechanism in the model calibration stage. Instead of passively accepting all data, it actively identifies and suppresses the negative impact of potential outliers during iteration. This allows the parameter optimization process to focus on capturing stable and reliable spectral and attribute relationships within the main data set, thereby obtaining a robust model with stronger generalization ability and insensitivity to local noise. This is equivalent to putting a protective armor on the model, ensuring reliable calibration results even with a small amount of uncontrollable perturbation in the data, providing a core guarantee for the final high-precision inversion.

[0070] In existing technologies, localized and calibrated models are used to perform pixel-by-pixel calculations on UAV multispectral imagery to directly obtain predicted values ​​for soil organic matter and available nitrogen. These raw predicted data layers are then typically used to generate distribution maps. This method assumes that the model inversion results are sufficiently reliable and does not consider random noise, local model uncertainties, or extreme values ​​at the pixel scale that may be introduced during the inversion process.

[0071] This results in soil fertility distribution maps often containing significant spatial noise, manifested as isolated, anomalous pixels scattered across the image that clash with their surroundings—the so-called "salt and pepper effect." Simultaneously, the data layer may contain non-physical outliers that are spatially discontinuous and abrupt. These noises and outliers, while not representing genuine spatial variations in soil properties, severely interfere with the visual interpretation of the distribution maps, reducing their credibility and practicality for precision agriculture decision-making.

[0072] To address the aforementioned issues, this invention, after completing pixel-by-pixel model inversion of the entire region's imagery to obtain preliminary predicted data layers for soil organic matter and available nitrogen, does not immediately create maps. Instead, it executes a set of intelligent post-processing procedures designed to improve spatial data quality on each of these two data layers. First, the ground control terminal scans the entire predicted data layer, automatically identifying and marking spatially isolated anomalous pixels. The judgment logic is as follows: Centering each pixel, the system observes its immediate three-by-three pixel window. The system calculates the difference between the predicted value of the central pixel and the median of the predicted values ​​of the other eight pixels within the window. Simultaneously, it calculates the standard deviation of all predicted values ​​within this nine-pixel window. If the difference exceeds twice the standard deviation, the central pixel is identified as a spatially isolated anomalous point. For each marked anomalous pixel, the system replaces its predicted value with the median of the predicted values ​​of all unmarked pixels within its three-by-three neighborhood. This step effectively eliminates obviously abrupt noise points.

[0073] After outlier correction, the system further performs edge-aware spatial smoothing to suppress residual random fluctuations. It divides the entire prediction data layer into numerous non-overlapping square blocks. Within each block, the system performs a local spatial weighted average calculation for each cell. During weighted averaging, other cells closer to the current cell have higher weights. Through this calculation, the value of each cell moderately approaches the values ​​of other cells in its spatial neighborhood, thus smoothing out small-scale irregular fluctuations. However, to preserve the true boundaries of soil properties, the system implements a protection mechanism: at the boundaries between blocks, if the gradient of cell value change exceeds a preset threshold, it is determined that a true property edge may exist there, and cells in this area will not undergo smoothing, thus maintaining clear boundaries.

[0074] Following the two-step post-processing steps of first correcting isolated outliers and then smoothing edges, the resulting soil organic matter and available nitrogen content data layers exhibit significantly enhanced spatial continuity, visual smoothness, and reliability. Random noise and obvious errors are effectively suppressed, while the true spatial distribution pattern and boundary characteristics of the soil are preserved. The resulting spatial distribution map of soil fertility is not only more aesthetically pleasing, but more importantly, it reflects more realistic and reliable spatial variability information, providing a high-quality data foundation for precise management and fertilization decisions for supplementary arable land.

[0075] A rapid assessment device for supplementary arable land soil fertility based on UAV multispectral imaging, used to implement the assessment method of the present invention, includes: Unmanned aerial vehicle (UAV) platform; A multispectral camera, mounted on a drone platform, is used to collect multispectral images of supplementary farmland. The positioning module, located on the drone platform, is used to record the center coordinates of each multispectral image; Data radio, set up on the drone platform, is used to transmit multispectral images to the ground in real time; The ground control terminal, which communicates with the data radio, is configured as follows: Plan the flight path of the unmanned aerial vehicle (UAV) platform and control its flight; It receives multispectral images and calculates the spatial variation coefficient of reflectance in the near-infrared band for each image in real time. When the spatial variation coefficient is greater than a preset threshold, additional sample point coordinates are generated based on the high variation region in the current image. Based on the position coordinates of all basic and additional sample points, the corresponding pixel reflectance is extracted from the multispectral image to form a pixel reflectance set. The set of pixel reflectance data is combined with all the corresponding field spectral data to construct a soil reflectance calibration dataset; The parameters of the preset soil reflectance model were calibrated using the soil reflectance calibration dataset to obtain the localized calibrated model. The localized and calibrated model was used to perform pixel-by-pixel calculations on all multispectral images to retrieve the content data of soil organic matter and soil available nitrogen. A spatial distribution map of soil fertility in the supplementary cultivated land is generated based on the content data; A handheld spectrometer, operated by a ground operator, is used to collect field spectral data based on the coordinates of additional sample points and to simultaneously transmit the coordinates and data back. The mobile terminal communicates with the ground control terminal and the handheld spectrometer to receive additional sample point coordinates and guide the ground operator, and to transmit back the field spectral data and its acquisition location coordinates.

[0076] Application examples In a land consolidation project in a certain province, a new plot of farmland with an area of ​​approximately 20 hectares was added. This plot was formed by backfilling and mixing topsoil stripped from multiple surrounding areas, resulting in a complex soil origin and an expected high spatial heterogeneity in its organic matter and available nitrogen content. The project required a rapid and accurate assessment of the spatial distribution of soil fertility in this plot to provide a basis for subsequent precision fertilization and soil improvement.

[0077] Implementation process 1. System Configuration and Basic Sampling: First, a complete evaluation system was configured, including a hexacopter UAV platform, an onboard five-band multispectral camera, a high-precision GNSS positioning module, a high-speed data transmission radio, and a ground control terminal (laptop) with dedicated control software. Before the UAV flight, based on the terrain and experience, the operators pre-defined 20 basic soil sample points roughly evenly throughout the entire plot. The operator used a handheld spectrometer to measure at each point, recording the on-site spectral data and its precise GNSS coordinates.

[0078] 2. Collaborative Flight Data Acquisition: An automated flight path is planned via a ground control terminal. Considering the parameters of the onboard multispectral camera (e.g., 3.75μm pixel size, 8mm focal length) and the accuracy requirements for subsequent point extraction, the UAV's flight altitude is set at 80 meters (at this altitude, the ground sampling distance GSD is approximately 5cm, which meets the centimeter-level calculation requirements for additional sampling point coordinates). To ensure the reliability of real-time spectral variation analysis, the forward overlap is set to 75%, and the lateral overlap to 60%, providing sufficient image redundancy and a stable, effective analysis area. The UAV flies automatically along the flight path, the multispectral camera continuously acquires images, the positioning module synchronously records the center point coordinates of each image, and the data is transmitted in real-time via a data radio.

[0079] 3. Real-time mutation analysis and dynamic sampling triggering: Soil pixel purification: After receiving each image with a precise UTC timestamp and center point coordinates, the ground control terminal immediately performs pixel-level analysis. The system filters out pure soil pixels according to preset rules (red / near-infrared band reflectance ratio > 0.7, green band reflectance < 0.4, blue band reflectance < 0.3), eliminating interference from vegetation, water stains, etc.

[0080] Coefficient of variation calculation: Based solely on the selected soil pixels, calculate the mean and standard deviation of their reflectance in the near-infrared band, and obtain the spatial coefficient of variation (standard deviation / mean). The system preset threshold is 0.22.

[0081] Precisely locate outliers: When the coefficient of variation of an image exceeds 0.22, it is determined that the area needs to be sampled again.

[0082] The system divides the image into multiple rectangular sub-blocks. To ensure that each sub-block represents a locally statistically significant micro-region (e.g., covering approximately 10-15 square meters of ground), based on the current image's GSD (calculated as 5cm above), the sub-block size is set to approximately 3.2 meters on each side of the corresponding ground surface, or approximately 64 × 64 pixels (3.2 meters / 0.05 meters / pixel ≈ 64 pixels). Within each sub-block, the two pixels with the highest factor values ​​are selected as local anomalies.

[0083] Coordinate Calculation and Task Distribution: The system calculates the latitude and longitude coordinates of these anomalous pixels using image georeferenced information, generating additional sampling points. These coordinates are then transmitted to the ground operator's mobile terminal via a wireless network (such as a 4G / 5G network or a local Wi-Fi bridge) between the ground control terminal and the mobile terminal. If the network signal in the operating area is poor, the mobile terminal can download a base map of the work area in advance. The coordinates generated by the system can be attached to the base map for relative position marking, or the operator can proceed to the approximate area based on the coordinate description and then use the mobile terminal's GNSS for precise positioning to ensure the sampling task can be executed.

[0084] 4. On-site supplementary sampling: The ground operator, following the navigation guidance on the tablet, goes to the newly assigned coordinate point, uses a handheld spectrometer to conduct on-site measurements, and immediately transmits the spectral data and confirmed coordinates back.

[0085] 5. Robust Spectral Data Extraction and Registration: After all flights and sampling are completed, the ground control terminal integrates the coordinates of all 25 valid sample points (20 base points + 5 dynamic additional points). Each point is verified according to the aforementioned validity criteria: For each point, the system calculates the initial buffer radius based on its coordinates, positioning error, and pixel size, and extracts all pixels within this buffer from the image. The soil pixel criteria are then applied again for filtering. If the number of pure soil pixels selected is no less than 3, the median is taken as the representative image reflectance value for that point, ultimately forming a pixel reflectance set.

[0086] 6. Construction of a High-Consistency Calibration Dataset: The system constructs a high-consistency calibration dataset. For each pixel reflectance value in the image reflectance set and its paired field spectral data, the system performs rigorous screening and correction according to the following steps to form a valid calibration data pair: Spatiotemporal consistency check: The system calculates the absolute difference between the timestamps of image acquisition and field measurement, retaining only data pairs with a time difference ≤ 5 minutes. For these data pairs, the system further calculates the solar altitude angle at the acquisition time by calling the built-in astronomical algorithm or querying the ephemeris, based on their respective UTC timestamps and geographic coordinates (the center point coordinates for images and the sampling point coordinates for field measurements), and retains only data pairs with a solar altitude angle difference ≤ 5 degrees to ensure that the illumination conditions are basically consistent.

[0087] Sensor spectral normalization: For data pairs that pass the spatiotemporal consistency test, the system uses a linear normalization model calibrated before flight by synchronous measurement on a standard diffuse reflectance reference plate to correct the reflectance values ​​extracted from the image. This eliminates the system spectral response differences between the UAV multispectral camera and the handheld spectrometer, ensuring that the corrected image reflectance values ​​are at a unified physical benchmark with the reference reflectance values ​​measured in the field.

[0088] Calibration dataset generation: The corrected image reflectance values ​​are paired with their corresponding field reference reflectance values ​​to form a calibration data unit. All calibration data units generated through this process are aggregated to form a highly consistent soil reflectance calibration dataset for model localization calibration. In this application example, a calibration dataset containing 25 such highly consistent data pairs was successfully constructed.

[0089] 7. Robust Model Localization Calibration: A weighted iterative re-estimation method is employed to optimize the parameters of the pre-defined soil organic matter-spectrum and alkaline nitrogen-spectrum inversion models using the aforementioned calibration dataset. During iteration, this algorithm automatically assigns low weights to potentially outlier data points that deviate from the mainstream trend (e.g., using a robust weighting function to give smaller weights to data points with exceptionally large prediction residuals), focusing parameter optimization on the main data and ultimately obtaining a localized calibration model for this site.

[0090] 8. Global Inversion and Intelligent Post-processing: A localized model was applied to perform pixel-by-pixel calculations on all UAV imagery to obtain preliminary distribution maps of soil organic matter and available nitrogen content. Two-step post-processing was then performed: First, spatially isolated outlier pixel values ​​were identified and replaced with the median of a 3x3 neighborhood; second, spatial smoothing with edge preservation was performed to suppress small-scale noise while retaining the true boundaries of abrupt content changes.

[0091] 9. Results Generation: The processed data is automatically rendered by the ground control terminal to generate clear and smooth spatial distribution maps of soil organic matter content and soil available nitrogen content, which are the soil fertility assessment results of the supplemented cultivated land.

[0092] To objectively evaluate the effectiveness of this method, a traditional fixed-grid sampling assessment method was implemented in parallel in the same supplementary cultivated land area as a control. The traditional method pre-determines 25 sampling points based on a regular grid before flight, resulting in a uniform spatial distribution, but this is entirely based on prior experience and does not consider the actual spectral information available during flight. The results show that the pre-determined grid points failed to effectively capture two existing high-variability soil fertility patches, leading to a spatial representativeness blind spot in the calibration dataset. The soil property inversion model calibrated based on this dataset performed poorly at independent validation points, with average relative errors of 18.5% and 22.1% for soil organic matter and available nitrogen, respectively. Furthermore, the generated distribution map exhibited significant noise and discontinuous anomalous patches, affecting the intuitiveness and reliability of decision-making.

[0093] In contrast, the evaluation method of this invention, using only 20 pre-set basic sampling points, dynamically supplemented 5 sampling points through real-time spectral variation analysis during UAV flight. These 5 points were automatically located by the system in localized areas with significant spectral variation, with most falling within highly heterogeneous patch regions missed by traditional methods, significantly enhancing the targeted coverage of key variation areas by the sampling. The resulting calibration dataset exhibited significantly enhanced spatial representativeness, leading to a substantial improvement in model calibration accuracy, with determination coefficients for soil organic matter and available nitrogen inversion models reaching 0.82 and 0.79, respectively. Validation at the same independent validation points showed that the average relative errors of the inversion values ​​decreased to 10.2% and 12.7%, respectively, demonstrating a significant improvement in accuracy.

[0094] The resulting spatial distribution map of soil fertility, after intelligent post-processing, exhibits good spatial continuity and effectively suppresses image noise (reducing the number of spatially isolated outliers by approximately 85% and the standard deviation of image noise by approximately 60%), clearly presenting the spatial gradient and boundaries of true soil properties. In summary, this invention, through intelligent sampling and end-to-end refined processing in a coordinated space-ground approach, achieves significantly enhanced spatial representativeness of sampling points with a comparable total number of sampling points (25 in total). This is because dynamic sampling points accurately cover highly variable areas (80% of dynamic points fall within the validated highly heterogeneous patches), resulting in higher assessment accuracy and reliability than traditional fixed-grid sampling methods (model determination coefficient R² increased to 0.82 and 0.79, and average relative error reduced to 10.2% and 12.7%). Simultaneously, this method achieves efficient and precise allocation of manpower and time resources for field sampling (reducing total operation time by approximately 40%, labor costs by approximately 25%, and overall project costs by approximately 15%), providing a practical and advanced technical means for the rapid and accurate assessment of farmland soil fertility.

[0095] The dynamic adaptive sampling and refined processing flow of this invention significantly improves the sampling representativeness, model calibration accuracy, inversion accuracy, and result map quality of soil fertility assessment in highly heterogeneous environments such as supplementary cultivated land, achieving a more reliable and rapid assessment with higher overall efficiency.

[0096] Although the embodiments of the present invention have been disclosed above, they are not limited to the applications listed in the specification and embodiments. They can be applied to various fields suitable for the present invention. For those skilled in the art, other modifications can be easily made. Therefore, without departing from the general concept defined by the claims and their equivalents, the present invention is not limited to the specific details.

Claims

1. A rapid assessment method for supplementary arable land soil fertility based on UAV multispectral imaging, characterized in that, Includes the following steps: The configuration includes a drone platform, a multispectral camera mounted on the platform, a positioning module, a data transmission radio, and a ground control terminal. The flight path is planned through the ground control terminal, and the UAV platform is controlled to fly and collect multispectral images of the farmland. At the same time, the positioning module records the center position coordinates of each image. Multispectral images are transmitted in real time to the ground control terminal via a data transmission radio. The ground control terminal calculates the spatial variation coefficient of reflectance of each image in the near-infrared band in real time. The spatial variation coefficient is the ratio of the standard deviation to the mean. When the spatial variation coefficient is greater than a preset threshold, the ground control terminal generates additional sample point coordinates based on the high variation area in the current image and sends them to the ground operator's mobile terminal. The preset threshold is 0.15-0.30; Ground operators use handheld spectrometers to collect on-site spectral data based on coordinates and simultaneously transmit the coordinates and data back. Before the drone platform takes off, 10-30 basic soil sample points are set up in the supplementary cultivated land, and the field spectral data and location coordinates of all basic sample points are collected. Based on the position coordinates of all basic and additional sample points, the ground control terminal extracts the corresponding pixel reflectance from the multispectral image to form a pixel reflectance set. The set of pixel reflectance data is combined with all the corresponding field spectral data to construct a soil reflectance calibration dataset; The parameters of the preset soil reflectance model were calibrated using the soil reflectance calibration dataset to obtain the localized calibrated model. The localized and calibrated model was used to perform pixel-by-pixel calculations on all multispectral images to retrieve the content data of soil organic matter and soil available nitrogen. Based on soil organic matter and available nitrogen content data, a spatial distribution map of soil fertility for supplementary cultivated land is generated.

2. The method for rapid assessment of supplementary arable land soil fertility based on UAV multispectral imaging according to claim 1, characterized in that, In the step of planning the flight route through the ground control terminal, based on the preset accuracy of generating additional soil sample point coordinates and the spatial resolution of the multispectral camera, the flight altitude of the UAV platform is set to 50-150 meters. According to the requirements of subsequent real-time calculation of spatial variation coefficient, the forward overlap is set to 65%-80% and the lateral overlap is set to 55%-70%, so as to ensure that the acquired multispectral images can not only meet the geometric positioning accuracy required to extract additional sample point coordinates from the images, but also provide image data with sufficient redundancy for real-time calculation of the spatial variation coefficient of near-infrared reflectance of each image.

3. The method for rapid assessment of supplementary arable land soil fertility based on UAV multispectral imaging according to claim 1, characterized in that, The specific steps involved in the ground control terminal calculating the spatial variation coefficient of reflectance in the near-infrared band for each image in real time include: After receiving the current multispectral image, soil pixel determination is performed on each pixel in the image. The determination criteria are: the ratio of red band reflectance to near-infrared band reflectance is greater than 0.7, green band reflectance is less than 0.4, and blue band reflectance is less than 0.

3. Pixels that simultaneously meet the above criteria are recorded as soil pixels. Based on all soil pixels, their reflectance values ​​in the near-infrared band are extracted, and the mean and standard deviation of these values ​​are calculated. The ratio obtained by dividing the standard deviation by the mean is the spatial variation coefficient of the image used to trigger additional sampling. The subsequent step of generating additional sample point coordinates is performed only if the spatial coefficient of variation is greater than a preset threshold.

4. The method for rapid assessment of supplementary arable land soil fertility based on UAV multispectral imaging according to claim 1, characterized in that, The steps for the ground control terminal to generate additional sample point coordinates based on highly variable areas in the current imagery specifically include: The ground control terminal performs image block processing on the current multispectral image with a spatial variation coefficient greater than a preset threshold, dividing it into several rectangular sub-regions; For each rectangular sub-region, the local anomaly factor value of the reflectance of all pixels within it in the near-infrared band is calculated based on the k-nearest neighbor distance, which is used to characterize the degree of deviation of the pixel reflectance from its local neighborhood. Within each rectangular sub-region, the local anomaly factor values ​​are sorted in descending order, and the top N pixels are selected, where N is an integer between 1 and 5. The row and column coordinates of each selected pixel are combined with the center coordinates of the image, the spatial resolution of the multispectral camera, and the imaging geometry model to calculate the corresponding geographic coordinates. The calculated geographic coordinates are used as the coordinates of additional sample points.

5. The rapid assessment method for supplementary arable land soil fertility based on UAV multispectral imaging according to claim 4, characterized in that, In the step of calculating the local anomaly factor value of the reflectance in the near-infrared band for all pixels within each rectangular sub-region, the k value on which the k-nearest neighbor distance used to calculate the local anomaly factor value depends is dynamically determined as follows: The ground control terminal first calculates the average and standard deviation of the near-infrared reflectance of all pixels within the rectangular sub-region; Based on the mean and standard deviation, a first reflectance threshold and a second reflectance threshold are set. The first reflectance threshold is the mean plus one standard deviation, and the second reflectance threshold is the mean minus one standard deviation. The number of pixels within the rectangular sub-region whose reflectance values ​​fall between the first reflectance threshold and the second reflectance threshold is counted and recorded as the core pixel count of the sub-region. The number of core pixels in the rectangular sub-region is compared with the preset minimum nearest neighbor number K. min and the number of maximum nearest neighbors K max Compare; K min For 5, K max It is an integer corresponding to 20% of the total number of cells in the rectangular sub-region, and does not exceed 50; If the number of core pixels is ≤ K min Then the value of k is determined to be K. min If the number of core pixels is ≥ K max Then the value of k is determined to be K. max If the number of core pixels is between K min With K max Between these values, the k value is determined to be equal to the number of core pixels; Based on the determined k value, the local anomaly factor value of each cell in the sub-region is calculated.

6. The method for rapid assessment of supplementary arable land soil fertility based on UAV multispectral imaging according to claim 3, characterized in that, The ground control terminal extracts the corresponding pixel reflectance from the multispectral image based on the position coordinates of all basic and supplementary sample points, forming a pixel reflectance set. This process specifically includes: For the location coordinates of each sample point, the ground control terminal performs the following operations: Based on the positioning accuracy parameters of the positioning module and the spatial resolution of the multispectral camera, the expected position uncertainty radius of the sample point is calculated. The expected position uncertainty radius is equal to the sum of the positioning accuracy error value and the size of the three pixels. Using the coordinates of the sample point as the center and the radius of the uncertain expected location as the radius, a circular neighborhood is defined on the corresponding multispectral image, and all pixels within the neighborhood are extracted. Based on the soil pixel identification criteria, pixels belonging to soil pixels are selected from the extracted pixels; If the number of selected soil pixels is ≥3, the median of the reflectance of these soil pixels in the near-infrared band is calculated as the pixel reflectance value representing the sample point. If the number of soil pixels selected is less than 3, the center coordinates of the circular neighborhood remain unchanged, the radius is expanded to 1.5 times its original value, the circular neighborhood is redefined and the selection is repeated. The radius expansion step is performed a maximum of 3 times. If the number of soil pixels selected after expanding the radius three times is still less than three, then the sample point is marked as an invalid sample point. Collect the pixel reflectance values ​​corresponding to all sample points that were not marked as invalid sample points to form a pixel reflectance set.

7. The method for rapid assessment of supplementary arable land soil fertility based on UAV multispectral imaging according to claim 1, characterized in that, The steps for constructing a soil reflectance calibration dataset by combining the pixel reflectance set with all corresponding field spectral data include: For each pixel reflectance value in the pixel reflectance set and its paired field spectral data, the ground control terminal performs the following operations to construct a valid data pair: Obtain the image acquisition timestamp and the field spectral acquisition timestamp. If the absolute difference between the two is less than or equal to 5 minutes, then mark them as valid time-synchronized data pairs. Acquire solar elevation angle data at the time of image acquisition and the time of on-site spectral acquisition. If the difference between the two is less than or equal to 5 degrees, then the time-synchronized valid data pair is further marked as a lighting condition matching data pair. For each pair of lighting condition matching data, the average reflectance value of the field spectral data in the corresponding spectral range of the near-infrared band of the multispectral camera is extracted as the field reference reflectance, and the corresponding pixel reflectance value is read as the image extracted reflectance. The field reference reflectance and the image extracted reflectance are input into a preset linear normalization model to correct the image extracted reflectance and obtain the normalized image reflectance value. The parameters of the linear normalization model are determined by controlling the UAV platform and the handheld spectrometer to perform synchronous spectral measurements and fitting on the same standard diffuse reflectance reference plate before each flight mission. The normalized image reflectance value is paired with its corresponding field reference reflectance value to form a calibration binary; All calibration pairs are aggregated to form a soil reflectance calibration dataset.

8. The method for rapid assessment of supplementary arable land soil fertility based on UAV multispectral imaging according to claim 1, characterized in that, The parameters of the pre-set soil reflectance model are calibrated using a soil reflectance calibration dataset to obtain a locally calibrated model, specifically including: The ground control terminal uses all normalized image reflectance values ​​in the soil reflectance calibration dataset as model input data and all corresponding field reference reflectance values ​​as expected model output data. The preset soil reflectance model contains a set of model parameters to be determined. The weighted iterative re-estimation method is used to optimize and calibrate the determined model parameters. The weighted iterative re-estimation method performs the following steps through multiple iterations: a) Based on the model parameters of the current iteration, calculate the deviation between the model-predicted reflectance value and its field reference reflectance value for each data point in the calibration dataset; b) Calculate the median of the absolute values ​​of all deviations as a measure of the dispersion of the current data deviation; c) Calculate a weight value for each data point based on the deviation dispersion scale, so that the data point with the greater deviation from the mainstream trend is assigned a lower weight; d) Construct a weighted optimization objective using the calculated weights. The optimization objective is to minimize the sum of squared weighted prediction biases for all data points. e) Solve for the weighted optimization objective and update the model parameters; f) Use the updated model parameters as the starting point for the next iteration, and repeat steps a) to e). Set an iteration termination condition: stop iterating when the change in model parameters between two consecutive iterations is less than a preset threshold or the total number of iterations reaches a preset upper limit. The model parameters obtained from the last iteration are used as the calibration results and substituted into the preset soil reflectance model to obtain the locally calibrated model.

9. The method for rapid assessment of supplementary arable land soil fertility based on UAV multispectral imaging according to claim 1, characterized in that, The localized and calibrated model was applied to perform pixel-by-pixel calculations on all multispectral images to retrieve soil organic matter and available nitrogen content data, specifically including: The ground control terminal inputs the reflectance values ​​of each pixel in each band of all multispectral images into the locally calibrated model pixel by pixel. The model outputs the corresponding predicted values ​​of soil organic matter content and soil available nitrogen content. For each image inversion layer of predicted data, the following post-processing steps are performed: First, identify and label spatially isolated outlier cells: Take a 3×3 cell window centered on each cell, calculate the absolute difference between the predicted value of the center cell and the median of the predicted values ​​of the other 8 cells in the window. If the difference exceeds twice the standard deviation of the predicted values ​​of all cells in the window, it is labeled as an outlier cell. Replace the value of each labeled outlier cell with the median of the predicted values ​​of all unlabeled cells in its 3×3 neighborhood. Secondly, edge-aware spatial smoothing is performed on the prediction data layer. The data layer is divided into several non-overlapping 32×32 pixel blocks. For each pixel in each block, the spatial weighted average of its predicted value and the predicted values ​​of other pixels in the block is calculated. The weights are calculated according to the inverse proportional function of the spatial Euclidean distance between pixels, with higher weights for closer pixels. The original predicted value of the pixel is replaced with this spatial weighted average, but areas where the gradient change at the edge between blocks is greater than a preset threshold are not smoothed. The predicted data layer after the above post-processing is used as the final inversion data for soil organic matter content and soil available nitrogen content.

10. A rapid assessment device for supplementary arable land soil fertility based on UAV multispectral imaging, used to implement the method as described in any one of claims 1-9, characterized in that, include: Unmanned aerial vehicle (UAV) platform; A multispectral camera, mounted on a drone platform, is used to collect multispectral images of supplementary farmland. The positioning module, located on the drone platform, is used to record the center coordinates of each multispectral image; Data radio, set up on the drone platform, is used to transmit multispectral images to the ground in real time; The ground control terminal, which communicates with the data radio, is configured as follows: Plan the flight path of the unmanned aerial vehicle (UAV) platform and control its flight; It receives multispectral images and calculates the spatial variation coefficient of reflectance in the near-infrared band for each image in real time. When the spatial variation coefficient is greater than a preset threshold, additional sample point coordinates are generated based on the high variation region in the current image. Based on the position coordinates of all basic and additional sample points, the corresponding pixel reflectance is extracted from the multispectral image to form a pixel reflectance set. The set of pixel reflectance data is combined with all the corresponding field spectral data to construct a soil reflectance calibration dataset; The parameters of the preset soil reflectance model were calibrated using the soil reflectance calibration dataset to obtain the localized calibrated model. The localized and calibrated model was used to perform pixel-by-pixel calculations on all multispectral images to retrieve the content data of soil organic matter and soil available nitrogen. A spatial distribution map of soil fertility in the supplementary cultivated land is generated based on the content data; A handheld spectrometer, operated by a ground operator, is used to collect field spectral data based on the coordinates of additional sample points and to simultaneously transmit the coordinates and data back. The mobile terminal communicates with the ground control terminal and the handheld spectrometer to receive additional sample point coordinates and guide the ground operator, and to transmit back the field spectral data and its acquisition location coordinates.