Regional grain supply service temporal and spatial change analysis method based on long time sequence NDVI (normalized difference vegetation index)

Abstract

The invention discloses a regional grain supply service time-space change analysis method based on long time sequence NDVI, and particularly relates to the field of remote sensing, and the method comprises the following steps: S100, obtaining NDVI data and LUCC data of a research region, collecting regional grain yield data in yearbook, and downloading various ecological and social economic grids related to influencing grain supply; and step S200, extracting a target year cultivated land area from the land utilization grids, and deleting other areas to extract NDVI grids of cultivated land in each year of the research area. According to the method, grain supply service is distributed and grid-scale spatialization is performed on the grain supply service based on long-time-sequence NDVI, so that the spatio-temporal evolution analysis process of the grain supply service is visually displayed, potential influences of cultivated land quality, water area distribution, soil and management and the like on grain production service are considered, and meanwhile, the potential influences are taken as the basis of spatial analysis of subsequent driving factors of the grain supply service; and key science and technology can be provided for policy formulators to identify ecological service function damaged areas.

Description

Technical Field

[0001] The present invention relates to the field of remote sensing, and more specifically, to a method for analyzing the spatiotemporal changes of regional food supply services based on long-term NDVI. Background Art

[0002] Food supply is an important ecosystem service that can provide the necessary material basis for human survival. Quantifying regional food supply services, analyzing their spatiotemporal evolution, and identifying their driving factors are important foundations for ecosystem management and regulation, and are of great significance for the efficient allocation of resources and the environment and the rational formulation of conservation policies.

[0003] At present, the analysis of the distribution of grain production is mainly based on administrative scales, such as total grain output at the county and city levels and crop yield data per unit of cultivated land. Although this is relatively simple and quick, it lacks an intuitive display of the spatial distribution of regional grain supply services.

[0004] In order to solve the above problems, a technical solution is now provided. Summary of the invention

[0005] In order to overcome the above-mentioned defects of the prior art, an embodiment of the present invention provides a method for allocating and spatializing food supply services at a grid scale based on long-term NDVI, which not only intuitively displays the spatiotemporal evolution analysis process, but also takes into account the potential impact of cultivated land quality, water distribution, soil and management on food production services. It also serves as the basis for spatial analysis of subsequent driving factors, and can provide key scientific and technological means for policymakers to identify areas with impaired ecological service functions, so as to solve the problems raised in the above-mentioned background technology.

[0006] To achieve the above object, the present invention provides the following technical solutions:

[0007] The spatial and temporal variation analysis method of regional food supply services based on long-term NDVI includes the following steps:

[0008] Step S100, obtaining NDVI data and LUCC data of the study area and collecting regional grain production data from yearbooks, and downloading various ecological and socio-economic grids that affect grain supply;

[0009] Step S200, extracting the cultivated land area of ​​the target year from the land use grid, deleting other areas, extracting the NDVI grid of the cultivated land in each year of the study area, and completing the allocation of NDVI to the cultivated land grid;

[0010] Step S300, allocating the grain yield of the study area obtained from the yearbook to the grid of the cultivated land NDVI to obtain spatialized grain supply functional quantity data;

[0011] Step S400, after estimating the food supply service function quantity in each year according to the above steps, the trend analysis method is used to evaluate the change trend of the food supply service, and the Mann-Kendall rank correlation test is used to perform a significance test on the trend of the food supply service;

[0012] Step S500, select ecological and social factors that may drive regional food supply services, such as regional temperature, potential evapotranspiration, rainfall, relative humidity, night lights and land use, obtain the trade-offs and synergies between driving forces and changes in food supply services, and evaluate their significance and driving force ranking.

[0013] In a preferred embodiment,

[0014] Step S100 specifically includes the following contents:

[0015] Step S111, processing the downloaded MOD3Q1 data including the study area, and producing an annual NDVI grid based on the maximum value synthesis method;

[0016] Step S112, using the downloaded remote sensing image data with the vector boundary of the study area as a mask, cropping all the grids involved, and unifying the coordinates and resolution;

[0017] Step S113, calculate and record the total grain output in the region in the target year.

[0018] In a preferred embodiment, step S200 specifically includes the following contents:

[0019] Step S211, importing LUCC data into ArcGIS software, extracting cultivated land and setting its value to 1;

[0020] Step S212, using the obtained cultivated land rasters of each year, further extract the NDVI raster of the cultivated land in the target year.

[0021] In a preferred embodiment, step S300 specifically includes the following contents:

[0022] The formula for allocating food production is:

[0023]

[0024] Among them, GP i is the food supply function of the ith cultivated land grid; GP t represents the total grain output of the entire study area, in tons; n represents the total number of cultivated land grids in the study area; NDVI i is the annual NDVI value of the i-th cultivated land grid; NDVI min and NDVI maxare the minimum and maximum annual NDVI values ​​of cultivated land in the study area, respectively.

[0025] In a preferred embodiment, step S400 specifically includes the following contents:

[0026] Step S411, using trend analysis method to perform linear regression analysis on the food supply service that changes over time, so as to determine its changing trend. The trend analysis calculation formula is:

[0027]

[0028] Among them, S l ope is the slope of the pixel regression equation, GP k is the average value of the food supply function in the kth year, and t is the length of the study time; when S l ope>0, it means that the GP of the pixel is increasing; when S l ope=0, it means that the GP trend of the pixel is basically unchanged; when S l ope<0, it means that the GP of the pixel is decreasing;

[0029] Step S412, after trend analysis, a trend result grid is obtained, and a Mann-Kendall non-parametric rank test is performed on the trend result grid to obtain the significance result of the change in food supply services, and firstly, a statistical test quantity S is calculated;

[0030]

[0031] Among them, the value rules of sgn() are:

[0032]

[0033] Next, calculate the test statistic Z:

[0034]

[0035] Where t is the length of time for the study, x k and x j They are the food supply service values ​​corresponding to the kth year and the jth year respectively, j>k. When the given significance level is α=0.05, and the absolute value of Z is greater than 1.65, 1.96 and 2.58, it means that the trend has passed the significance test with confidence levels of 90%, 95% and 99% respectively.

[0036] In a preferred embodiment, step S500 specifically includes the following contents:

[0037] Step S511, processing temperature, potential evapotranspiration, rainfall, relative humidity, etc. into annual data;

[0038] Step S512, calculate the land use degree, which is one of the factors, and assign values ​​to different land use modes based on LUCC. The formula is as follows:

[0039]

[0040] Among them, La is the land use degree, A i is the index of the i-th land use type; C i is the area percentage of the i-th land use type. Different types of land are assigned different index values, reflecting the relative differences in the degree of utilization of different types of land. 耕地 =3, A 林地 =A 草地 =A 水域 =2, A 建设用地 =4, A 未利用地 =1;

[0041] Step S513, using spatial correlation analysis to analyze the correlation between the driving factors at the pixel scale and the changes in food supply services. The spatial correlation analysis formula is as follows:

[0042]

[0043] Among them, R xy is the correlation coefficient of variables x and y, x is the food supply service function, y is the driving force, t is the length of the study, and x i and i are the values ​​corresponding to the i-th year, and They represent the average values ​​of food supply service function and driving force respectively;

[0044] Step S514, apply terra and Hmisc packages in the RStud io software in turn to obtain the size of the driving coefficient, identify the overall ranking of the driving force according to its absolute value, and judge the overall significance of the driving force according to the P value, that is, when P<0.05, it is a significant trend feature.

[0045] The technical effects and advantages of the method for analyzing the spatiotemporal variation of regional food supply services based on long-term NDVI of the present invention are as follows:

[0046] By allocating and spatializing food supply services at a grid scale based on long-term NDVI, we can not only intuitively display the spatiotemporal evolution analysis process, but also consider the potential impact of arable land quality, water distribution, soil and management on food production services. It also serves as the basis for spatial analysis of subsequent driving factors, and can provide key scientific and technological means for policymakers to identify areas where ecological service functions are impaired. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1 It is a flow chart of the method for analyzing the spatiotemporal changes of regional food supply services based on long-term NDVI of the present invention.

[0048] Figure 2 This is the spatial distribution map of Qinling grain supply services in 2000, 2010 and 2020.

[0049] Figure 3 This is a chart showing the trends and significance of Qinling grain supply services.

[0050] Figure 4 It is a graph of the spatial correlation results and significance results of changes in ecological and social driving factors and food supply services. DETAILED DESCRIPTION

[0051] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0052] Example 1

[0053] Figure 1 The present invention provides a method for analyzing the spatiotemporal changes of regional food supply services based on long-term NDVI, which specifically includes the following steps:

[0054] Step S100, obtaining NDVI data, LUCC (land use / land cover) data of the study area and regional grain production data collected in the yearbook, and downloading various ecological and socio-economic grids that affect food supply;

[0055] Step S200, extracting the cultivated land area of ​​the target year from the land use grid, deleting other areas, extracting the NDVI grid of the cultivated land in each year of the study area, and completing the allocation of NDVI to the cultivated land grid;

[0056] Step S300, allocating the grain yield of the study area obtained from the yearbook to the grid of the cultivated land NDVI to obtain spatialized grain supply functional quantity data;

[0057] Step S400, after estimating the food supply service function quantity in each year according to the above steps, the trend analysis method is used to evaluate the change trend of the food supply service, and the Mann-Kendall rank correlation test is used to perform a significance test on the trend of the food supply service;

[0058] Step S500, selecting ecological and social factors that may drive regional food supply services, such as regional temperature, potential evapotranspiration, rainfall, relative humidity, night lights and land use, obtaining the trade-off and synergy results between driving forces and changes in food supply services, and evaluating their significance and driving force ranking;

[0059] Step S100 specifically includes the following contents:

[0060] Step S111, processing the downloaded MOD3Q1 data including the study area, and producing an annual NDVI grid based on the maximum value synthesis method;

[0061] Step S112, using the downloaded remote sensing image data with the vector boundary of the study area as a mask, cropping all the grids involved, and unifying the coordinates and resolution;

[0062] Step S113, calculate and record the total grain output in the region in the target year.

[0063] Step S200 specifically includes the following contents:

[0064] Step S211, importing LUCC data into ArcGIS software, extracting cultivated land and setting its value to 1;

[0065] Step S212, using the obtained cultivated land rasters of each year, further extract the NDVI raster of the cultivated land in the target year.

[0066] Step S300 specifically includes the following contents:

[0067] The formula for allocating food production is:

[0068]

[0069] Among them, GP i is the food supply function of the ith cultivated land grid; GP t represents the total grain output of the entire study area, in tons; n represents the total number of cultivated land grids in the study area; NDVI i is the annual NDVI value of the i-th cultivated land grid; NDVI min and NDVI max are the minimum and maximum annual NDVI values ​​of cultivated land in the study area, respectively.

[0070] Step S400 specifically includes the following contents:

[0071] Step S411, using trend analysis method to perform linear regression analysis on the food supply service that changes over time, so as to determine its changing trend. The trend analysis calculation formula is:

[0072]

[0073] Among them, S l ope is the slope of the pixel regression equation, GP k is the average value of the food supply function in the kth year, and t is the length of the study time; when S l ope>0, it means that the GP of the pixel is increasing; when S l ope=0, it means that the GP trend of the pixel is basically unchanged; when S l ope<0, it means that the GP of the pixel is decreasing;

[0074] Step S412, after trend analysis, a trend result grid is obtained, and a Mann-Kendall non-parametric rank test is performed on the trend result grid to obtain the significance result of the change in food supply services, and firstly, a statistical test quantity S is calculated;

[0075]

[0076] Among them, the value rules of sgn() are:

[0077]

[0078] Next, calculate the test statistic Z:

[0079]

[0080] Where t is the length of time for the study, x k and x j They are the food supply service values ​​corresponding to the kth year and the jth year respectively, j>k. When the given significance level is α=0.05, and the absolute value of Z is greater than 1.65, 1.96 and 2.58, it means that the trend has passed the significance test with confidence levels of 90%, 95% and 99% respectively.

[0081] Step S500 specifically includes the following contents:

[0082] Step S511, processing temperature, potential evapotranspiration, rainfall, relative humidity, etc. into annual data;

[0083] Step S512, calculate the land use degree, which is one of the factors, and assign values ​​to different land use modes based on LUCC. The formula is as follows:

[0084]

[0085] Among them, La is the land use degree, A i is the index of the i-th land use type; C i is the area percentage of the i-th land use type. Different types of land are assigned different index values, reflecting the relative differences in the degree of utilization of different types of land. 耕地 =3, A林地 =A 草地 =A 水域 =2, A 建设用地 =4, A 未利用地 =1;

[0086] Step S513, using spatial correlation analysis to analyze the correlation between the driving factors at the pixel scale and the changes in food supply services. The spatial correlation analysis formula is as follows:

[0087]

[0088] Among them, R xy is the correlation coefficient of variables x and y, x is the food supply service function, y is the driving force, t is the length of the study, and x i and i are the values ​​corresponding to the i-th year, and They represent the average values ​​of food supply service function and driving force respectively;

[0089] Step S514, apply terra and Hmisc packages in RStud io software in turn to obtain the size of the driving coefficient, identify the overall ranking of the driving force according to its absolute value, and judge the overall significance of the driving force according to the P value, that is, when P<0.05, it is a significant trend feature. It should be noted that, in general, when the P value is less than the given significance level (common significance levels are 0.05, 0.01, etc.): reject the null hypothesis and believe that the tested driving force is statistically significant, which means that the driving force has a real impact on the changes in food supply services and is not the result of random factors.

[0090] For example, if the P value is less than 0.05, you can usually say that the driving force is significant at a 95% confidence level.

[0091] When the P value is greater than or equal to the significance level: the null hypothesis is not rejected, that is, the driving force cannot be considered statistically significant. This may mean that the driving force has no obvious impact on the change of food supply services, or its impact may be masked by other factors, or it may be a result of the randomness of the sample, and it cannot be determined whether the driving force has an actual impact on food supply services that can be distinguished from random fluctuations.

[0092] In short, by comparing the P value with the pre-set significance level, we can judge the overall significance of the driving force, thereby determining which factors have played an important role in the changes in food supply services, and provide a basis for further analysis and decision-making.

[0093] Generally, if the P value is less than a pre-set significance level (such as 0.05), the driving factor can be considered statistically significant, that is, the factor has a practical impact on the phenomenon being studied.

[0094] (1) This embodiment selects remote sensing image data of the Qinling region in 2000, 2010 and 2020, processes and downloads MOD3Q1 data including the study area, and produces an annual NDVI grid based on the maximum synthesis method after splicing, projection and conversion;

[0095] (2) Download all temperature (mean temperature TMP, maximum temperature TMX, minimum temperature TMN), rainfall (PRE), potential evapotranspiration (PET), relative humidity (RH), and night light data (PANDA) for 2000, 2010, and 2020, use the vector boundary of the Qinling region as a mask, clip all involved rasters (including NDVI), and unify the coordinates to WGS_1984 with a resolution of 30 meters;

[0096] (3) Find the total grain output of the Qinling region in 2000, 2010 and 2020 from the yearbook and record it, with the unit converted to t;

[0097] (4) Import the LUCC data into ArcGIS software, use the raster calculator function of map algebra, use the conditional function Con() to extract cultivated land and set its value to 1;

[0098] (5) Using the extracted cultivated land grids of each year, the mask extraction function of ArcGIS extraction analysis was used to obtain the NDVI grids of cultivated land of each year;

[0099] (6) The grain yield of the study area obtained from the yearbook was allocated to the NDVI grid of cultivated land using the grid calculator function in ArcGIS through the grain yield allocation formula, thereby obtaining the spatialized food supply functional quantity data based on pixels. The unit of food supply service in each year was converted to t / ha according to the grid size, see Figure 2 ;

[0100] (7) The trend analysis formula was used to perform linear regression analysis on the food supply service that changes over time, so as to determine its changing trend. The spatial trend map was reclassified in ArcGIS, with the parts greater than 0 set to dark colors and the parts less than 0 set to light colors. Figure 3 a;

[0101] (8) The Mann-Kendall nonparametric rank test was performed on the raster data of the above trend results to obtain the significance results of the changes in food supply services. The spatial trend map was reclassified according to the value of Z in ArcGIS, see Figure 3b, judge their significance respectively;

[0102] (9) Select ecological and social factors that may drive the food supply service in the Qinling region, such as temperature, potential evapotranspiration, rainfall, relative humidity, night lights, and land use in the Qinling region in 2000, 2010, and 2020. Use spatial correlation analysis to analyze the correlation characteristics between the driving factors at the pixel scale and the changes in the food supply service in the Qinling region, identify the trade-offs and synergies between the driving forces and the changes in food supply services, and evaluate their spatial significance. Based on the driving coefficient between the driving forces and the changes in food supply services, identify the overall ranking and significance of the driving forces, see Figure 4 ;

[0103] (10) Taking various factors as driving forces, spatial correlation analysis was carried out using MATLAB software to identify their spatial driving effects on the grain supply service in Qinling Mountains, where positive coefficients were identified as synergistic relationships and negative coefficients were identified as trade-off relationships;

[0104] (11) Figure 4 The results show that La, PRE, and TMN are spatially synergistic with the food supply services in the Qinling Mountains, while PANDA and TMX are spatially trade-offs with the food supply services in the Qinling Mountains. PET, RH, and TMP are spatially trade-offs in the eastern Qinling Mountains and synergistic in the western Qinling Mountains. Among the above eight driving forces, TMN and TMX have a larger significance range in spatial driving. Therefore, temperature is identified as a more important driving force for the spatial change of food supply services in the Qinling Mountains.

[0105] (12) The overall ranking and significance of the driving forces of food supply services were completed using RStud io software. First, the terra and Hmisc packages were installed and called to read the raster data of the study area, merge the rasters and convert them into matrices for subsequent analysis;

[0106] (13) After deleting the na values ​​in the matrix, the Hmisc package was executed and the relationship between driving forces, the relationship between driving forces and regional food supply services, and the significance test were completed in the RStud io software. The driving factors were ranked according to the absolute value of the driving coefficient. The overall driving effect ranking was TMX (0.50) > TMP (0.47) > TMN = PRT (0.44) > PRE (0.40) > PANDA = RH (0.36) > La (0.33). The overall significance of the driving effect was P < 0.001, which was consistent with the results of spatial driving effect. Therefore, it was identified that temperature change was the most important driving force for the spatial change of food supply services in the Qinling region.

[0107] The above formulas are all dimensionless and numerical calculations. The formula is a formula for the most recent real situation obtained by collecting a large amount of data and performing software simulation. The preset parameters and thresholds in the formula are set by technicians in this field according to actual conditions.

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

[0109] Those of ordinary skill in the art will appreciate that the units and algorithm steps of each example described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are performed in hardware or software depends on the specific application and design constraints of the technical solution. Professional and technical personnel can use different methods to implement the described functions for each specific application, but such implementation should not be considered to be beyond the scope of this application.

[0110] In the several embodiments provided in the present application, it should be understood that the disclosed methods can be implemented in other ways. For example, the device embodiments described above are only schematic, for example, the division of the units is only a logical function division, and there may be other division methods in actual implementation. Another point is that the mutual coupling or direct coupling or communication connection shown or discussed may be through some interfaces, indirect coupling or communication connection of devices or units, which may be electrical, mechanical or other forms.

[0111] The units described as separate components may or may not be physically separated, and the components shown as units may or may not be physical units, and may be located in one place or distributed on multiple network units. Some or all of the units may be selected according to actual needs to achieve the purpose of the solution of this embodiment.

[0112] In addition, each functional unit in each embodiment of the present application may be integrated into one processing unit, or each unit may exist physically separately, or two or more units may be integrated into one unit.

[0113] If the functions are implemented in the form of software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present application, or the part that contributes to the prior art or the part of the technical solution, can be embodied in the form of a software product, which is stored in a storage medium and includes several instructions for a computer device (which can be a personal computer, a server, or a network device, etc.) to perform all or part of the steps of the methods described in each embodiment of the present application. The aforementioned storage medium includes: various media that can store program codes, such as a USB flash drive, a mobile hard disk, a read-only memory (ROM), a random access memory (RAM), a magnetic disk or an optical disk.

[0114] The above is only a specific implementation of the present application, but the protection scope of the present application is not limited thereto. Any person skilled in the art who is familiar with the present technical field can easily think of changes or substitutions within the technical scope disclosed in the present application, which should be included in the protection scope of the present application. Therefore, the protection scope of the present application should be based on the protection scope of the claims.

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

Claims

1. Analysis method of spatiotemporal variation of regional food supply services based on long-term NDVI, It is characterized in that The steps include: Step S100, obtaining NDVI data and LUCC data of the study area and collecting regional grain production data from yearbooks, and downloading various ecological and socio-economic grids that affect grain supply; Step S200, extracting the cultivated land area of ​​the target year from the land use grid, deleting other areas, extracting the NDVI grid of the cultivated land in each year of the study area, and completing the allocation of NDVI to the cultivated land grid; Step S300, allocating the grain yield of the study area obtained from the yearbook to the grid of the cultivated land NDVI to obtain spatialized grain supply functional quantity data; Step S400, after estimating the food supply service function quantity in each year according to the above steps, the trend analysis method is used to evaluate the change trend of the food supply service, and the Mann-Kendal l rank correlation test is used to perform a significance test on the trend of the food supply service; Step S500, select ecological and social factors that may drive regional food supply services, such as regional temperature, potential evapotranspiration, rainfall, relative humidity, night lights and land use, obtain the trade-offs and synergies between driving forces and changes in food supply services, and evaluate their significance and driving force ranking.

2. The method for analyzing the spatiotemporal variation of regional food supply services based on long-term NDVI according to claim 1, Features: Step S100 specifically includes the following contents: Step S111, processing the downloaded MOD3Q1 data including the study area, and producing an annual NDVI grid based on the maximum value synthesis method; Step S112, using the downloaded remote sensing image data with the vector boundary of the study area as a mask, cropping all the grids involved, and unifying the coordinates and resolution; Step S113, calculate and record the total grain output in the region in the target year.

3. The method for analyzing spatiotemporal changes of regional food supply services based on long-term NDVI according to claim 1, Features: Step S200 specifically includes the following contents: Step S211, importing LUCC data into ArcGIS software, extracting cultivated land and setting its value to 1; Step S212, using the obtained cultivated land rasters of each year, further extract the NDVI raster of the cultivated land in the target year.

4. The method for analyzing spatiotemporal changes of regional food supply services based on long-term NDVI according to claim 1, Features: Step S300 specifically includes the following contents: The formula for allocating food production is: Among them, GP i is the food supply function of the ith cultivated land grid; GP t represents the total grain output of the entire study area, in tons; n represents the total number of cultivated land grids in the study area; NDVI i is the annual NDVI value of the i-th cultivated land grid; NDVI min and NDVI max are the minimum and maximum annual NDVI values ​​of cultivated land in the study area, respectively.

5. The method for analyzing spatiotemporal changes of regional food supply services based on long-term NDVI according to claim 1, Features: Step S400 specifically includes the following contents: Step S411, using trend analysis method to perform linear regression analysis on the food supply service that changes over time, so as to determine its changing trend. The trend analysis calculation formula is: Among them, Slope is the slope of the pixel regression equation, GP k is the average value of the food supply function in the kth year, and t is the length of the study; Step S412: After trend analysis, a trend result grid is obtained, and a Mann-Kendal l non-parametric rank test is performed on the trend result grid to obtain the significance result of the change in food supply services.

6. The method for analyzing spatiotemporal changes of regional food supply services based on long-term NDVI according to claim 1, Features: Step S500 specifically includes the following contents: Step S511, processing temperature, potential evapotranspiration, rainfall, relative humidity, etc. into annual data; Step S512, calculate the land use degree, which is one of the factors, and assign values ​​to different land use modes based on LUCC. The formula is as follows: Among them, La is the land use degree, A i is the index of the i-th land use type; C i is the area percentage of the i-th land use type. Different types of land are assigned different index values, reflecting the relative differences in the degree of utilization of different types of land. 耕地 =3, A 林地 =A 草地 =A 水域 =2, A 建设用地 =4, A 未利用地 =1; Step S513, using spatial correlation analysis to analyze the correlation between the driving factors at the pixel scale and the changes in food supply services. The spatial correlation analysis formula is as follows: Among them, R xy is the correlation coefficient of variables x and y, x is the food supply service function, y is the driving force, t is the length of the study, and x i and i are the values ​​corresponding to the i-th year, and They represent the average values ​​of food supply service function and driving force respectively; Step S514, in RStudio software, apply terra and Hmisc packages in turn to obtain the size of the driving coefficient, identify the overall ranking of the driving force according to its absolute value, and judge the overall significance of the driving force according to the P value, that is, when P<0.05, it is a significant trend feature.

Images