A large-scale farmland biomass monitoring method based on unmanned aerial vehicle and satellite image

Abstract

The application discloses a large-scale farmland biomass monitoring method based on unmanned aerial vehicle and satellite image and belongs to the technical field of agricultural supervision.The application aims to solve the problem that the precision of an existing biomass estimation model based on satellite remote sensing data is limited by the mismatch between satellite data and ground measurement in spatial scale, and provides a biomass estimation precision modeling method based on satellite remote sensing data by taking unmanned aerial vehicle data as an intermediate bridge and combining with a proposed SS-GMI sampling method.The method is a general method and is suitable for modeling of crop biomass estimation based on satellite data; moreover, by combining with the proposed SS-GMI sampling method suitable for farmland geographic space, more representative and independent sample point sets can be obtained while unmanned aerial vehicle and satellite image data are fully utilized, a more reliable basis is provided for analysis and modeling of geographic space data, and the precision and reliability of crop biomass estimation are further improved.

Description

Technical Field

[0001] This invention belongs to the field of agricultural monitoring technology, specifically relating to a large-scale farmland biomass monitoring method based on drones and satellite imagery. Background Technology

[0002] With population growth and increasing extreme weather events, food security faces challenges. Crop biomass is one of the key factors determining yield. Accurately estimating crop biomass can help farmers and researchers better understand crop growth, yield potential, and environmental impact. Accurate crop biomass estimation can guide growers to optimize planting density, fertilization programs, and harvest timing to maximize crop yield and quality. As a crucial factor in yield forecasting, accurate crop biomass estimation can also help farms and processing plants predict raw material supply and processing schedules.

[0003] Currently, remote sensing has become the primary data source for estimating aboveground biomass over large areas. Unmanned aerial vehicle (UAV) remote sensing platforms, characterized by high temporal and spatial resolution, can carry multiple sensors to obtain high-resolution multimodal information at the canopy scale. Combined with ground surveys, they have been widely applied in aboveground biomass estimation in agriculture and forestry. Compared to UAV remote sensing, satellite remote sensing can more efficiently measure vegetation cover over large areas, playing a crucial role in large-scale aboveground biomass estimation. However, limited by resolution and revisit cycles, satellite remote sensing struggles to accurately estimate crop biomass over the entire crop growth cycle at the farmland scale. Currently, with the continuous updating and improvement of satellite sensor technology, the spatial resolution of satellite remote sensing data has significantly improved, as seen in satellites such as the GeoEye series, Quickbird, Ikonos, Planet series, and the Chinese Gaofen series. However, because aboveground biomass surveys primarily involve labor-intensive and time-consuming destructive sampling, the sampling units for most available farmland-level sampling point data are relatively small. Therefore, matching high-resolution satellite remote sensing with ground surveys still presents challenges. Discrepancies in matching ground reference data with satellite data can significantly impact the accuracy of satellite remote sensing-based models for estimating aboveground biomass. Therefore, a method is needed to accurately match ground sampling and satellite imagery to facilitate the construction of farmland-level biomass estimation models based on satellite data. Unmanned aerial vehicles (UAVs), with their high-resolution mapping capabilities, have the potential to serve as an intermediary bridge in matching ground sampling and satellite data.

[0004] Currently, machine learning models are commonly used to estimate aboveground biomass based on remote sensing methods. These models achieve high-precision prediction capabilities through training on large-scale data. However, due to their "black box" nature, the internal workings of these models are difficult to interpret, making it challenging to understand the specific prediction process in some cases. Furthermore, the robustness and scalability of machine learning models typically rely on extensive ground-based measurement data for training and validation. Incomplete or biased ground-based observation data can negatively impact the accuracy and applicability of machine learning models. In contrast, allometric growth models offer advantages because they are based on the biophysical processes of vegetation growth, simulating energy acquisition and utilization during growth. These models typically estimate vegetation biomass based on other crop parameters, making them easier to interpret and understand. Currently, research on estimating crop biomass using remote sensing data combined with allometric growth models is relatively limited, and studies utilizing more features from remote sensing data, such as vegetation spectral characteristics, for biomass allometric growth models are also insufficient. Furthermore, while random sampling is a widely used traditional sampling method, its application in farmland geospatial data suffers from issues such as data concentration, uneven coverage, and potential spatial autocorrelation between samples. In contrast, systematic sampling is limited by the choice of step size, may also exhibit spatial autocorrelation between samples, and is less efficient. Therefore, in the analysis and modeling of farmland geospatial data, a new and suitable sampling method is needed to fully utilize UAV and satellite data from farmland, providing a reliable foundation for geospatial data analysis and modeling. Summary of the Invention

[0005] The purpose of this invention is to address the limitation of existing biomass estimation models based on satellite remote sensing data in terms of accuracy, which is caused by the mismatch between the spatial scale of satellite data and ground measurements. This invention provides a method for accurate biomass estimation modeling based on satellite remote sensing data, which uses UAV data as an intermediate bridge and combines it with the proposed SS-GMI sampling method.

[0006] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0007] A method for large-scale farmland biomass monitoring based on UAVs and satellite imagery, the method being:

[0008] Step 1: Use drones to collect digital and multispectral images of multiple crop test areas at multiple time points, and perform preprocessing;

[0009] Step 2: Collect satellite imagery data covering the areas surveyed by drones multiple times in Step 1, and perform preprocessing;

[0010] Step 3: Measure the aboveground dry weight and plant height of all sugar beet varieties in the experimental areas;

[0011] Step 4: Calculate the growth degree days and relative growth degree days for each experimental area based on meteorological data;

[0012] Step 5: Based on the UAV preprocessing results, determine the common phenological zone of sugar beets in all test areas and estimate the sugar beet plant height, and verify the estimation accuracy using actual measured plant height;

[0013] Step Six: After selecting the better variables and models using one of the test areas, combine phenological periods and heat index, and evaluate the performance of the model when applied to different test areas;

[0014] Step 7: Biomass estimation by matching field sampling data with drone and satellite data:

[0015] UAV data was collected, and a biomass estimation model was constructed based on allometric growth relationships. Then, aboveground biomass mapping was performed based on the UAV data. To match the scale of satellite imagery, the biomass map based on UAV data needed to be downsampled to match the 3m and 10m resolutions of Planet Scope and Sentinel-2A. The downsampled UAV biomass map was then sampled (see step eight below for details) to construct a sample set for building the biomass estimation model based on satellite data, recording the latitude and longitude coordinates and biomass values ​​of the sampling points. Since the current satellite data resolution cannot accurately estimate crop height within 1 meter, satellite spectral data was used as a variable to construct the biomass estimation model based on satellite imagery. The reflectance values ​​of each satellite band corresponding to the sample point location and the aboveground biomass values ​​of the sampling points were used as the training set to construct the model.

[0016] Step 8: The SS-GMI sampling method, which combines systematic sampling, global Moran index and optimal sampling strategy, is as follows: Set the initial step size s, and start sampling points from the starting point (0, 0) of the image with image size (m, n) to form a point set X;

[0017]

[0018] Where s is the initial step size, i and j are the number of points horizontally and vertically in an image of size (m, n), respectively; the number N of the point set X is:

[0019]

[0020] After obtaining the point set X, calculate its global Moran index:

[0021]

[0022] Where N is the number of points in the point set, and x represents the value of a point. W represents the average value of the points. ij It is a spatial weight matrix;

[0023] When calculating the weight matrix W ij (Eq.7) uses the K-Nearest Neighbor (KNN) weighting method, and the expression for the weight matrix is ​​as follows:

[0024]

[0025] The advantages of this invention over existing technologies are as follows: This method uses UAV data as an intermediate bridge to solve the matching error problem that exists when directly matching small-scale ground measurements with coarse-resolution satellite data. This method is a general approach suitable for modeling crop biomass estimation based on satellite data. Furthermore, by combining it with the proposed SS-GMI sampling method applicable to farmland geospatial data, it can obtain a more representative and independent sample point set while fully utilizing both UAV and satellite imagery data, providing a more reliable foundation for geospatial data analysis and modeling, and further improving the accuracy and reliability of crop biomass estimation. Attached Figure Description

[0026] Figure 1 Layout diagram for Experiment 1;

[0027] Figure 2 Layout diagram for Experiment 2;

[0028] Figure 3 Layout diagram for Experiment 3;

[0029] Figure 4 A diagram illustrating foreground and background segmentation methods for UAV imagery;

[0030] Figure 5 A power function graph showing the relationship between the aboveground dry weight of sugar beets and plant height.

[0031] Figure 6 A power function graph showing the relationship between the aboveground dry weight of sugar beets and the Dreg exponent;

[0032] Figure 7 A graph showing the power function relationship between the aboveground dry weight of sugar beets and plant height * Dreg;

[0033] Figure 8 The graph shows the verification of the power function model of plant height and aboveground dry weight in Experiment 2;

[0034] Figure 9 The graph shows the verification of the power function model of plant height and aboveground dry weight in Experiment 3;

[0035] Figure 10 The graph shows the validation of the power function model of Dreg index and aboveground dry weight in Experiment 3;

[0036] Figure 11 The graph shows the validation of the power function model of plant height * Dreg versus aboveground dry weight in Experiment 3.

[0037] Figure 12 This is a validation plot of the partial least squares regression model in Experiment 2 when plant height is used as a variable;

[0038] Figure 13 This is a validation plot of the random forest regression model in Experiment 2 when plant height is used as a variable;

[0039] Figure 14 This is a validation plot of the support vector machine regression model in Experiment 2 when plant height is used as a variable;

[0040] Figure 15 This is a validation plot of the partial least squares regression model in Experiment 3 when plant height is used as a variable;

[0041] Figure 16 This is a validation plot of the random forest regression model in Experiment 3 when plant height is used as a variable;

[0042] Figure 17 This is a validation plot of the support vector machine regression model in Experiment 3 when plant height is used as a variable;

[0043] Figure 18 The partial least squares regression model was validated in Experiment 3 when Dreg was used as a variable;

[0044] Figure 19 The validation plot of the random forest regression model on Experiment 3 when Dreg is used as a variable;

[0045] Figure 20 This is a validation plot of the support vector machine regression model in Experiment 3 when Dreg is used as a variable;

[0046] Figure 21 This is a validation plot of the partial least squares regression model in Experiment 3 when plant height * Dreg is used as a variable;

[0047] Figure 22 This is a validation plot of the random forest regression model in Experiment 3 when plant height * Dreg is used as a variable;

[0048] Figure 23 This is a validation plot of the support vector machine regression model on Experiment 3 when plant height *Dreg is used as a variable;

[0049] Figure 24 The verification graph of the power function model that combines phenology and heat index with plant height as a variable in Experiment 2;

[0050] Figure 25The verification graph of the power function model that combines phenology and heat index with plant height as a variable in Experiment 3;

[0051] Figure 26 The validation graph of the power function model that combines phenology and heat index with plant height * Dreg as the variable in Experiment 3;

[0052] Figure 27 A map of aboveground biomass of sugar beets as of July 20, 2023, based on power-law relationships of plant height, using UAV data.

[0053] Figure 28 A map of aboveground biomass of sugar beets, drawn on August 7, 2023, using UAV data, based on a power function relationship of plant height.

[0054] Figure 29 A map of aboveground biomass of sugar beets as of August 20, 2023, based on power-law relationships of plant height, using UAV data.

[0055] Figure 30 This is a spatial distribution map of points sampled from the UAV ground biomass map on July 20, 2023, based on the SS-GMI sampling method.

[0056] Figure 31 This is a frequency distribution map of samples taken from the UAV ground biomass map on July 20, 2023, based on the SS-GMI sampling method.

[0057] Figure 32 This is a spatial distribution map of points sampled from the UAV ground biomass map on August 7, 2023, based on the SS-GMI sampling method.

[0058] Figure 33 This is a frequency distribution map of samples taken from the UAV ground biomass map on August 7, 2023, based on the SS-GMI sampling method.

[0059] Figure 34 This is a spatial distribution map of points sampled from the UAV ground biomass map on August 20, 2023, based on the SS-GMI sampling method.

[0060] Figure 35 This is a frequency distribution map of samples taken from the UAV ground biomass map on August 20, 2023, based on the SS-GMI sampling method.

[0061] Figure 36 The graph shows the Moran's I, Geary's C, and Getis's G indices and their significance test for a sample set from ten samplings on July 20, 2023, based on random sampling.

[0062] Figure 37The graph shows the Moran's I, Geary's C, and Getis's G indices and their significance test for a sample set from ten samplings on August 7, 2023, based on random sampling.

[0063] Figure 38 The graph shows the Moran's I, Geary's C, and Getis's G indices and their significance test for a sample set from ten samplings on August 20, 2023, based on random sampling.

[0064] Figure 39 A power function graph showing the relationship between the red band of PS imagery and aboveground biomass.

[0065] Figure 40 The graph shows the power function relationship between the red band of PS imagery and aboveground biomass on the validation dataset.

[0066] Figure 41 Biomass mapping based on PS imagery on 2023-07-21;

[0067] Figure 42 Biomass mapping based on PS imagery on 2023-08-08;

[0068] Figure 43 Biomass mapping based on PS imagery on August 21, 2023;

[0069] Figure 44 A diagram of a biomass estimation model based on S2A imagery;

[0070] Figure 45 Spatial mapping of biomass in the sugar beet growing area of ​​Suqin Farm based on S2A imagery created on July 16, 2023;

[0071] Figure 46 Spatial mapping of biomass in the sugar beet growing area of ​​Suqin Farm based on S2A imagery drawn on August 5, 2023;

[0072] Figure 47 Spatial mapping of biomass in the sugar beet growing area of ​​Suqin Farm based on S2A imagery drawn on August 21, 2023. Detailed Implementation

[0073] The technical solution of the present invention will be further described below with reference to the accompanying drawings and embodiments, but it is not limited thereto. Any modifications or equivalent substitutions to the technical solution of the present invention that do not depart from the spirit and scope of the technical solution of the present invention should be covered within the protection scope of the present invention.

[0074] Example 1:

[0075] (1) Experimental Design

[0076] This model was validated in the field using sugar beets as a representative example. Experiment 1 (Exp. 1) was conducted in 2021 at the Liangcheng Experimental Station in Inner Mongolia (40.5021N, 112.1452E). A total of 132 sugar beet varieties were selected for breeding trials, including 2 from the United States, 18 from the United Kingdom, 30 from Germany, and 82 from the Netherlands. A randomized block design was used, with one variety planted in each plot and one replicate for each variety in different plots. Sowing was carried out on May 18, 2021. Each plot was 1.2m × 2m in size, with a row spacing of 0.25m and a density of 15 plants / m². 2 Field management methods were consistent with local conventional agricultural production practices. Ground control points (GCPs) were evenly distributed in the experimental area before sowing, and their locations were measured for geolocation correction and registration during subsequent image processing. The experimental plots were arranged as follows: Figure 1 As shown.

[0077] Experiment 2 was conducted in 2021 at Suqin Farm in Sanhe Hui Township, Hulunbuir City, Inner Mongolia Autonomous Region (50.4823N, 119.8837E). The experimental area was 60m × 50m, and 10 high-quality commercial sugar beet varieties suitable for growth in high-latitude northern regions were used as experimental materials: KN1387 (KN), KUHN1357 (KUHN), KW2314 (KW), LN17101-1210 (LN1), LN17101 (LN2), LN17101-1805 (LN3), MK4044 (MK1), MK4085 (MK2), SR411 (SR), and SV1588 (SV). The experimental area consisted of 30 plots, with one variety planted in each plot. Each variety had three replicate plots, and each plot measured 1.5m × 15m. Three rows of sugar beets were planted in each plot, with a row spacing of 50cm and a plant spacing of 20cm, resulting in a planting density of 10 plants / m². 2 Field management methods were carried out according to local conventions. Ground control points (GCPs) were evenly distributed in the experimental area before sowing, and their locations were measured for geolocation correction and registration during later image processing. The experimental plots were arranged as follows: Figure 2 As shown.

[0078] Experiment 3 was conducted in a field production area in Sanhe Hui Township, Hulunbuir City, Inner Mongolia Autonomous Region (50.462N, 119.855E). The sugar beet variety planted in this field was LN1701 (LN2), sown on May 9, 2023, with a row spacing of 50cm, a plant spacing of 19cm, and a planting density of 10.5 plants / m². 2 Fertilization, irrigation, and field management were all carried out according to local conventional methods. The experimental setup was as follows: Figure 3 As shown.

[0079] Experiment 1 is located in a mid-temperate continental climate zone, with an average annual temperature of 2-5℃, a frost-free period of 80-125 days, an average annual precipitation of 350-450 mm, an annual evaporation of 1938 mm, and an average sunshine duration of 3026.7 hours. Experiments 2 and 3 are located in a mid-temperate continental grassland climate zone, with an average annual temperature of -5 to 3℃, an average annual precipitation of 358.6 mm, sunshine duration of 2500-3000 hours, and a frost-free period of approximately 89 days.

[0080] (2) Data Acquisition

[0081] (2.1) UAV data acquisition and preprocessing

[0082] Digital and multispectral images from drones were collected at multiple time points (see Table 1). The RGB and multispectral images of the Experiment 1 area were captured by a DJI Hannover 4 RTK drone (equipped with a 1-inch CMOS sensor with 20 million effective pixels) and a Parrot Sequoia multispectral camera carried by a Parrot Bluegrass drone (Parrot, Paris, France), respectively. The RGB images for Experiment 2 were acquired using a DJI Hannover 4 RTK drone. The RGB and multispectral images of the Experiment 3 area were acquired using a DJI Hannover 4 Multispectral (P4M) drone (DJI, Shenzhen, China). This drone integrates six 1 / 2.9-inch CMOS sensors, including one color sensor for visible light imaging and five monochrome sensors for multispectral imaging, each with 2.08 million effective pixels.

[0083] All drone data was collected under clear, cloudless conditions at midday (11:00 AM to 2:00 PM) with wind speeds below level 3. Pre-set flight paths were used, with flight altitudes set at 30m and heading and lateral overlap rates set to ≥80%.

[0084] The processing of the acquired image data, including image feature extraction, geometric correction, 3D reconstruction, point cloud generation, reflectivity correction, and surface model construction, was performed using Pix4Dmapper software (v4.4.10, Pix4D, Switzerland). Images acquired by the DJI Phantom 4 multispectral drone were processed using the accompanying DJI Terra software (v4.0.1, DJI, Shenzhen, China). The resolutions of the final digital orthophoto mosaics (DOM), digital surface models (DSM), and multispectral orthophoto mosaics (MOM) for Experiments 1, 2, and 3 are shown in Table 1.

[0085] Table 1. Sampling Dates of UAV Data and Measured Data

[0086]

[0087] (2.2) Satellite imagery data acquisition and preprocessing

[0088] Sentinel-2A (S2A) and PlanetScope (PS) satellite data covering the study area of ​​Field Experiment 3 were acquired on the dates closest to the UAV and ground measurements of Experiment 3. S2AMSI imagery was publicly released to the Google Earth Engine platform by the European Space Agency (ESA). In Google Earth Engine, the Sentinel-2A imagery data underwent preprocessing, including atmospheric correction, radiometric calibration, and surface reflectance calculation. This preprocessed data allows users to directly access high-quality remote sensing imagery and facilitates subsequent scientific research and application analysis. Image downloads and cloud cover removal were performed by accessing the Google Earth Engine Code Editor. PlanetScope spectral imagery was provided by Planet's Education and Research (E&R) program (Planet Team, 2017), and all bands were calibrated to a 3-meter spatial resolution.

[0089] (2.3) Plant Measurement

[0090] In Experiment 1, except for the sampling on August 22, 2021, which only measured 50 varieties (see Table 1), the aboveground fresh weight (AFW) and plant height of all sugar beet varieties were measured after each UAV flight. In Experiment 2, the aboveground fresh weight and plant height of all varieties were measured after the UAV survey. In Experiment 3, after each UAV survey, 15 points were randomly selected, and a 1-meter radius area centered on each point was taken as the sampling area to measure the aboveground fresh weight and plant height of the sugar beets. The latitude and longitude of each sampling were recorded using a smart RTK system (CHCNAV-i90, Shanghai, China). The plant height measurement method was to randomly select three sugar beets from the center of each plot / sampling area, and measure the height from the bottom to the highest point of the plant in its natural state. The specific sampling method for aboveground fresh weight was as follows: three sugar beets were randomly selected from the center of each plot / sampling area, and the aboveground parts were cut below the green leaf scar and weighed separately. The samples were then brought back and dried, and the average dry weight of the three samples was taken as the measured aboveground dry weight (ADW) of the plot / sampling area. The measured aboveground biomass (AGB) of the plot / sampling area was expressed as ADW multiplied by the planting density.

[0091] The plant height was extracted by subtracting the DSM from the bare land period from the digital surface model (DSM) generated by stitching together images from the above sampling periods. The pixel value of each point in the DSM image represents the elevation value of that point in the geographic coordinate system. Since no images of the bare land period were acquired in Experiments 1 and 2, the elevation values ​​of the bare land portion of the DSM image from the earlier sampling period were used as the ground elevation.

[0092] (2.4) Meteorological data collection

[0093] Meteorological data were obtained from the National Centers for Environmental Information (NCEIR) in the United States, extracting daily maximum and minimum temperatures from stations near the study area. Daily average temperatures were calculated by averaging these maximum and minimum temperatures. Then, based on the daily average temperatures, growing degree days (GDD) and relative growing degree days (RGS) were calculated.

[0094] (3) BBCH phenology of beets

[0095] In experiments 2 and 3, due to irrigation limitations and soil and topographical factors, sugar beets at different locations exhibited different phenological stages at the same time. In experiment 1, due to varietal varietal varietal variations, different sugar beet varieties showed different phenological stages at the same time. After obtaining images containing only the foreground using image processing techniques (background and vegetation segmentation based on a random forest regression classification model)... Figure 4 Coverage was calculated using the formula below, and then the BBCH (Biologische Bundesantalt, Bundessortenamt and CHemische Industrie, 0–99) codes were looked up to determine the phenology of sugar beets in each plot / sampling area. To ensure data consistency in the subsequent construction and validation of the biomass estimation model and to determine the model's applicability, the BBCH range (31–39) shared by all three datasets was considered, and data outside this range were removed to exclude them from model construction and validation.

[0096]

[0097] Where VF represents vegetation cover, and P A It is the total number of pixels, P V SR is the number of vegetation pixels and SR is the germination rate.

[0098] (4) Estimation of aboveground biomass based on allometric growth relationship

[0099] Allometric growth models typically include three mathematical forms: exponential growth models, logarithmic growth models, and power function growth models. Extensive research has shown a strong correlation between plant height and aboveground biomass for both crops and trees. Therefore, this study investigated the performance of the three allometric growth models using plant height estimated by the UAV in Experiment 1 as a variable, aiming to identify the optimal model. Furthermore, the study explored allometric growth models using spectral vegetation indices obtained by the UAV (Table 2) as variables, employing the coefficient of determination as an evaluation metric to select the optimal indices and models. Data from Experiments 2 and 3 were used to validate the models and evaluate their performance when extended to other habitats. The models were also compared with data-driven models such as random forest regression, partial least squares regression, and support vector machine regression.

[0100] After determining the optimal parameters and model, we further combined BBCH and the thermal index RGS to evaluate the performance of the model that combines phenology and thermal index, and extended its application to Experiments 2 and 3.

[0101] y = a * e b*V

[0102] y = a + b * ln(V)

[0103] y = a * V b

[0104] Where V represents plant height, vegetation index, and the product of plant height and vegetation index.

[0105] Table 2 shows the extracted canopy vegetation indices.

[0106]

[0107]

[0108]

[0109] Note: G, R, RE, and NIR are the reflectances of green, red, red edge, and near-infrared in their respective multispectral bands, respectively.

[0110] (5) Biomass estimation by matching field sampling data with UAV and satellite data

[0111] UAV data was collected and a biomass estimation model was constructed based on allometric growth relationships. Then, aboveground biomass mapping based on UAV data was performed. To match the scale of satellite imagery, the biomass map based on UAV data needed to be downsampled to match the 3m and 10m resolutions of Planet Scope and Sentinel-2A. Then, the downsampled UAV biomass map was sampled (see Section (6) for details) to construct a sample set for model construction, recording the latitude and longitude coordinates and biomass values ​​of the sampling points. Since the height of crops within 1 meter cannot be accurately estimated under the existing satellite data resolution, satellite spectral data was used as variables to construct a biomass estimation model based on satellite imagery. The reflectance values ​​of each satellite band corresponding to the sample point location and the biomass values ​​of the sampling points were used as training sets to construct the model.

[0112] Considering the limited area surveyed by the UAV in Experiment 3, a biomass estimation model based on Planet Scope imagery was first constructed. Then, the biomass and PS spectral values ​​at the ground measurement points in Experiment 3 were used as a validation set to verify the model's accuracy. Subsequently, based on the PS biomass estimation model, a biomass map based on PS imagery was generated for the sugar beet planting area outside the UAV's observation range. Then, using the same method, a biomass estimation model based on Sentinel-2A multispectral imagery was constructed based on this PS imagery biomass map. Background areas such as soil in the satellite imagery were removed in areas where the Renormalized Difference Vegetation Index (RDVI) was <= 0.33.

[0113] (6) SS-GMI, a method applicable to geospatial sampling

[0114] Due to resolution issues, ground-based sampling data is difficult to match satellite data, resulting in limited available sample data when modeling based on satellite data. After obtaining UAV-based biomass maps using the methods described above, sampling can be performed to expand the sample size. The most commonly used sampling method is random sampling. However, random sampling has several drawbacks: (1) data bias occurs due to concentration in certain areas; (2) random sampling may not adequately cover the entire geographic space, leading to a lack of representativeness; and (3) the samples exhibit spatial autocorrelation. Therefore, this study proposes the SS-GMI sampling method, which combines systematic sampling, Global Moran's I, and the optimal sampling strategy.

[0115] Systematic sampling, also known as interval sampling, is a sampling method that arranges N individuals in a population in a certain order, determines a random starting point according to rules, and then draws sample units one by one at regular intervals. Geospatial data is arranged in an orderly manner according to geographical location, so systematic sampling can avoid the sample being concentrated in certain areas and achieve coverage of the entire geographic space. Then, the step size is adjusted by combining the global Moran's index to ensure that the obtained sample set does not have spatial autocorrelation problems. The specific process of the SS-GMI sampling method is as follows: First, set the initial step size s, and start sampling points from the starting point (0,0) of the image with image size (m,n) to form a point set X.

[0116]

[0117] The number N of point set X is:

[0118]

[0119] After obtaining the point set X, calculate its global Moran index:

[0120]

[0121] Where N is the number of points in the point set, x represents the value of a point, x' represents the average value of the point values, and W' ... ij It is a spatial weight matrix.

[0122] When calculating the weight matrix W ij Eq.7 employs the K-Nearest Neighbors (KNN) weighted method, which offers advantages such as relative simplicity, high flexibility, and consideration of local spatial relationships. This method can capture the local spatial relationships around each point, making it suitable for many spatial analysis tasks that require consideration of local influences.

[0123]

[0124] To verify the significance of the global Moran's index, we check if its p-value is greater than 0.1. If p ≤ 0.1, we adjust the step size S = S * lc by the change coefficient lc, re-sample and recalculate the global Moran's index, updating I and p-values ​​until a step size satisfying p > 0.1 is obtained. Within the step size interval satisfying p > 0.1, we use a smaller change coefficient sc to find the minimum step size that satisfies p > 0.1, achieving the maximum number of samples, and output the final point set X. Spatial autocorrelation is caused by sampling neighboring points; therefore, the change coefficient lc can be set to a larger number to improve sampling efficiency.

[0125] After obtaining the final point set using the SS-GMI method, the spatial heterogeneity of the point set was determined using independent Geary's C and Getis's G indices. Based on a predetermined number of samples, a random sampling method was used to obtain the point set, and Geary's C, Getis's G, and global Moran's index were also used as evaluation metrics. To ensure fair comparison, ten random samplings were performed.

[0126] (7) Evaluation metrics for model testing

[0127] Select the coefficient of determination (R²) 2 The root mean square error (RMSE) and relative root mean square error (rRMSE) are used to evaluate the performance of the regression model.

[0128]

[0129]

[0130]

[0131] (8) Results Explanation

[0132] (8.1) Allometric relationship between plant height, vegetation index and aboveground dry weight

[0133] According to the data in Table 3, the allometric growth relationship between plant height and aboveground dry weight in Experiment 1 showed high accuracy in all three mathematical forms (R0). 2 ≥0.7). Among them, the power function has the highest accuracy, with a coefficient of determination of R. 2 =0.76, with a relative root mean square error (rRMSE) of 28%. Secondly, there is an exponential relationship, R... 2 =0.73. When vegetation index is used instead of plant height as the variable, the power function relationship between Dreg and aboveground dry weight shows the best performance, with a coefficient of determination of R. 2 =0.752; however, the accuracy decreased slightly when using the exponential relationship. Furthermore, the accuracy of the model was further evaluated when the product of plant height and Dre was used as a variable. The results showed that combining plant height and Dre improved the accuracy of all three allometric growth relationships, with the power function growth relationship performing best, achieving a coefficient of determination of R0. 2 =0.79.

[0134] Table 3 models the allometric growth relationship between plant height, vegetation index and its combination extracted from Experiment 1 and aboveground dry weight.

[0135]

[0136]

[0137] Note: PH represents plant height, and Dreg represents the vegetation index, which is the difference between red and green light.

[0138] Figure 5 , Figure 6 and Figure 7 Power function models of plant height, Dreg, and PH*Dreg versus aboveground dry weight are presented. As shown in the figures, a saturation problem exists when using plant height and Dreg as variables in the later stages of growth; that is, aboveground dry weight still increases significantly even when plant height and Dreg no longer increase significantly. This is particularly pronounced when using Dreg as a variable, possibly due to spectral saturation caused by canopy closure in the later stages of growth. Using the product of plant height and Dreg as a variable improves the saturation phenomenon observed when these two variables are used individually.

[0139] The best-performing power-law growth model from Experiment 1 was extended to beet experiments 2 and 3 for model validation on independent datasets. Figures 8-11 Since field trial 2 lacked spectral data, performance in trial 3 was only evaluated when Dreg and PH*Dreg were used as variables. Figure 10 , Figure 11 The allometric growth model of plant height and aboveground fresh weight derived from Experiment 1 performed well in Experiment 3, with R... 2 =0.87, rRMSE=20.04% Figure 9 However, this allometric growth model did not perform well in Experiment 2, R 2 =0.66, rRMSE=34.18% ( Figure 8 When using Dreg and PH*Dreg as variables to fit the allometric growth relationship with aboveground fresh weight, the performance in field trial 2 was not as good as the model using plant height as a variable. The estimation accuracy when using PH*Dreg as a variable was close to that of PH, with R0... 2 =0.87, rRMSE=20.59% Figure 11 ).

[0140] Further comparisons were made with machine learning regression models built using the same training set from Experiment 1 (Table 4). These trained machine learning regression models all performed well on the training set. Among them, the random forest regression model, regardless of whether PH, Dreg, or PH*Dreg were used as variables, significantly outperformed the partial least squares and support vector machine regression models (R0, R1, R2, R3, R4, R5, R6, R7, R8, R9 ... 2 ≥0.87). Similarly, these three machine learning regression models were also validated on independent datasets in Experiments 2 and 3 ( Figures 12-23When plant height was used as a variable, the predicted values ​​deviated significantly from the measured values ​​in both Experiment 2 and Experiment 3 during the early stages of sugar beet growth. Partial least squares regression, in particular, produced negative estimates in the early stages of sugar beet growth, severely deviating from reality, while its performance in the later stages of growth resembled a power function growth model. In Experiment 3, the model performance was generally poor when using Dreg or PH*Dreg as variables, except for the partial least squares regression model using the PH*Dreg variable, where the R-squared value between the predicted and measured values ​​in Experiment 3 was low. 2 The rRMSE reached 0.87, and the rRMSE was 20.24%. Figure 21 The performance of the power function growth model using PH*Dreg as a variable in Experiment 3 is close to that of the model. However, its accuracy is still slightly lower than that of the power function growth model using plant height as a variable. Figure 9 ).

[0141] Table 4 shows the performance of partial least squares, random forest, and support vector machine on the training set in machine learning regression models.

[0142]

[0143] Note: PH is plant height, and Dreg is the vegetation index representing the difference between red and green light.

[0144] (8.2) Power function growth model combining BBCH and relative growth days

[0145] Based on the aforementioned results, the model using the power function to estimate aboveground dry weight exhibits stronger robustness and scalability, outperforming overall the exponential and logarithmic allometric growth relationships and the three machine learning regression models. Therefore, this invention aims to investigate the impact of combining the BBCH phenological index and relative growth days on the robustness and scalability of the power function growth relationship. As shown in Table 5, after adding BBCH and relative growth days, the accuracy of the power function allometric growth model on the training set did not show a significant change, except that its accuracy decreased significantly when Dreg was used as a variable, R0. 2 The value was only 0.511, while the rRMSE was as high as 40%. Therefore, the model using Dreg as a variable was not extended to the independent datasets of Experiments 2 and 3 for validation. Figures 24-26 The results of using plant height and power functions with BBCH and relative growth days as variables were similar in Experiments 2 and 3, with R... 2 The values ​​are 0.8 and 0.83 respectively. (By comparing...) Figure 8 and Figure 9 It can be seen that the estimation accuracy of Experiment 2 was significantly improved, while the accuracy of Experiment 3 decreased slightly. The performance of Experiment 3 using PH*Dreg as a variable combined with BBCH and relative growth rate day (R0) is also shown. 2=0.82) is slightly less accurate than when only plant height is used.

[0146] Table 5 shows the allometric growth relationship between plant height, vegetation index and its combination, and aboveground dry weight after adding BBCH and relative growth rate.

[0147]

[0148] Note: PH is plant height, Dreg is the vegetation index of red edge and green light difference, BBCH is the phenological index, and RGS is the relative growth degree day.

[0149] (8.3) Mapping of aboveground biomass by UAV and sampling analysis of SS-GMI in Experiment 3

[0150] Based on the known plant density, the power function growth relationship between plant height and aboveground dry weight, y = 0.056 * PH, was used to determine the optimal growth rate for the three experimental results. 1.7866 The aboveground biomass (kg / m³) of the area surveyed by the UAV in Experiment 3 was obtained. 2 Spatial mapping Figures 27-29 As shown in the figure, the growth of sugar beets in this area is uneven, with the western side growing significantly better than the eastern side (see the supplementary direction in the figure). Furthermore, the areas with poorer growth at the beginning failed to catch up with the better-growing areas later on. Overall, the sugar beet growth in this area is not ideal.

[0151] Subsequently, the three UAV aboveground biomass maps were downsampled to match the PS satellite imagery (3m, 3m). The downsampled aboveground biomass maps were then sampled using the proposed SS-GMI sampling method. Figures 30-35 As shown in Table 6, the Moran's indices for all three sample sets are close to 0, and the p-values ​​of the significance t-tests are all greater than 0.1, indicating that spatial autocorrelation is not significant. Geary's C and Getis's G indices were further calculated and significance t-tests were performed to determine the spatial heterogeneity of the obtained sample point sets. The results show that the p-values ​​are all greater than 0.05, and the Geary's C values ​​are all close to 1, further indicating the independence of the data.

[0152] Table 6 Spatial autocorrelation assessment based on the SS-GMI sampling dataset

[0153]

[0154]

[0155] Note: p is the probability value.

[0156] Using a random sampling method, the same number of points were selected on the UAV ground biomass mapping, and this was repeated ten times. The global Moran index, Geary's C, and Getis's G index were then calculated. Figures 36-38 The results showed that the values ​​of the three indices and the p-values ​​obtained from the ten random samples across the three periods varied considerably. Regardless of the period, the p-values ​​for the significance tests of these three indices were consistently less than 0.05, indicating that the sample point sets obtained using the random sampling method exhibited a considerably high degree of spatial autocorrelation under these conditions.

[0157] (8.4) Modeling of aboveground biomass estimation based on satellite imagery

[0158] Using the aboveground biomass values ​​of sugar beets obtained by sampling using the SS-GMI method and the reflectance values ​​of each band of Planet Scope at the corresponding sampling points as the training set, a power function relationship between each band and the aboveground biomass values ​​of sugar beets was fitted. Among these relationships, the power function relationship between the red band of Planet Scope and the aboveground biomass of sugar beets showed the best performance, with R0. 2 Reached 0.83 ( Figure 39 Subsequently, we used the aboveground biomass values ​​of A at the actual ground sampling points in Experiment 3 and the Planet Scope spectral values ​​as a validation set to verify the accuracy of the model. R 2 =0.92, rRMSE is 12.89% ( Figure 40 Based on this, a map of aboveground sugar beet biomass based on PS imagery was generated in known sugar beet growing areas outside the observation range of nearby drones. Figures 41-43 As can be clearly seen from the figure, due to the limitations of the irrigation method, the sugar beets in the area covered by irrigation water, i.e., the circular irrigation area in the figure, grow significantly better than those in other areas.

[0159] Subsequently, using the SS-GMI method, aboveground biomass maps based on PlanetScope imagery were sampled to construct an aboveground biomass estimation model based on S2A imagery. Due to the date difference between the available Sentinel-2A imagery and the aboveground biomass maps based on PlanetScope imagery, the accuracy of the constructed model is generally R0. 2 =0.73 ( Figure 44 Based on a geospatial map of approximately 150 hectares of sugar beet growing areas provided by local farms, an aboveground biomass estimation model based on S2A imagery was used to invert the aboveground biomass of sugar beets at three time points in 2023 at Suqin Farm. Figures 45-47 As shown in the figure, there are low-value areas of aboveground biomass in both the early and late stages, and these areas are all located outside the irrigated area.

Claims

1. A method for large-scale farmland biomass monitoring based on UAVs and satellite imagery, characterized in that: The method is as follows: Step 1: Use drones to collect digital and multispectral images of multiple crop test areas at multiple time points, and perform preprocessing; Step 2: Collect satellite imagery data covering the areas surveyed by drones multiple times in Step 1, and perform preprocessing; Step 3: Measure the aboveground dry weight and plant height of all sugar beet varieties in the experimental areas; Step 4: Calculate the growth degree days and relative growth degree days for each experimental area based on meteorological data; Step 5: Based on the UAV preprocessing results, determine the common phenological zone of sugar beets in all test areas and estimate the sugar beet plant height, and verify the estimation accuracy using actual measured plant height; Step Six: After selecting the better variables and models using one of the test areas, combine phenological periods and heat index, and evaluate the performance of the model when applied to different test areas; Step 7: Biomass estimation by matching field sampling data with drone and satellite data: UAV data was collected, and a biomass estimation model was constructed based on allometric growth relationships. Then, aboveground biomass mapping was performed using the UAV data. The biomass map based on the UAV data was downsampled to match the 3m and 10m resolutions of PlanetScope and Sentinel-2A. The downsampled UAV biomass map was then sampled, as described in step eight, to construct a sample set for building a biomass estimation model based on satellite data. The latitude and longitude coordinates and biomass values ​​of the sampling points were recorded. Satellite spectral data was used as variables to construct a biomass estimation model based on satellite imagery. The reflectance values ​​of each satellite band corresponding to the sample point location and the aboveground biomass values ​​of the sampling points were used as a training set to construct the model. Step 8: The SS-GMI sampling method, which combines systematic sampling, the global Moran's index, and the optimal sampling strategy, is implemented by setting the initial step size. and from the image size Points are taken starting from the image's origin (0, 0) to form a point set. ; in, The initial step size, and In the image size The number of points in the horizontal and vertical directions of the image; dot set Number of for: Get point set Then calculate its global Moran exponent: in, Let x be the number of points in the point set, and let x represent the value of a point. This represents the average value of the points. It is a spatial weight matrix; When calculating the weight matrix The K-nearest neighbor weighting method was adopted, and the expression of the weight matrix is ​​as follows: Check if the p-value is greater than 0.

1. If p ≤ 0.1, adjust the step size S=S using the change coefficient lc. lc, re-select points and calculate the global Moran index, update I and p values ​​until a step size that satisfies the condition p > 0.1 is obtained; within the step size interval that satisfies p > 0.1, use a small change coefficient sc to find the minimum step size that satisfies p > 0.1, so as to achieve the maximum number of samples, and finally output the final point set X.

2. The method for large-scale farmland biomass monitoring based on UAV and satellite imagery according to claim 1, characterized in that: In step one, all UAV data were collected at noon on a clear, cloudless day with wind speeds less than level 3. The flight mission was carried out using a pre-set flight path, with the flight altitude set at 30 m and the heading and lateral overlap set at ≥ 80%.

3. The method for large-scale farmland biomass monitoring based on UAV and satellite imagery according to claim 1, characterized in that: In step one, the processing of the acquired image data includes image feature extraction, geometric correction, 3D reconstruction, point cloud generation, reflectivity correction, and surface model construction.

4. The method for large-scale farmland biomass monitoring based on UAV and satellite imagery according to claim 1, characterized in that: In step two, the preprocessing includes atmospheric correction, radiometric calibration, and surface reflectance calculation.

5. The method for large-scale farmland biomass monitoring based on UAV and satellite imagery according to claim 1, characterized in that: In step seven, UAV data was collected and a biomass estimation model was constructed based on allometric growth relationships. Then, aboveground biomass mapping was performed based on UAV data. To match the scale of satellite imagery, the biomass map based on UAV data needed to be downsampled to match the 3m and 10m resolutions of Planet Scope and Sentinel-2A. Then, the downsampled UAV biomass map was sampled to construct a sample set for building the biomass estimation model based on satellite data. First, a biomass estimation model based on Planet Scope imagery was constructed. Then, the aboveground biomass at actual ground measurement points and Planet Scope spectral values ​​were used as a validation set to verify the model's accuracy. Next, based on the Planet Scope biomass estimation model, aboveground biomass maps based on Planet Scope imagery were generated in sugar beet planting areas outside the UAV observation range. Then, using the same method, aboveground biomass estimation models based on Sentinel 2A and Landsat 8 multispectral imagery were constructed based on this Planet Scope imagery aboveground biomass map. Soil background areas in satellite imagery were removed in areas where the renormalized vegetation index (RDVI) <= 0.

33.

6. The method for large-scale farmland biomass monitoring based on UAV and satellite imagery according to claim 1, characterized in that: In step eight, the coefficient of determination is selected. Root mean square error and relative root mean square error Evaluate the performance of the regression model; 。

Images