Crop per unit yield statistical data downscaling method based on comprehensive agricultural condition index mean deviation

Abstract

The invention discloses a crop per unit yield statistical data downscaling method based on a comprehensive agricultural condition index standardization mean deviation. The method comprises the following steps: A1, data preparation; a2, performing agricultural condition parameter inversion; a3, comprehensive agricultural condition index construction; key phenological period comprehensive agricultural condition indexes are constructed, correlation between key phenological period crop agricultural condition parameters and statistical per unit yield is analyzed, key phenological period comprehensive agricultural condition indexes based on a normalized regression coefficient weight determination method are constructed, and precision verification is performed on the key phenological period comprehensive agricultural condition indexes based on per unit yield correlation; and A4, downscaling and precision analysis of per unit yield statistical data. The method comprehensively considers the correlation between the crop agricultural condition parameters in the key growth period and the crop per unit yield change, and also considers the influence of the spatial heterogeneity between the pixels of the comprehensive agricultural condition parameters on the regional per unit yield change at the regional level. The crop growth state spatial heterogeneity expression method based on the comprehensive agricultural condition parameter spatial-temporal dynamic weight is of great significance in improving the simulation capability of pixel scale crop unit yield spatialization change.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of software and information technology, and particularly relates to a crop yield statistical data downscaling method based on standardized deviation from mean of comprehensive agricultural condition indexes. BACKGROUND

[0002] Traditional crop yield statistical data is mostly tabular data with administrative units as statistical units. The statistical data can only reflect the differences between administrative regions, and the spatial differences of internal factors are hidden, and the real distribution and spatial variation of crop yield per unit area in geographical space have not been truly revealed. This is very unfavorable for the comprehensive spatial analysis and data fusion of crop yield statistical data with other natural, geographical, ecological and other spatial data on the spatial grid, thus greatly restricting the wide application of crop yield statistical data in geospatial analysis (Khan et al., 2010; Hu Yunfeng et al., 2011). Although the traditional method of directly assigning crop yield statistical values to administrative units can obtain crop yield statistical data spatialization information in the form of map patches, there is still the problem of homogeneity of crop yield statistical information within the map patches, and the contradiction between spatial heterogeneity of agricultural production factors and homogeneity of statistical data within the region, landscape scale of geographical-ecological processes and administrative unit scale of statistical data cannot be solved (Liu Zhong and Li Baoguo, 2012; Fan Yida et al., 2004). Therefore, with the continuous improvement of quantitative and precision requirements of modern scientific and technological research, the crop yield statistical data spatialization information obtained by traditional methods has been unable to continue to meet the precision requirements of the influence of climate change, resources and environment, natural disasters and other factors on food production system, and it is urgent to develop new technologies and methods for crop yield statistical data spatialization, so as to solve the bottleneck problem of obtaining high-precision spatialization information of agricultural production statistical data.

[0003] In addition, with the continuous application of artificial intelligence technology in the field of agricultural remote sensing in recent years, there have been many research results on crop yield estimation technology based on administrative unit crop yield statistics and machine / deep learning algorithms. However, with the development of science and technology, the administrative unit scale yield estimation results cannot meet the increasing demand of agricultural management departments for spatial resolution of crop yield estimation results, and it is urgent to carry out research on pixel scale large-scale crop yield prediction technology. However, in the pixel scale crop yield prediction and estimation based on machine / deep learning, in order to ensure the accuracy level of the crop yield prediction model, sufficient and complete crop yield sample data are needed as model training data to realize large-scale high-precision crop yield prediction at pixel scale. Because the number of yield samples collected on the ground is limited, it is difficult to meet the large sample requirements of machine / deep learning models, therefore, sample generation has become one of the key bottlenecks restricting the development of yield prediction technology to some extent. However, generating grid crop yield spatial distribution information from crop yield statistical data is one of the effective and important ways to obtain crop yield sample generation, therefore, it is urgent to carry out research on crop yield statistical data downscaling technology to meet the development needs of crop yield sample expansion based on pixel scale crop yield prediction and estimation technology.

[0004] Compared with the existing research results of spatialization of socio-economic statistical data (such as population, GDP, etc.), the spatialization of crop yield (yield per unit, total yield) statistical data, especially the spatialization of yield per unit statistical data, is relatively less in the research of spatialization of agricultural statistical data. The common method for spatialization of crop yield per unit statistical data is to directly assign the value, that is, to link the crop yield per unit statistical data with the administrative division data to form vector data, and then convert it into a certain scale of grid data. However, this method still cannot completely solve the problem of homogeneity of crop yield statistical data within the statistical unit. With the development of computer technology, information technology and spatial technology, in recent years, scholars at home and abroad have made some progress in the spatialization of crop yield per unit statistical data. For example, You et al. (2006) established a crop spatial allocation model in Brazil based on cross-entropy theory with the support of multi-source data such as land use, agricultural irrigation, crop rotation, arable land suitability, population density and crop distribution map, and obtained the spatial distribution information of the total yield of 9 crops such as wheat and corn at the scale of 10km*10km grid. Among them, the research on the spatialization of yield per unit statistical data is to use the proportion of the actual yield of cultivated land in the grid cell to the regional average potential crop yield to realize the allocation of crop yield per unit in the spatial grid of administrative unit, but this method of spatialization of crop yield per unit statistical data still mainly relies on experience or statistical model. In China, Shi Shuqin et al. (2011) used the method of combining zoning and statistical modeling to conduct spatialization research on crop yield per unit statistical data, that is, first establish a multivariate statistical model relationship between crop yield per unit statistical data and natural environmental factors such as terrain, climate and soil, and then realize the quantitative distribution of crop yield per unit statistical information in the spatial unit by using the values of each related factor in the statistical standard grid cell (such as 1km*1km grid), thereby providing data basis for the integration with geographic spatial data. However, this method is limited by the simulation accuracy of multivariate regression analysis method for simulating crop yield per unit spatial distribution information, which not only affects the accuracy of grid scale crop yield per unit simulation results, but also cannot effectively control the regional accuracy of crop yield per unit spatialization, so it is necessary to forcibly correct the crop yield per unit regression simulation results to meet the requirements of regional accuracy of crop yield per unit spatialization. Liu Zhong et al. (2012) distributed provincial grain yield to spatial grid based on population density as the dependent variable and land use data as auxiliary factors, and obtained the 1km*1km grid map of grain yield in China in 2000. Among them, the method for spatialization of provincial grain yield per unit statistical data mainly uses the linear model between provincial population density and unit arable land area grain yield to realize the spatialization of provincial grain yield per unit statistical data. Zhao Xueqing et al. (2023) used multivariate linear regression analysis and error correction method to realize the spatialization simulation of spring corn yield per unit in a large area of Northeast China based on 13 influencing factors such as soil, terrain and climate.

[0005] In general, there are few research results on the downscaling and spatialization of crop yield statistical data, especially the spatialization of crop yield statistical data. In the process of downscaling and spatialization of crop yield statistical data, the natural environmental conditions such as land use, arable suitability, terrain, climate, soil, and population density are mainly considered, and the crop growth parameters which have important influence on crop yield are not considered or only a single crop growth parameter (such as vegetation index) is considered to establish a correlation model between the statistical yield and the single crop growth parameter, so as to realize the downscaling and spatialization of regional crop yield statistical data. However, the single crop growth parameter has limited expression ability for the formation rule of crop yield, and it is difficult to fully represent the complex process of crop growth and development, which limits the accuracy of the downscaling method of crop yield statistical data to some extent. Therefore, in order to overcome the above shortcomings, it is necessary to carry out research on the downscaling of crop yield statistical data based on comprehensive agricultural indicators (CAI) constructed by comprehensively utilizing multiple growth periods and multiple agricultural indicators, so as to improve the accuracy and reliability of the spatial expression of crop yield statistical data. In addition, the existing research on the spatialization of crop yield does not consider the influence of the spatial heterogeneity of crop growth conditions within the administrative unit on the formation and yield level of regional crop yield, which hinders the further improvement of the spatialization accuracy of crop yield to some extent. SUMMARY

[0006] The technical problem to be solved by the present application is to provide a downscaling method for crop yield statistical data based on comprehensive agricultural indicators standardized deviation, which considers the temporal variation and spatial heterogeneity of crop growth conditions, and carries out research on the downscaling method for crop yield statistical data based on the standardized deviation of comprehensive agricultural indicators, so as to further improve the accuracy and level of spatialization of crop yield statistical data.

[0007] The technical solution of the present application is as follows:

[0008] A downscaling method for crop yield statistical data based on comprehensive agricultural indicators standardized deviation, comprising the following steps:

[0009] A1: data preparation; preparing ground measured data, remote sensing data and auxiliary data; the ground measured data includes crop growth parameters and yield ground measured data, crop leaf area index, crop coverage, crop aboveground biomass, crop canopy water content, crop canopy chlorophyll and crop ground measured yield data; the remote sensing data is Sentinel-2 multispectral satellite data; the auxiliary data includes crop distribution data and research area range vector data;

[0010] A2: Crop parameter inversion; firstly, crop parameter inversion is carried out based on the SNAP model; then, NPP is obtained based on MODIS and Sentinel-2 remote sensing data; finally, the accuracy of remote sensing inversion of crop parameters is verified by using ground measured crop parameters;

[0011] A3: Construction of comprehensive crop condition index; a comprehensive crop condition index at a key phenological period is constructed, and the correlation between crop condition parameters and statistical yield at the key phenological period is analyzed, a comprehensive crop condition index at the key phenological period based on a normalized regression coefficient weight determination method is constructed, and the comprehensive crop condition index at the key phenological period based on the yield correlation is verified in accuracy;

[0012] A4: Statistical data of yield per unit area is reduced in scale and accuracy is analyzed;

[0013] A4.1: Calculate the spatial and temporal comprehensive weight of the comprehensive crop condition index at the key phenological period;

[0014] A4.2: Calculate the spatial difference index of the comprehensive crop condition in the growth season;

[0015] A4.3: Statistical data of regional crop yield per unit area is reduced in scale;

[0016] A4.4: Accuracy verification.

[0017] In the method, in step A3, the correlation between the crop condition parameters at the key phenological period and the statistical yield is analyzed, and the correlation between the crop condition parameters at the key phenological period and the statistical yield is analyzed.

[0018] A monadic regression equation between each crop condition parameter at the key phenological period and the statistical yield of the crop is established, and a regression coefficient is determined as a weight measurement index of the crop condition parameter;

[0019] Y′=α i P i +b (5)

[0020] In the formula, Y' represents the ground measured yield per unit area (kg / ha); alpha i represents a monadic regression coefficient of the i-th crop condition parameter at a certain phenological period and the ground measured yield Y'; P i is a crop condition parameter at different phenological periods; b is a constant; alpha i The greater the value of alpha

[0021] In the method, in step A3, the comprehensive crop condition index at the key phenological period based on the normalized regression coefficient weight determination method is constructed, and the comprehensive crop condition index at the key phenological period based on the normalized regression coefficient weight determination method is constructed.

[0022] The normalized regression coefficient weight determination method is used to construct the comprehensive crop condition index of the key phenophase, that is, on the basis of the regression coefficient of the monadic regression equation between each crop condition parameter of the key phenophase and the statistical yield per unit, the regression coefficient and the comprehensive crop condition parameter are normalized respectively, so that the corresponding weight coefficient of each normalized crop condition parameter of the key phenophase is obtained, and then the above normalized regression coefficient is used as the weight to construct the comprehensive crop condition index of the phenophase.

[0023]

[0024] In the formula, P ik is the pixel value of the kth pixel of the ith crop condition parameter; P ik is the normalized pixel value of the kth pixel of the ith crop condition parameter; min(P ik ) is the minimum pixel value of the ith crop condition parameter; max(P ik ) is the maximum pixel value of the ith crop condition parameter; P i is the normalized ith crop condition parameter; F represents the comprehensive crop condition index of a certain phenophase; and m is the total number of crop condition parameters.

[0025] In the method, in step A4, the step 4.1 specifically includes the following steps.

[0026] A4.1.1 determining the regression coefficient of the comprehensive crop condition index of the key phenophase;

[0027] A4.1.2 calculating the difference coefficient of the comprehensive crop condition index;

[0028] A4.1.3 calculating the spatial and temporal comprehensive weight of the key phenophase.

[0029] In the method, the step 4.1.1 specifically includes the following steps: establishing a monadic regression equation between the comprehensive crop condition index of the key phenophase of the crop and the statistical yield per unit, determining the regression coefficient as the time weight measurement index of the comprehensive index of the crop in the phenophase;

[0030] Y′=β j F j +c (8)

[0031] In the formula, F j is the comprehensive crop condition index of the jth key phenophase; β j represents the regression coefficient of the comprehensive crop condition index of the jth phenophase and the ground measured yield per unit Y′, and the greater the value of β j , the more significant the influence of the comprehensive crop condition index on the ground measured yield per unit, and the higher the corresponding time weight coefficient; and c is a constant.

[0032] The method, the step 4.1.2 specifically includes the following steps: the difference coefficient of the comprehensive farmland condition index is calculated according to the following formula:

[0033]

[0034] d j =1-σ j (11)

[0035] In the formula, F jk The kth pixel value of the jth key phenological period comprehensive farmland condition index; F' jk The proportion of the kth pixel of the jth key phenological period comprehensive farmland condition index in the research area; n is the total number of crop pixels in the research area; sigma j The entropy value of the jth key phenological period comprehensive farmland condition index; d j The difference coefficient of the comprehensive farmland condition index in the research area.

[0036] The method, in step A4, the step 4.2 specifically includes the following steps: first, the standardized deviation from the mean of the comprehensive farmland condition index of each key phenological period is calculated on the pixel scale; then, the normalized spatiotemporal dynamic weight of the comprehensive farmland condition index of the key phenological period is multiplied by the standardized deviation from the mean to obtain the comprehensive farmland condition spatial difference index in the growth season.

[0037] According to the method of claim 1, in step A4, the step 4.2 specifically includes the following steps:

[0038] a) calculating the regional mean of the comprehensive farmland condition index of the key phenological period

[0039]

[0040] In the formula, The regional mean of the jth key phenological period comprehensive farmland condition index (j is 1, 2, 3, 4); F jk The value of the kth pixel of the jth key phenological period comprehensive farmland condition index; n represents the total number of crop pixels in the region to be spatialized;

[0041] b) standard deviation from the mean of the comprehensive farmland condition index

[0042]

[0043] In the formula, F jk The pixel value of the kth pixel of the jth key phenological period comprehensive farmland condition index; G jk The standardized deviation from the mean of the kth pixel of the jth key phenological period comprehensive farmland condition index;

[0044] c) based on the comprehensive farmland condition spatial difference index in the growth season

[0045]

[0046] In the formula, H k is the growth season comprehensive agricultural condition spatial difference index of the kth pixel.

[0047] In the method, in step A4, the step 4.3 specifically comprises the following steps: firstly, the present application assigns the yield per unit to each pixel; then, the regional average yield per unit is multiplied by the growth season comprehensive agricultural condition spatial difference index to obtain a yield per unit variation corresponding to the growth season comprehensive agricultural condition spatial difference index, and the regional average yield per unit is added to the corresponding yield per unit variation to obtain a yield per unit statistical data downscaling spatialization result, so that the conversion of the corn yield per unit statistical data from the administrative unit scale to the pixel scale is realized.

[0048] The beneficial effects of the present application are as follows:

[0049] (1) The present application proposes a corn yield per unit statistical data downscaling spatial expression method based on a spatiotemporal dynamic weighted comprehensive agricultural condition parameter construction. The method is based on the obtained multi-source remote sensing data in the main growth period, ground measured agricultural condition parameters, ground measured corn yield per unit and regional yield per unit statistical data, and utilizes the correlation between the key phenological period agricultural condition parameters and the yield per unit to construct a comprehensive agricultural condition index. On this basis, the spatial heterogeneity of the crop growth state is fully considered, the spatiotemporal dynamic weight of the agricultural condition index is calculated, and the yield per unit at the pixel scale relative to the regional statistical yield per unit is represented by the standardized deviation of the comprehensive agricultural condition index at the key phenological period, so that the construction of the comprehensive agricultural condition spatial difference index in the growth season is realized, and the downscaling method research of the regional yield per unit statistical data is completed. The method has important significance for promoting the development of crop statistical data downscaling spatialization.

[0050] (2) The present application utilizes the correlation between the key phenological period agricultural condition parameters and the yield per unit to complete the construction of the key phenological period comprehensive agricultural condition index based on the normalized regression coefficient weight, obtains the key phenological period comprehensive agricultural condition index, and provides important basic data for the yield per unit statistical data downscaling spatial expression research. The construction of the comprehensive agricultural condition spatial difference index in the growth season is carried out, the spatial heterogeneity expression of the crop growth state is completed through the determination of the spatiotemporal dynamic weight of the comprehensive agricultural condition parameter, the comprehensive agricultural condition spatial difference index in the growth season is realized, which has important significance for the research of the corn yield per unit statistical data downscaling spatial expression method. Finally, the yield per unit statistical data downscaling spatialization information of the corn is realized based on the comprehensive agricultural condition spatial difference index, and the conversion of the crop yield per unit statistical data from the regional administrative unit scale to the pixel scale is realized.

[0051] (3) The application carried out an application research of corn yield statistical data downscaling spatial expression based on spatio-temporal dynamic weighted comprehensive agricultural condition parameters in Hailun City, Heilongjiang Province. Through ground investigation data and crop yield statistical data verification, it is known that the regional accuracy of the crop yield statistical data downscaling method is close to 100%, the pixel scale accuracy, the ground measured yield and the spatialized pixel yield R 2 is 0.82, RMSE is 516.08 kg / ha, NRMSE is 13.15%, MRE is 4.75%, the above verification accuracy can meet the accuracy requirements of large-scale crop yield statistical data downscaling, proving that the crop yield statistical data downscaling spatialization method of the application has rationality and feasibility, and also showing that the yield downscaling spatialization method of the application has potential for application in a larger range. BRIEF DESCRIPTION OF DRAWINGS

[0052] Figure 1 is the technical roadmap of the application;

[0053] Figure 2 is the remote sensing inversion result of key agricultural condition parameters of the study area (2023);

[0054] Figure 3 is the remote sensing inversion accuracy verification result of key agricultural condition parameters of the study area (2023);

[0055] Figure 4 is the correlation between main agricultural condition parameters and crop yield at key corn phenological stages;

[0056] Figure 5 is the comprehensive agricultural condition index construction result at key corn phenological stages (2023);

[0057] Figure 6 is the comprehensive agricultural condition index accuracy verification at key corn phenological stages;

[0058] Figure 7 is the spatial distribution result of comprehensive agricultural condition spatial difference index of corn growth season (2023);

[0059] Figure 8 is the corn yield statistical data downscaling result in Hailun City (2023);

[0060] Figure 9 is the corn yield statistical data downscaling accuracy verification based on ground measured yield data. DETAILED DESCRIPTION

[0061] The application will be described in detail below in combination with specific embodiments.

[0062] REFERENCE Figure 1Under the support of sky-ground integrated observation technology, on the basis of obtaining multi-source remote sensing data in the main growth period, ground measured crop condition parameters, measured corn yield per unit and regional statistical yield per unit data, through the normalized regression coefficient weight determination method, considering the correlation between multiple crop condition parameters and yield per unit, a comprehensive crop condition index which can fully represent the crop growth and development status is constructed. Then, the influence of the dynamic change of crop growth status in time and space on the correlation between the comprehensive crop condition index and yield per unit is calculated to obtain the time and space dynamic weight of the comprehensive crop condition index at each key growth stage. After that, in order to more accurately describe the change of the pixel yield per unit relative to the regional average, the present application takes the pixel standardized deviation of the comprehensive crop condition index at the key growth stage as the research basis, and takes the time and space dynamic weight of the comprehensive crop condition index at the key growth stage as the main reference basis, and constructs a comprehensive crop growth condition spatial difference index (CSDICGC) in the growth season. On this basis, combined with the regional yield per unit average, the crop yield per unit value of each pixel in the administrative unit is obtained. Finally, the regional yield per unit statistical data is realized to be expressed in the form of spatial downscaling.

[0063] The following detailed description of the steps of the method of the present application:

[0064] Step one: data preparation;

[0065] 1.1 Ground measured data

[0066] 1.1.1 Crop growth parameters and yield per unit ground measured data

[0067] In order to carry out crop parameter remote sensing inversion and precision verification of crop yield per unit spatialization results, the present application carries out crop ground observation and investigation. The observation indexes mainly include crop leaf area index (LAI), canopy chlorophyll (CCC), canopy water content (CWC), fractional vegetation coverage (FVC), aboveground biomass (AGB), crop phenology information and crop yield per unit measured data, etc. Finally, a total of 5 ground investigations are carried out, and the main ground observation time and the main crop measured parameters are shown in Table 1.

[0068] Chlorophyll Content, CCC), canopy water content (Canopy Water Content, CWC), fractional vegetation coverage

[0069] (Fractional Vegetation Coverage, FVC), aboveground biomass (Aboveground Biomass, AGB), crop phenology information and crop yield per unit measured data, etc. Finally, a total of 5 ground investigations are carried out, and the main ground observation time and the main crop measured parameters are shown in Table 1.

[0070] Table 1 Ground observation time and main observation parameters

[0071]

[0072] 1.1.2 Crop leaf area index (LAI)

[0073] The present application uses LAI-2200 canopy analyzer to determine the leaf area index (LAI), which is a common optical equipment for measuring leaf area index. The device is portable and easy to operate, suitable for field environment. The measurement points are set at the lower part of the canopy, with the center of the row as the reference, and the displacement observation points are set at 1 / 4, 1 / 2 and 3 / 4 row distance along the vertical row. Each point is measured three times independently to eliminate random errors. The instrument calculates the canopy transmittance by measuring the top and bottom radiation values of the canopy, and inversely calculates the LAI according to the relationship between radiation attenuation rate and LAI and canopy structure. Finally, the arithmetic mean of the three measurement LAI data of each point is taken as the final determination value, which ensures the accuracy and reliability of the data.

[0074] 1.1.3 Crop coverage (FVC)

[0075] The present application uses unmanned aerial vehicle remote sensing technology to estimate the corn vegetation coverage (FVC). A quadcopter unmanned aerial vehicle platform is selected, which carries a visible light sensor. The relative flight height is set to 10 meters, and the vertical downward photography mode is used to obtain ground image data. The image acquisition time is selected at 10:00 to 14:00 local time, which can minimize the interference of terrain shadow and ensure the uniformity of light conditions in space and time. The unmanned aerial vehicle has the ability to collect corn canopy spectral information and can clearly present the corn vegetation coverage. Then, the obtained unmanned aerial vehicle image data is combined with algorithms such as pixel bisection model to divide the pixels and distinguish the vegetation and non-vegetation parts, so as to realize the accurate estimation of corn vegetation coverage. A large number of field tests have verified that this technology can quickly obtain the coverage information of large-area corn fields, and its precision is highly consistent with the requirements of agricultural production monitoring, providing a solid data support for corn growth evaluation and field management decision-making.

[0076] 1.1.4 Crop aboveground biomass (AGB)

[0077] In the present application, destructive sampling method is used for measuring aboveground biomass (AGB). In each sampling area, the part above the base of the corn stem is cut off as a sample, and is quickly sealed in a self-sealing bag, and then transported to the laboratory for subsequent analysis. In the laboratory, the stems, leaves, ears and other parts of the corn plants are separated one by one, and are respectively packed into paper bags and accurately weighed. After the fresh weight measurement is completed, each sample is placed in an oven, first treated at 105℃ for 2 hours, then the temperature is adjusted to 75℃ to dry to a constant weight state (mass change less than 0.5%), and the dry weight data is recorded in detail. Finally, combined with the planting density of the sampling area, the biomass per unit area is calculated to obtain accurate corn aboveground biomass data, which provides support for the accuracy verification of the subsequent key crop condition parameter net primary productivity (NPP).

[0078] 1.1.5 Crop canopy water content (CWC)

[0079] The present application adopts hierarchical sampling method to obtain canopy water data, the specific steps are as follows: in the target sample plot, randomly arrange sample points (≥3 repeats), respectively collect 3 representative leaves of the canopy. The sample is immediately sealed and stored in an aluminum foil bag after collection, and the sample point number and time are recorded. After the fresh weight (FW) is measured by using an electronic balance (accuracy 0.01g), the dry weight (DW) is obtained by first treating at 105℃ for 2 hours, then adjusting the temperature to 75℃ to dry to a constant weight state, and calculating the canopy water content (CWC=(FW-DW) / DW).

[0080] 1.1.6 Crop canopy chlorophyll (CCC)

[0081] In the field data collection, the present application selects corn plants with certain representative growth vigor at each sampling point, and respectively collects fresh leaf samples from the upper and middle parts, each group twice. Spectrophotometric method is used for chlorophyll content determination, that is, the leaf samples are stored in the dark after liquid nitrogen quick freezing, then extracted by ethanol-acetone mixed solution (volume ratio 1:1), the pigment extract is separated by centrifugation, and the absorbance values at 663nm and 645nm wavelengths are measured by ultraviolet visible spectrophotometer, and the chlorophyll a, b and total amount (Chl(a+b)) are calculated based on Arnon formula. This method verifies the data accuracy through the extinction coefficient and the standard curve, but attention should be paid to the sample grinding homogenization degree and the light protection treatment of the extract to avoid photodegradation error. The above CCC determination results provide ground truth data for subsequent key crop condition parameter CCC remote sensing inversion.

[0082] 1.1.7 Crop ground measured yield data

[0083] The application lays out observation sample points in Helen City on October 2, 2023 to carry out corn maturity period yield data collection work. In the crop yield measurement, first select 75 yield measurement sample points (of which 53 are used for comprehensive agricultural condition index precision verification, and 22 are used for crop yield statistical data downscaling pixel precision verification) in the field, select 3-5 corns at each sample point, dry naturally after taking the corn ear, and record the ear number, ear grain number and hundred grain weight. In addition, a 3m x 5m area is selected by tape measure, and the number of corn plants in each plot is measured to calculate the planting density. The crop yield measurement result of the measured point is calculated through the planting density, ear number and hundred grain weight.

[0084] 1.2 Remote sensing data acquisition and preprocessing

[0085] The application adopts the Sentinel-2 multispectral satellite data provided by the European Copernicus Data Center (https: / / scihub.copernicus.eu / ), which includes 13 band information from visible light to short wave infrared band, and the highest spatial resolution is 10m. The constellation is composed of two multispectral imagers (MSI) satellites, and the cooperative work of the double satellites can realize a global revisit period of 5 days. The Red-Edge Bands carried by the Sentinel-2 satellite provide reliable multispectral data support for dynamic monitoring and simulation of crop critical phenological period growth conditions. The advantages of high temporal and spatial resolution and the carried Red-Edge Bands provide an ideal data source for crop critical phenological period growth condition simulation (Clerici et al., 2017). The application processes the Sentinel-2 time series data in 2023 based on the Google Earth Engine (GEE) cloud platform. In the application, the key phenological period images with cloud cover less than 30% are screened, 1 scene of the best quality remote sensing image data is obtained for each key phenological period, and the QA band provided by the GEE platform is used to process the cloud mask of the obtained Sentinel-2 image (Hemmerling et al., 2021). At the same time, the linear interpolation method is used to fill the masked pixels (Kandasamy et al., 2013). The band and resolution of the Sentinel-2 remote sensing image used in the application are shown in Table 2.

[0086] Table 2 Main bands of Sentinel-2 applied in the application

[0087]

[0088]

[0089] 1.3 Auxiliary data

[0090] In the present application, the auxiliary data for the study of crop yield per unit statistical data downscaling spatialization includes crop distribution data, research area range vector data, etc. Among them, the research area Sentinel-2 cloud mask data is the QA band data provided by the GEE platform. The 30-meter corn crop distribution data in the research area comes from the National Ecological Science Data Center (https: / / doi.org / 10.57760 / sciencedb.08490.https: / / cstr.cn / 31253.11.sciencedb.08490.), which is mainly generated based on Landsat / Sentinel-2 NDVI fusion data set and time weighted dynamic time warping algorithm (TWDTW), and the total number of samples used is 54281. After verification, the overall accuracy is 80.06% on average; at the county level, the correlation coefficient (R2) between the identified area and the statistical area is between 0.657-0.903, and this data meets the precision requirements of crop mask data for crop yield statistical data spatialization.

[0091] Step two: crop condition parameter inversion;

[0092] 2.1 Crop condition parameter inversion

[0093] 2.1.1 Crop condition parameter inversion based on SNAP model

[0094] The present application uses the BiophysicalO module of SNAP (Sentinel Application Platform, SNAP) to carry out research on crop condition parameter inversion at key phenological stages. The model effectively combines the advantages of PROSAIL radiation transfer model and neural network algorithm, realizes the purpose of high-precision inversion of key crop condition parameters such as crop canopy chlorophyll, leaf area index, canopy water content and canopy coverage at multiple phenological stages, and uses ground measured data for precision verification, providing basic data for the construction of comprehensive crop condition parameters and the downscaling method of yield statistical data. The main crop condition parameter inversion and image date are shown in Table 3.

[0095] Table 3 Remote sensing acquisition of crop condition parameters based on Sentinel-2

[0096]

[0097] 2.1.2 NPP acquisition based on MODIS and Sentinel-2 remote sensing data

[0098] The present application is based on the existing biomass estimation method, adopts MODIS and Sentinel-2 remote sensing data to carry out 10m resolution NPP acquisition. Among them, through the fusion of 8-day PSNnet data of MODIS and 8-day maximum synthesis data of LAI remote sensing data inversion of Sentinel-2, the fine estimation of 8-day synthesis NPP with spatial resolution of 10m in crop growth period is realized, and the ground measured aboveground biomass data is used for accuracy verification. The main calculation formula is as follows:

[0099] NPP = PSNnet - Leaf_GR - Froot_GR (1)

[0100] PSNnet = GPP - Leaf_GR - Froot_GR (2)

[0101] Leaf_GR = 8day_leaf_mass_max * 8day_turnover_proportion * Leaf_gr_base (3)

[0102] Froot_GR = Leaf_GR * Leaf_gr_ratio (4)

[0103] In the formula, PSNnet has 8d data product in MOD17A2 product, Leaf_GR and Froot_GR in the application are leaf growth respiration and root growth respiration (kg C / 8day -1 ) respectively; 8day_leaf_mass_max is 8-day maximum leaf area index, 8day_turnover_proportion is a conversion coefficient, which is known from MOD17 product guide (https: / / lpdaac.usgs.gov / products / mod17a2hv006 / ) that the value of this coefficient is 1; Leaf_gr_base is the base value of leaf growth respiration, and the value is 0.30; Leaf_gr_ratio is the ratio of root growth respiration to leaf growth respiration, and the value is 2.0.

[0104] 2.2 Precision verification of crop condition parameter inversion

[0105] In the application, the ground measured crop condition parameters are used for precision verification of remote sensing inversion of crop condition parameters, in the application, the determination coefficient (Determination Coefficient, R 2Root mean square error (RMSE), Normalized root mean square error (NRMSE) and Mean relative error (MRE) are taken as the precision verification indexes of the agricultural condition parameters, and the specific calculation is shown in formulas (19) to (22) in 4.4.1.

[0106] In the present application, the number of ground survey samples for precision verification of the results of the agricultural condition parameters is 53; the reference results are the precision verification of the comprehensive agricultural condition indexes of the key phenological stages in 5.2.4.

[0107] Step three: construction of the comprehensive agricultural condition index;

[0108] 3.1 Construction of the comprehensive agricultural condition index of the key phenological stage

[0109] In the present application, four key phenological stages of the corn growth season, i.e., the jointing stage, the tasseling stage, the milk stage and the dough stage, are screened, canopy chlorophyll content (CCC), leaf area index (LAI), canopy water content (CWC), fractional vegetation coverage (FVC) and net primary production (NPP) are taken as the main agricultural condition parameters, and different weight coefficients are assigned to the agricultural condition parameters according to the correlation between the different agricultural condition parameters and the yield per hectare at each phenological stage, so as to construct the comprehensive agricultural condition index of the key phenological stage.

[0110] 3.2 Correlation analysis between the agricultural condition parameters of the crops at the key phenological stage and the statistical yield per hectare

[0111] A monadic regression equation between each agricultural condition parameter at the key phenological stage and the statistical yield per hectare of the corn is established, and the regression coefficient is determined as the weight measurement index of the agricultural condition parameter. In the present application, the number of ground survey samples for the correlation analysis between the agricultural condition parameters of the crops at the key phenological stage and the statistical yield per hectare is 53.

[0112] Y′=α i P i +b (5)

[0113] In the formula, Y' represents the ground measured yield per hectare (kg / ha); α i represents the monadic regression coefficient of the i-th agricultural condition parameter at a certain phenological stage and the ground measured yield per hectare Y'; P i is the agricultural condition parameter at different phenological stages; and b is a constant. iThe greater the value is, the stronger the influence of the crop condition parameter on the ground measured yield is, and the higher the corresponding weight coefficient is.

[0114] 3.3 Construction of the comprehensive crop condition index of the key phenophase based on the normalized regression coefficient weight determination method

[0115] The normalized regression coefficient weight determination method (NRCWDM) is proposed in the present application to construct the comprehensive crop condition index of the key phenophase, that is, on the basis of the regression coefficients of the monadic regression equation between each crop condition parameter and the statistical yield at the key phenophase, the regression coefficients and the comprehensive crop condition parameters are normalized respectively, so as to obtain the corresponding weight coefficients of the normalized crop condition parameters at each key phenophase, and then the normalized regression coefficients are used as the weights to construct the comprehensive crop condition index of the phenophase.

[0116]

[0117] In the formula, P ik is the pixel value of the kth pixel of the ith crop condition parameter; P ik is the normalized pixel value of the kth pixel of the ith crop condition parameter; min(P ik ) is the minimum pixel value of the ith crop condition parameter; max(P ik ) is the maximum pixel value of the ith crop condition parameter. P i ' represents the normalized ith crop condition parameter; F represents the comprehensive crop condition index of a certain phenophase; m is the total number of crop condition parameters, and m = 5 in the present application.

[0118] 3.4 Precision verification of the comprehensive crop condition index of the key phenophase based on the yield correlation

[0119] Since the comprehensive crop condition index in the present application is a comprehensive parameter representing the growth and development status of crops constructed by using multiple crop condition parameters in the main growth period, the construction result of the comprehensive parameter cannot be obtained in reality. Considering that there is a certain correlation between a single crop condition parameter (such as the leaf area index, the canopy water content, the canopy chlorophyll, the canopy coverage and the NPP) and the crop yield (Yao et al., 2022), the comprehensive crop condition index obtained from the above crop condition parameters should also have a certain correlation with the crop yield. Therefore, the correlation between the comprehensive crop condition index and the ground measured crop yield data is used as the basis for the indirect precision verification of the comprehensive crop condition index in the present application. Generally, the greater the correlation between the comprehensive crop condition index of the key phenophase and the measured crop yield is, the stronger the ability of the comprehensive crop condition index to describe the change of the crop yield is.

[0120] In the application, determination coefficient (Determination Coefficient, R 2 ), root mean square error (Root mean square error, RMSE), normalized root mean square error (Normalized root mean square error, NRMSE) and mean relative error (Mean relative error, MRE) are selected as the accuracy verification index of comprehensive agricultural condition index, and the specific calculation is shown in formula (19) to formula (22) of 4.4.1.

[0121] The application is used for comprehensive agricultural condition index verification of ground yield survey points, and the number is 53; the reference results and analysis in 5.2.1 correlation analysis between key phenological period agricultural condition parameters and crop yield.

[0122] Step four: yield statistical data downscaling and precision analysis

[0123] 4.1 Time and space comprehensive weight calculation of key phenological period comprehensive agricultural condition index

[0124] The application calculates the time weight by the correlation between key phenological period comprehensive agricultural condition index and yield, then calculates the space weight by the space information entropy of key phenological period comprehensive agricultural condition index. Finally, the time weight and space weight are normalized to obtain the time and space comprehensive weight of key phenological period comprehensive agricultural condition index.

[0125] 4.1.1 Regression coefficient determination of key phenological period comprehensive agricultural condition index

[0126] A regression equation between crop key phenological period comprehensive agricultural condition index and statistical yield is established, and the regression coefficient is determined as the time weight measurement index of the comprehensive index of the crop in the phenological period. In the application, the number of yield ground survey points for correlation analysis between single phenological period crop comprehensive agricultural condition index and statistical yield is 53.

[0127] Y′=β j F j +c (8)

[0128] In the formula, F j is the comprehensive agricultural condition index of the jth key phenological period; β j indicates the regression coefficient of the jth phenological period comprehensive agricultural condition index and ground measured yield Y′, and the greater the value of β j , the more significant the influence of the comprehensive agricultural condition index on the ground measured yield, and the higher the corresponding time weight coefficient; c is a constant.

[0129] 4.1.2 Comprehensive agricultural condition index difference coefficient calculation

[0130] To accurately reflect the spatial distribution characteristics of yield per unit area in the study area, the spatial weight coefficient of the comprehensive agricultural indicators was quantified based on the size of the difference coefficient in the Entropy Weight Method (EWM) to accurately characterize the spatial heterogeneity of yield distribution. Among them, the difference coefficient represents the spatial differences in regional crop growth conditions. The larger the difference coefficient, the greater the degree of uneven crop growth conditions within the region and the stronger its ability to express the spatial distribution of crop yield per unit area. Therefore, the greater the spatial weight represented. Conversely, the smaller the difference coefficient, the smaller the spatial weight represented. The calculation formula for the difference coefficient of the comprehensive agricultural indicators is as follows:

[0131]

[0132] d j =1-σ j (11)

[0133] Where, F jk The kth pixel value of the comprehensive agricultural index in the jth key phenological period; F′ jk is the proportion of the kth pixel of the comprehensive agricultural index in the jth key phenological period in the study area; n is the total number of crop pixels in the study area; σ j is the entropy value of the comprehensive agricultural index in the jth phenological period; d j is the difference coefficient of comprehensive agricultural indicators in a certain phenological period in the study area.

[0134] 4.1.3 Calculation of spatiotemporal weights of key phenological periods

[0135] a) Calculation of the time weight of key phenological periods

[0136]

[0137] Where C j is the time weight of the jth key phenological period; t is the total number of key phenological periods, and in the present invention, t=4.

[0138] b) Calculation of spatial weights of key phenological periods

[0139] Spatial weights are calculated based on the spatial difference coefficients of comprehensive agricultural indicators during key phenological periods. The larger the difference coefficient, the greater the unevenness of regional crop growth conditions, the stronger the ability to express the spatial distribution of crop yields, and the greater the corresponding spatial weight. The formula is as follows:

[0140]

[0141] Where D j is the spatial weight of the jth key phenological period.

[0142] c) Calculation of spatiotemporal comprehensive weights of key phenological periods

[0143]

[0144] In the formula, E j is the spatio-temporal normalized weight of the jth key phenological period.

[0145] 4.2 Calculation of the spatial difference index of the comprehensive crop condition in the growth season

[0146] The present application calculates the standardized deviation from the mean of the comprehensive crop condition index of the key phenological period by using the difference between the comprehensive crop condition index of the key phenological period and the regional mean. On this basis, the spatial difference index of the comprehensive crop condition in the growth season is obtained by multiplying the above-mentioned results of the four key phenological periods and the spatio-temporal comprehensive weight obtained in 4.1.3.

[0147] Under the influence of the crop growth environment, the crop growth and development conditions have certain spatial heterogeneity within the region. In order to accurately describe the difference of the crop growth and development conditions at the pixel scale and more accurately describe the spatial heterogeneity of the spatial distribution of crop yield, the present application proposes the concept of the spatial difference index of the comprehensive crop condition in the growth season. Firstly, the standardized deviation from the mean of the comprehensive crop condition index of each key phenological period is calculated at the pixel scale. Then, the spatial difference index of the comprehensive crop condition in the growth season is obtained by multiplying the normalized spatio-temporal dynamic weight of the comprehensive crop condition index of the key phenological period and the standardized deviation from the mean.

[0148] a) Calculation of the regional mean of the comprehensive crop condition index of the key phenological period

[0149]

[0150] In the formula, F represents the regional mean of the comprehensive crop condition index of the jth key phenological period (j∈1, 2, 3, 4); F jk represents the value of the kth pixel of the comprehensive crop condition index of the jth key phenological period; and n represents the total number of crop pixels in the region to be spatialized.

[0151] b) Calculation of the standardized deviation from the mean of the comprehensive crop condition index

[0152]

[0153] In the formula, F jk represents the pixel value of the kth pixel of the comprehensive crop condition index of the jth key phenological period; and G jk represents the standardized deviation from the mean of the kth pixel of the comprehensive crop condition index of the jth key phenological period.

[0154] c) Calculation based on the spatial difference index of the comprehensive crop condition in the growth season

[0155]

[0156] In the formula, Hk The growth season comprehensive agricultural condition spatial difference index of the kth pixel.

[0157] 4.3 Regional crop yield per unit statistical data downscaling

[0158] Since the regional yield per unit statistical data is the average yield per unit of the region, i.e., the average value of all pixel yield per units, first, the present application assigns the statistical yield per unit to each pixel. Then, the regional average yield per unit is multiplied by the growth season comprehensive agricultural condition spatial difference index to obtain the yield per unit variation corresponding to the growth season comprehensive agricultural condition spatial difference index, and the regional average yield per unit is added to the corresponding yield per unit variation to obtain the yield per unit statistical data downscaling spatialization result, thereby realizing the conversion of the corn yield per unit statistical data from the administrative unit scale to the pixel scale.

[0159] Finally, the ground yield per unit survey data is used to verify the accuracy of the crop yield per unit statistical data downscaling result.

[0160]

[0161] In the formula, Y is the statistical yield per unit of the study area (kg / ha); Y k Yk is the yield per unit of the kth pixel (kg / ha).

[0162] 4.4 Accuracy verification

[0163] 4.4.1 Accuracy verification based on ground measured data

[0164] In the present application, the ground measured yield per unit data is used to verify the accuracy of the yield per unit statistical data downscaling result. In the present application, the determination coefficient (Determination Coefficient, R 2 ), root mean square error (Root mean square error, RMSE), normalized root mean square error (Normalized root mean square error, NRMSE), and mean relative error (Mean relative error, MRE) are selected as the accuracy verification indexes of the crop yield per unit statistical data downscaling spatialization result, and the specific calculation is shown in formula (19) to formula (22).

[0165]

[0166] In the formula, Q r is the simulation value of the crop yield per unit statistical data downscaling result verification point; q r is the ground measured value corresponding to the crop yield per unit statistical data downscaling result verification point; is the average value of the simulation results of all verification points; The mean value of the corresponding all ground measured values; s represents the total number of measured data. In the precision verification index of the present study, the smaller the RMSE, the higher the precision. When NRMSE and MRE are less than 10%, it is judged that the precision of the simulation result is excellent, when NRMSE and MRE are greater than 10% and less than 20%, the simulation result is good, when NRMSE and MRE are greater than 20% and less than 30%, the simulation result is medium, and when NRMSE and MRE are greater than 30%, the simulation result is poor (Michele et al., 2003). When the precision evaluation is carried out, the precision evaluation standard gives priority to NRMSE, and then MRE. In addition, the consistency of the simulation value and the observation value can be measured by the determination coefficient (R 2 ) value, and the closer the R 2 value to 1, the better the consistency of the simulation value and the observation value, otherwise.

[0167] The number of ground yield survey points for crop yield spatialization result precision verification in the present application is 22.

[0168] 4.4.2 Regional verification based on crop yield statistical data

[0169] In addition to using ground measured yield data to verify the results of crop yield statistical data downscaling at the pixel scale, the present application also uses crop yield statistical data to verify the results of crop yield statistical data downscaling at the regional scale. Based on the construction of comprehensive agricultural condition parameters based on time and space dynamic weighting to obtain the results of statistical data downscaling of corn yield, the regional precision analysis and evaluation of the spatial distribution results of crop yield statistical data are carried out by using administrative unit crop yield statistical data.

[0170] In the present application, the administrative unit crop yield statistical data is used as the "mean true value" of the crop yield spatialization result verification The regional mean value (R) of crop yield obtained by statistical analysis of the results of spatialization of statistical data downscaling of corn yield based on the construction of comprehensive agricultural condition parameters based on time and space dynamic weighting is compared with , and the regional precision (w) in the administrative unit of the study area is obtained:

[0171]

[0172] Among them, R is the regional mean value result of the spatialization of crop yield in the administrative unit of the study area; Y k (kg / ha) is the yield value of the kth pixel of the results of spatialization of statistical data downscaling of crop yield; n is the total number of crop pixels in the study area; w is the regional precision of the results of spatialization of crop yield statistical data in the regional range; (kg / ha) is the crop yield statistical data in the administrative unit of the study area.

[0173] 5. Results and analysis

[0174] With the support of integrated sky-ground observation technology, the present invention utilizes multi-source remote sensing data, ground-based agricultural parameters, measured corn yields, and regional statistical yield data to propose a method for downscaling yield statistics based on comprehensive agricultural indicators of key phenological periods. In the present invention, comprehensive agricultural indicators of key phenological periods are first constructed, and the spatiotemporal weights of the comprehensive agricultural indicators of each key phenological period are calculated based on a spatiotemporal dynamic weighting method. Subsequently, the standardized mean deviations of the comprehensive agricultural indicators of key phenological periods are calculated at the pixel scale, and a weighted summation is used to construct a comprehensive agricultural spatial difference value index for the growing season. Finally, combined with regional statistical yield data, the spatial difference value of the comprehensive agricultural conditions during the growing season is used to infer crop yield information at the pixel scale, thereby obtaining the downscaling results of regional corn yield statistics.

[0175] 5.1 Inversion results and accuracy verification of main growth period parameters

[0176] In the study of corn yield prediction, agricultural parameters such as chlorophyll content (CCC), canopy cover (FVC), canopy water content (CWC), net primary productivity (NPP), and canopy chlorophyll content are important variables that characterize crop growth conditions. They can reflect the physiological state and biomass accumulation of crops from different angles, and play an important role in the spatialization of yield statistical data. At the same time, the remote sensing inversion accuracy of these parameters directly affects the reliability and accuracy of the yield statistical data downscaling results. The present invention uses the BiophysicalO module of SNAP and Sentinel-2 remote sensing images to invert agricultural parameters of crops in key phenological periods, and obtains agricultural parameters such as corn canopy chlorophyll, leaf area index, canopy water content, canopy cover and NPP in the key phenological periods of the study area. The specific results are as follows. Figure 2 As shown, this provides important agricultural parameter data for the construction of comprehensive agricultural parameters and the study of downscaling methods for yield statistics.

[0177] Depend on Figure 2 It can be seen that during the jointing stage, the values ​​of agricultural parameters were at a low level, with the average LAI value of 2.09m 2 / m 2 ; CCC value is 0-325.68μg / cm 2 The average value was 120.44 μg / cm 2 Although FVC reached 0.98 in a few areas in southern Hailun City, most areas were still at a low level, with an average FVC of 0.60. For CWC, the value of CWC in more than 90% of the areas was less than 0.1g / cm 2 Affected by the low temperature in June, the 8-day cumulative NPP value in some areas was low, and about 60% of the area was lower than 237.31gC / m 2LAI, CCC, FVC and CWC all reached the maximum value, and LAI of about 70% of the area reached 6.5 m 2 / m 2 ; FVC of about 90% of the area reached more than 0.9, and CCC of about 70% of the area reached 500 μg / cm 2 ; and CWC of about 40% of the area reached the average value of 0.15 g / cm 2 ; and the 8-day cumulative NPP of the area was at a relatively high level, with an average value of 272.52 g C / m 2 ; at the milking stage, the values of LAI, CCC and FVC were slightly lower than those at the tasseling stage, with average values of 2.42 m 2 / m 2 , 151.35 μg / cm 2 and 0.52, respectively; due to the influence of the previous rainfall, the value of CWC was at a relatively high level, with an average value of 0.14 g / cm 2 ; at the same time, due to the influence of high temperature in the first half of August, the 8-day cumulative NPP of Hailun City reached the highest value, and the 8-day cumulative NPP of about 90% of the area reached 228 g C / m 2 ; after reaching the ripening stage, part of the corn plants in the area entered the senescence stage, and the average value of LAI decreased to 1.74 m 2 / m 2 , the average value of CCC decreased to 107.54 μg / cm 2 , the average value of FVC decreased to 0.44 due to the influence of leaf senescence and dehydration, and the average value of CWC decreased to 107 g / cm 2 ; at the same time, due to the gradual decrease of temperature in September, the average value of 8-day cumulative NPP of part of the area in the central plain of Hailun City decreased to 239.8 g C / m 2 and below; the key phenological parameters obtained by remote sensing data inversion fully reflect the crop growth changes and spatial differences in the study area, and provide reliable data basis for the spatialization research of yield statistical data.

[0178] The present application uses the ground measured key agricultural parameters to verify the accuracy of the above-mentioned remote sensing inversion results, and the verification results are shown in Figure 3 and Table 4. Among them, the inversion accuracy R 2 of leaf area index is 0.66-0.80, and the RMSE is 0.30-0.69; the inversion accuracy R 2 of canopy water content is 0.66-0.78, and the RMSE is 0.01-0.02 g / cm 2 ; and the inversion accuracy R 20.68-0.89, RMSE 23.46-78.62 pg / cm 2 ; canopy coverage inversion accuracy R 2 0.62-0.79, RMSE 0.03-0.12; plant net primary productivity inversion accuracy R 2 0.66-0.78, RMSE 14.87-41.24 g C / m 2 ; from the NRMSE index, the NRMSE of the five key agricultural parameters of the application is 3.81%-26.60%, and the average NRMSE is 15.48%. It can be seen that the inversion accuracy of the overall agricultural parameters has reached a high level, and the accuracy of the agricultural parameter inversion result meets the subsequent application requirements.

[0179] Table 4 Remote sensing inversion accuracy verification results of key agricultural parameters

[0180]

[0181]

[0182] Note: The leaf area index unit is m 2 / m 2 , the canopy water content unit is g / cm 2 , the canopy coverage is dimensionless, and the canopy chlorophyll unit is pg / cm 2 , and the unit of NPP is g C / m 2 . In addition, ** indicates that it reaches a very significant level at P<0.01.

[0183] 5.2 Construction of comprehensive agricultural index of key phenological period

[0184] Crop growth and yield formation is a complex physiological and ecological process, which is affected by multiple factors and has obvious time sequence characteristics. Although the agricultural parameters inverted based on remote sensing data have certain correlation with crop yield per unit, the correlation between agricultural parameters at different phenological periods and yield per unit is different, and it is difficult to fully represent the spatial distribution of yield per unit by relying on a single agricultural parameter. Compared with a single agricultural parameter, the comprehensive agricultural index established by combining multiple agricultural parameters at key phenological periods can more accurately express the spatial distribution of yield per unit, thereby further improving the accuracy of crop yield per unit statistical data.

[0185] 5.2.1 Correlation analysis between key phenological period agricultural parameters and crop yield per unit

[0186] To determine the correlation between the different phenological parameters and the yield per hectare, the present application selects the key phenological periods, such as the jointing stage, the tasseling stage, the milk stage, the dough stage, and the like, to establish a linear regression model between the main agricultural parameters (such as the crop leaf area index, the canopy coverage, the canopy chlorophyll, the canopy water content, and the NPP, etc.) and the yield per hectare of corn. The regression equation and the accuracy are shown in Table 5. Figure 4 2 The present application uses the determination coefficient (R 2 ), the root mean square error (RMSE), the normalized root mean square error (NRMSE), and the average relative error (MRE) as the accuracy verification indexes of the linear regression results. The correlation between the key phenological period comprehensive agricultural parameters and the crop yield per hectare is shown in Table 5.

[0187] Table 5 Correlation analysis results between the main agricultural parameters of the key phenological period of corn and the crop yield per hectare

[0188]

[0189]

[0190] Note: The unit of the leaf area index is m 2 / m 2 , the unit of the canopy water content is g / cm 2 , the unit of the canopy coverage is dimensionless, the unit of the canopy chlorophyll is μg / cm 2 , and the unit of the net primary productivity is g C / m 2 . In addition, ** indicates that it reaches a very significant level at P < 0.01.

[0191] From Table 5, it can be seen that in the jointing stage, the correlation between the comprehensive agricultural indicators and the yield per hectare is NPP > LAI > FVC > CCC > CWC; in the tasseling stage, the correlation between the comprehensive agricultural indicators and the yield per hectare is LAI > NPP > CCC > FVC > CWC; in the milk stage, the correlation between the comprehensive agricultural indicators and the yield per hectare is LAI > NPP > CCC > FVC > CWC; and in the dough stage, the correlation between the comprehensive agricultural indicators and the yield per hectare is NPP > LAI > FVC > CCC > CWC. Among them, the correlation between the key phenological period agricultural parameters NPP and LAI and the yield per hectare of corn is better than the other three agricultural parameters, and the correlation of CWC is the weakest. At the same time, in the corn growing season, the correlation between the agricultural parameters LAI, NPP, FVC, CCC and CWC and the crop yield per hectare shows that the tasseling stage > the milk stage > the dough stage > the jointing stage, that is, the correlation between the agricultural parameters and the crop yield per hectare from the tasseling stage to the milk stage is the highest.

[0192] 5.2.2 Weight determination of the key phenological period agricultural parameters

[0193] To comprehensively compare the correlation between each agricultural condition parameter and the ground measured yield data, the comprehensive agricultural condition index is constructed. The application first normalizes each agricultural condition parameter in the key phenological period, and then allocates weight values according to the coefficient size of the monadic regression equation between the normalized agricultural condition parameters and the measured yield data, to obtain the weight coefficient of each agricultural condition parameter in the phenological period. The weight coefficient of the agricultural condition parameter in the key phenological period is shown in Table 6. From the table, it can be seen that the correlation coefficient between the agricultural condition parameter and the crop yield in different phenological periods (such as the jointing stage, the male flowering stage, the milk ripening stage, and the wax ripening stage) is higher, and the weight coefficient is higher, which indicates that the correlation between the agricultural condition parameter and the crop yield level in the corresponding phenological period is stronger. Among them, the minimum weight of LAI in the jointing stage is 0.1976, and the maximum weight in the male flowering stage is 0.2183; with the growth of crops, the weight of CWC reaches the highest 0.1894 in the milk ripening stage, and then decreases in the wax ripening stage. The weight of FVC is the highest in the jointing stage, which is 0.2317, and then gradually decreases, and the weight in the wax ripening stage is 0.1741, which reflects the importance of FVC in the early growth stage, and with the aging of leaves in the later growth stage, its contribution to the yield gradually decreases. The male flowering stage is the most active stage of photosynthesis, so the weight of CCC reaches the highest 0.2110 in the male flowering stage, and the lowest 0.1479 in the wax ripening stage; the weight of NPP reaches the highest 0.2882 in the wax ripening stage, and the weights in the male flowering stage and the milk ripening stage are relatively low, which are 0.2108 and 0.2126 respectively, which shows that the contribution of NPP to the formation of crop yield is more important in the later period, especially in the wax ripening stage.

[0194] Table 6 Weight determination results of main agricultural condition parameters in key phenological period of corn

[0195]

[0196] 5.2.3 Construction of comprehensive agricultural condition index in key phenological period

[0197] The agricultural condition parameters in each key phenological period are weighted and summed according to the corresponding weight coefficients in Table 6 to obtain the comprehensive agricultural condition index in the key phenological period. The specific results of the comprehensive agricultural condition index in the key phenological period are as follows: Figure 5It can be seen from the figure that in the jointing stage, the overall corn comprehensive farmland condition index in Hailun City is at a low level, with a regional mean of 0.39, mainly due to the low temperature in the northern part of Hailun City in June, which leads to a low value of the comprehensive farmland condition index in the northern part of Hailun City in the jointing stage, while the temperature in the southern part of Hailun City is normal and high in June, so the comprehensive farmland condition index in the southern part of Hailun City in the jointing stage is normal and high. In the tasseling stage, the temperature, precipitation and other conditions in the region from July to early August are suitable, so the comprehensive farmland condition index of corn in more than 90% of the regions in Hailun City reaches the highest value. In the milking stage, due to the influence of high temperature and rainy weather in the middle and late August in some regions, the comprehensive farmland condition index value in some regions begins to decrease. In the ripening stage, due to the influence of the gradual dehydration and aging of the leaves of corn in some regions in the middle and western parts of Hailun City, the comprehensive farmland condition index value in these regions decreases rapidly, and the comprehensive farmland condition index in some regions is at a low level.

[0198] 5.2.4 Precision verification of comprehensive farmland condition index in key phenological stage

[0199] The present application utilizes the correlation between the comprehensive farmland condition index in the key phenological stage and the ground measured yield data to establish a linear regression model, and adopts the determination coefficient (R 2 ), root mean square error (RMSE), normalized root mean square error (NRMSE) and mean relative error (MRE) as the precision verification indexes of the linear regression result, so as to represent the construction precision of the comprehensive farmland condition index in the key phenological stage. The correlation analysis results between the comprehensive farmland condition index in the key phenological stage and the ground measured yield data are shown in Table 6 and Table 7. Figure 6

[0200] As can be seen from Table 7, in the corn growing season, the correlation between the comprehensive farmland condition index in the key phenological stage and the crop yield reaches a highly significant level, and the correlation overall presents tasseling stage > milking stage > ripening stage > jointing stage. Among them, the correlation between the comprehensive farmland condition index in the tasseling stage and the milking stage and the ground measured yield data is the highest, the determination coefficient R 2 is 0.68 and 0.63 respectively, the RMSE is 587.39 kg / ha and 626.38 kg / ha respectively, the NRMSE is 7.29% and 7.77% respectively, and the MRE is 6.33% and 6.48% respectively, while the comprehensive farmland condition index in the jointing stage shows the weakest correlation. Compared with the correlation between the single farmland condition parameter (such as leaf area index, crown water content, crown chlorophyll, crown coverage, NPP, etc.) in the key phenological stage and the crop yield, the correlation between the comprehensive farmland condition index in the key phenological stage and the crop yield is higher. Figure 4 ​As shown in Table 5, the determination coefficient between the single agricultural condition parameter at the jointing stage and the crop yield per unit is between 0.13 and 0.33, and the determination coefficient between the comprehensive agricultural condition index of the key growth stage constructed by the application and the ground measured yield data is 0.37; the determination coefficient between the single agricultural condition parameter at the tasseling stage and the crop yield per unit is between 0.41 and 0.66, and the determination coefficient between the comprehensive agricultural condition index of the key growth stage constructed by the application and the ground measured yield data is 0.68; the determination coefficient between the single agricultural condition parameter at the milking stage and the crop yield per unit is between 0.40 and 0.59, and the determination coefficient between the comprehensive agricultural condition index of the key growth stage constructed by the application and the ground measured yield data is 0.63; the determination coefficient between the single agricultural condition parameter at the dough stage and the crop yield per unit is between 0.16 and 0.44, and the determination coefficient between the comprehensive agricultural condition index of the key growth stage constructed by the application and the ground measured yield data is 0.55. It can be seen that, compared with using a single agricultural condition parameter, using a comprehensive agricultural condition index can significantly improve the expression ability of the spatial distribution of crop yield per unit, and can better explain the dynamic change rule of the yield per unit at the key growth stage, and the model established by using the comprehensive agricultural condition index to represent the yield change is more reasonable.

[0201] Table 7 Correlation analysis between comprehensive agricultural condition index of key growth stage of corn and crop yield per unit

[0202]

[0203] 5.3 Calculation based on comprehensive agricultural condition spatial difference index in growth season

[0204] Under the influence of the change of crop growth environment, the spatial distribution of crop yield per unit has significant spatial heterogeneity. In order to accurately describe the time dynamic change in the process of forming crop yield per unit, and considering the spatial heterogeneity of crop yield distribution in the research area, the correlation between the comprehensive agricultural condition index and the measured yield data is taken as the time weight measurement index, and the difference coefficient of the comprehensive agricultural condition index in the research area is taken as the spatial weight measurement index of the yield distribution, on this basis, the standardized deviation of the comprehensive agricultural condition index at each key growth stage is calculated. Finally, the comprehensive agricultural condition spatial difference value index in the growth season is constructed.

[0205] 5.3.1 Time weight calculation of comprehensive agricultural condition index at key growth stage

[0206] Based on the normalized regression coefficient weight determination method, the regression coefficient of the monadic regression equation between the comprehensive agricultural condition index at the key growth stage in Table 7 and the measured yield data is normalized as the time weight of the comprehensive agricultural condition index at each growth stage, and the specific comprehensive agricultural condition index time weight is shown in Table 8. It can be seen that the time weight at the milking stage is the highest, which is 0.3334, indicating that the correlation between the comprehensive agricultural condition index at this growth stage and the yield per unit is the strongest, which further verifies the key role of the comprehensive agricultural condition index in the yield formation process.

[0207] Table 8 Time weight of comprehensive agricultural condition index of key phenological stage of corn

[0208]

[0209] 5.3.2 Calculation of spatial weight of comprehensive agricultural condition index of key phenological stage

[0210] In order to fully consider the spatial heterogeneity of regional crop yield level, the difference coefficient is calculated based on the establishment of the comprehensive agricultural condition index of key phenological stage, and the difference coefficient is used as the spatial weight measurement index of different key phenological stages. The greater the difference coefficient, the greater the unevenness of the growth of crops in the region, the stronger the expression ability of the spatial distribution of crop yield, and the greater the spatial weight represented. On the contrary, the smaller the difference coefficient, the smaller the spatial weight represented. Finally, the spatial weight of the comprehensive agricultural condition index of key phenological stage is shown in Table 9. It can be seen that the difference coefficient of the comprehensive agricultural condition index of the tasseling stage is the largest, reaching 0.0770, which indicates that the comprehensive agricultural condition index of the tasseling stage has greater heterogeneity in space, and therefore, the spatial weight of the comprehensive agricultural condition index of the tasseling stage is larger, reaching 0.4200. The difference coefficient of the comprehensive agricultural condition index of the jointing stage is the smallest, only 0.0257, which indicates that the comprehensive agricultural condition index of the jointing stage is more uniform in spatial distribution, and therefore, the spatial weight of the comprehensive agricultural condition index of the jointing stage is the smallest, only 0.1401.

[0211] Table 9 Spatial weight of comprehensive agricultural condition index of key phenological stage of corn

[0212]

[0213] 5.3.3 Calculation of spatio-temporal dynamic weight coefficient of comprehensive agricultural condition index of key phenological stage

[0214] In the downscaling of yield statistical data, accurately describing the correlation between the comprehensive agricultural condition index and the yield and the spatial heterogeneity of the yield distribution in the study area is the key to improving the spatialization accuracy of the downscaling of yield statistical data. The spatio-temporal dynamic weight coefficient of the comprehensive agricultural condition index of key phenological stage is finally obtained by comprehensively considering the time weight and the spatial weight of the comprehensive agricultural condition index of key phenological stage in the study area. The specific spatio-temporal dynamic weight coefficient of the comprehensive agricultural condition index obtained according to the time weight in Table 8 and the spatial weight in Table 9 is shown in Table 10. Through the calculation of the spatio-temporal dynamic weight of the comprehensive agricultural condition index of key phenological stage, it can be seen that there is spatio-temporal difference in the expression ability of the comprehensive agricultural condition index of different phenological stages to the spatial distribution of yield. Among them, the spatio-temporal dynamic weight shows that the tasseling stage has the largest contribution to the expression of the spatial distribution of yield, and the spatio-temporal weight value is 0.3356; the expression of the spatial distribution of yield by the milking stage and the waxing stage is at an intermediate level, and the spatio-temporal weight values are 0.2871 and 0.2000, respectively; the expression of the spatial distribution of yield by the jointing stage is the smallest, and the spatio-temporal weight value is 0.1773.

[0215] Table 10 Temporal and spatial dynamic weight coefficients of comprehensive agricultural condition indicators of key corn phenological stages

[0216]

[0217] 5.3.4 Calculation based on spatial difference index of comprehensive agricultural condition in growth season

[0218] Under the influence of crop growth environment, crop growth conditions have certain spatial heterogeneity within the region. In order to effectively capture the difference of crop growth conditions at pixel scale and more reasonably represent the difference of spatial distribution of crop yield, the present application assigns reasonable weights to comprehensive agricultural condition indicators at different phenological stages by constructing temporal and spatial dynamic weights. Then, the change amount of pixel scale yield relative to the regional average statistical yield is obtained by using the standardized deviation of comprehensive agricultural condition indicators at each key phenological stage, and finally the spatial difference index of comprehensive agricultural condition in growth season is obtained. The specific spatial difference index calculation result is shown in Figure 7 From the figure, it can be seen that the spatial difference index of comprehensive agricultural condition in growth season can directly reflect the spatial difference of crop growth in the study area. The overall difference of corn growth conditions in Hailun City in 2023 is small. Among them, it is known by statistics that about 63% of the crop growth conditions in the study area are higher than the average level of regional crop growth conditions, representing that the crop yield in the above 63% area will increase; about 37% of the crop growth conditions in the study area are lower than the average level of regional crop growth, representing that the crop yield in the above area will decrease.

[0219] 5.4 Descent and precision analysis of regional corn yield statistical data

[0220] In order to accurately describe the spatial distribution of regional yield, the present application obtains the yield change amount corresponding to the index by using the regional average yield and the spatial difference index of comprehensive agricultural condition in growth season. Then, the regional average yield and the yield change amount are added, so as to realize the conversion of corn yield statistical data from administrative statistical unit scale to pixel level spatial scale. The specific yield descent scale result is shown in Figure 8 Finally, based on the ground yield investigation sample data of the study area, the precision verification analysis of the spatial distribution result of crop yield descent scale is carried out, and the specific verification result is shown in Figure 9

[0221] As can be seen from Figure 8 , the spatial distribution of corn yield in Hailun City in 2023 is generally uniform. Among them, the corn yield in the south is generally high, while the corn yield in the west is slightly low. It is known by statistics that the overall level of corn yield in the region is high, and the corn yield in 63% of the region is higher than the average level of regional corn yield, and only 37% of the region is lower than the average level of regional yield. As can be seen from Figure 9 ​It can be seen that by comparing the corn yield statistics of the study area with the regional mean of the crop yield downscaling results in this study area, the regional mean of the spatialized results of corn yield in the study area in 2023 is 7791kg / ha, while the corn yield statistics of Hailun City in 2023 are 7761kg / ha. It can be seen that the regional accuracy is 99.6%, close to 100%. In addition, the spatialized results of the yield statistics downscaling of the present invention are verified by the ground-measured corn yield survey data. It can be seen that the determination coefficient R between the pixel values ​​of the ground yield verification points and the yield downscaling spatial results of the present invention is 100%. 2 It reaches 0.82, RMSE is 516.08kg / ha, NRMSE is 13.15%, and MRE is 4.75%. It can be seen that the spatial expression method of downscaling corn yield statistical data based on spatiotemporal dynamic weighted comprehensive agricultural parameters proposed in the present invention has certain effectiveness and feasibility, and can obtain a large-scale and high-precision spatial distribution result of downscaling corn yield statistical data. The above results can not only meet the demand for spatial information of crop yield statistical data in the study of the impact of climate change, resource environment, natural disasters, etc. on food production systems, but also provide effective basic yield spatial information support for the establishment and accuracy verification of large-scale crop yield estimation and prediction models, and also have important application value in the field of precision agriculture.

[0222] It should be understood that those skilled in the art can make improvements or changes based on the above description, and all such improvements and changes should fall within the scope of protection of the appended claims of the present invention.

Claims

1. A method for downscaling crop yield statistics based on the mean deviation of comprehensive agricultural indicators, characterized in that: The following steps are involved: A1: Data preparation; Prepare ground truth data, remote sensing data, and ancillary data; A2: Agricultural parameter inversion: First, agricultural parameters are inverted based on the SNAP model; then, NPP is obtained based on MODIS and Sentinel-2 remote sensing data; finally, the accuracy of remote sensing inversion agricultural parameters is verified using ground-based agricultural parameters; A3: Construction of comprehensive agricultural indicators: Construct comprehensive agricultural indicators for key phenological periods, analyze the correlation between agricultural parameters of crops in key phenological periods and statistical yields, construct comprehensive agricultural indicators for key phenological periods based on the normalized regression coefficient weighting method, and verify the accuracy of comprehensive agricultural indicators for key phenological periods based on yield correlation; A4: Downscaling and accuracy verification of yield statistics; A4.1 Calculate the spatiotemporal weights of comprehensive agricultural indicators during key phenological periods; A4.2 Calculate the spatial difference index of comprehensive agricultural conditions during the growing season; A4.3 Downscaling of regional crop yield statistics; A4.4 Accuracy Verification.

2. The method according to claim 1, characterized in that In step A3, the analysis of the correlation between the key phenological period crop parameters and the statistical yield is specifically as follows: Establish a univariate regression equation between each agricultural parameter in the key phenological period and the statistical yield of the crop, and determine the regression coefficient as the weight measurement indicator of the agricultural parameter; Y ′ = a i P i +b (5) Where Y ′ represents the measured yield per unit area (kg / ha); α i Indicates the relationship between the i-th agricultural parameter of a certain phenological period and the ground-measured yield Y ′ The univariate regression coefficient of i is the agricultural parameter of different phenological periods; b is a constant; α i The larger the value, the stronger the impact of agricultural parameters on the measured yield on the ground, and the higher the corresponding weight coefficient.

3. The method according to claim 1, characterized in that In step A3, the construction of the key phenological period comprehensive agricultural condition index based on the normalized regression coefficient weight determination method is specifically as follows: The normalized regression coefficient weight determination method is used to construct the comprehensive agricultural condition index of the key phenological period. That is, based on the regression coefficient of the univariate regression equation between each agricultural condition parameter of the key phenological period and the statistical yield, the regression coefficient and the comprehensive agricultural condition parameter are normalized respectively to obtain the corresponding weight coefficient of the normalized agricultural condition parameter of each key phenological period. Then, the above normalized regression coefficient is used as the weight to construct the comprehensive agricultural condition index of the phenological period. Where, P ik is the pixel value of the kth pixel of the i-th agricultural parameter; P ik ′ is the normalized pixel value of the kth pixel of the i-th agricultural parameter; min(P ik ) is the minimum pixel value of the i-th agricultural parameter; max(P ik ) is the pixel maximum value of the i-th agricultural parameter; P i ′ represents the i-th agricultural parameter after normalization; F represents the comprehensive agricultural index of a certain phenological period; m is the total number of agricultural parameters.

4. The method according to claim 1, wherein In step A4, step 4.1 specifically includes the following steps: A4.1.1 Determine the regression coefficients of comprehensive agricultural indicators during key phenological periods; A4.1.2 Calculate the coefficient of variation of comprehensive agricultural indicators; A4.1.3 Calculate the spatiotemporal weights of key phenological periods.

5. The method according to claim 4, characterized in that The step 4.1.1 specifically includes the following steps: establishing a univariate regression equation between the comprehensive agricultural index of the key crop phenological period and the statistical yield, and determining the regression coefficient as a time weight measurement index for the comprehensive index of the crop in that phenological period; Y ′ =b j F j +c (8) Where, F j is the comprehensive agricultural index of the jth key phenological period; β j represents the comprehensive agricultural index of the jth phenological period and the ground-measured yield Y ′ The regression coefficient, β j The larger the value of , the more significant the impact of the comprehensive agricultural indicators on the measured yield on the ground, and the higher the corresponding time weight coefficient; c is a constant.

6. The method according to claim 4, characterized in that The step 4.1.2 specifically includes the following steps: The calculation formula for the comprehensive agricultural index difference coefficient is as follows: d j =1-s j (11) Where, F jk The kth pixel value of the comprehensive agricultural index in the jth key phenological period; F j ′ k is the proportion of the kth pixel of the comprehensive agricultural index in the jth key phenological period in the study area; n is the total number of crop pixels in the study area; σ j is the entropy value of the comprehensive agricultural index in the jth phenological period; d j is the difference coefficient of comprehensive agricultural indicators in a certain phenological period in the study area.

7. The method according to claim 1, characterized in that In step A4, step 4.2 specifically includes the following steps: first, calculating the standardized mean deviation of the comprehensive agricultural indicators of each key phenological period at the pixel scale; then, multiplying the normalized spatiotemporal dynamic weight of the comprehensive agricultural indicators of the key phenological period by the standardized mean deviation to obtain the spatial difference index of the comprehensive agricultural conditions in the growing season.

8. The method according to claim 1, characterized in that In step A4, step 4.2 specifically includes the following steps: a) Calculate the regional mean of comprehensive agricultural indicators during key phenological periods Where, represents the regional mean of the comprehensive agricultural index in the jth key phenological period (j∈1, 2, 3, 4); F jk represents the value of the kth pixel of the comprehensive agricultural index of the jth key phenological period; n represents the total number of crop pixels in the area to be spatialized; b) Calculation of standardized deviation from the mean of comprehensive agricultural indicators Where, F jk is the pixel value of the kth pixel of the comprehensive agricultural index in the jth key phenological period; G jk is the standardized deviation from the mean of the kth pixel of the comprehensive agricultural index in the jth key phenological period; c) Calculation based on the spatial difference index of comprehensive agricultural conditions during the growing season Where H k is the spatial difference index of comprehensive agricultural conditions in the k-th pixel during the growing season.

9. The method according to claim 1, characterized in that In step A4, the step 4.3 specifically includes the following steps: first, the present invention assigns the statistical yield to each pixel; then, the regional average yield is multiplied by the comprehensive agricultural spatial difference index of the growing season to obtain the yield change corresponding to the comprehensive agricultural spatial difference index of the growing season, and then the regional average yield is added to the corresponding yield change to obtain the spatial downscaling result of the yield statistical data, thereby realizing the conversion of corn yield statistical data from the administrative unit scale to the pixel scale.