Multidimensional Assessment of Food Provisioning Ecosystem Services Using Remote Sensing and Agricultural Statistics

: With the increasing global population, human demands for natural resources continue to grow. There is a critical need for the sustainable use and development of natural resources. In this context, ecosystem services have attracted more and more attention, and ecosystem services assessment has proven to be useful for guiding research, policy formulation, and management implementation. In this paper, we attempted to assess ecosystem services more comprehensively from various perspectives. We used food provisioning ecosystem services in Minnesota as a case study and proposed two new concepts for assessing ecosystem services: e ﬃ ciency and trend. We designed a multidimensional assessment framework, analyzed the total output, e ﬃ ciency, and trend temporally based on both area and space with Exploratory Spatial Data Analysis (ESDA). We also identiﬁed major inﬂuencing factors based on remote sensing images in Google Earth Engine and explored the quantitative inﬂuence on each assessment dimension. We found that: (1) Food provisioning ecosystem service in Minnesota has generally been improving from 1998 to 2018. (2) We identiﬁed food provisioning ecosystem services in Minnesota as superior zones, mixed zones, and inferior zones with a ‘sandwich geo-conﬁguration’. (3) The total output tends to be stable while the e ﬃ ciency is disturbed by some natural disasters. Simultaneously, the trend index has been improving with slight ﬂuctuations. (4) Agricultural disaster ﬁnancial support has a stronger impact on stabilizing the total output of food provisioning than the other two dimensions. (5) Soil moisture, diurnal temperature di ﬀ erence, and crop growth are the three main inﬂuencing aspects of food provisioning ecosystem services, and the order of the inﬂuential density is: the Perpendicular Drought Index (PDI), Normalized Di ﬀ erence Vegetation Index (NDVI), Rainfall (RF), Daytime Temperature (DT), and Diurnal Temperature Di ﬀ erence (DIF).


Introduction
Ecosystem services are the various benefits provided to humans by the natural environment. With the increasing global population, human demands for natural resources continue to grow. There is a critical need for the sustainable use and development of natural resources. Various definitions of ecosystem services have been proposed in the literature. For example, Costanza et al. [1] defined ecosystem services as natures direct or indirect contribution to human welfare and classified them into 17 categories. Daily et al. [2] proposed that ecosystem services are the conditions and processes through which natural ecosystems, and the species that make them up, sustain and fulfill human life.
This study aimed to assess food provisioning ecosystem services comprehensively from various perspectives, building a framework in which the output, efficiency, and trend could be reflected simultaneously. Based on the theory of material quality evaluation, we believe that the amount of harvested crops can reflect the food provisioning ecosystem services. Compared with previous studies, the highlights of this paper are as follows: (1) We proposed two new concepts in ecosystem services assessment. One is the ecosystem service efficiency, which is the output of a unit ecosystem area, reflecting the quality of the ecosystem. The other is the ecosystem trend, reflecting the trend of ecosystem efficiency. The trend index is a vector whose direction and magnitude represent the trend and rate of change, respectively. (2) We designed an assessment framework to evaluate ecosystem services from three different aspects: total output, efficiency, and trend. (3) We analyzed the temporal changes of each dimension from the perspectives of area and geographic patterns. (4) We identified the main influencing aspects and analyzed the quantitative relationship between remote sensing indices and each assessment dimension using the Google Earth Engine cloud computing platform.

Study Area
We chose the state of Minnesota ( Figure 1) as our study area. Minnesota is the 12th largest in area, the 22nd most populous [53], and the 5th largest in agricultural production among all U.S. states [54]. Agriculture is the second largest industry in Minnesota, accounting for USD 57.5 billion in sales and more than 147,000 jobs [54].  Minnesota consists of 87 counties, covering 225,163 square kilometers. Given that the data of agricultural statistics are county-based, we performed the analysis at the county level.

Data Sources
We used statistical data from the U.S. Department of Agriculture (USDA) [55] to calculate multidimensional assessment indicators at the county level. Remote sensing images from the Google Earth Engine data catalog [56][57][58] were used to extract influencing factors. We selected corn, wheat, oat, and soybean as the main crops in Minnesota according to the Minnesota State Agriculture Overview [59].

Methods
We designed a multidimensional assessment framework and applied it to the 87 counties. Exploratory Spatial Data Analysis (ESDA) was adopted to explore the geospatial pattern of the results. We extracted remote sensing indices via Google Earth Engine and analyzed their quantitative relationships with the multidimensional assessment results. There were some statistical analysis methods used in this paper. The following sections describe the main methods we used.

A Multidimensional Assessment Framework
We proposed a multidimensional assessment framework ( Figure 2) that uses the total output (P), efficiency (Q), and trend (D) to uniquely characterize ecosystem services. It defines eight assessment spaces as shown in Table 1 Minnesota consists of 87 counties, covering 225,163 square kilometers. Given that the data of agricultural statistics are county-based, we performed the analysis at the county level.

Data Sources
We used statistical data from the U.S. Department of Agriculture (USDA) [55] to calculate multidimensional assessment indicators at the county level. Remote sensing images from the Google Earth Engine data catalog [56][57][58] were used to extract influencing factors. We selected corn, wheat, oat, and soybean as the main crops in Minnesota according to the Minnesota State Agriculture Overview [59].

Methods
We designed a multidimensional assessment framework and applied it to the 87 counties. Exploratory Spatial Data Analysis (ESDA) was adopted to explore the geospatial pattern of the results. We extracted remote sensing indices via Google Earth Engine and analyzed their quantitative relationships with the multidimensional assessment results. There were some statistical analysis methods used in this paper. The following sections describe the main methods we used.

A Multidimensional Assessment Framework
We proposed a multidimensional assessment framework ( Figure 2) that uses the total output (P), efficiency (Q), and trend (D) to uniquely characterize ecosystem services. It defines eight assessment spaces as shown in Table 1. In Figure 2b, objects 'a' and 'd' are located in Space V, 'b' is located in Space III, and 'c' is located in Space II. That is to say, 'a' and 'd' share the properties of Space V, and 'b' and 'c' share the properties of Space III and Space II, respectively. However, 'a' has stronger properties compared to 'd' because it has a larger assessment cuboid.  In Figure 2b, objects 'a' and 'd' are located in Space V, 'b' is located in Space III, and 'c' is located in Space II. That is to say, 'a' and 'd' share the properties of Space V, and 'b' and 'c' share the properties of Space III and Space II, respectively. However, 'a' has stronger properties compared to 'd' because it has a larger assessment cuboid.
The food humans obtain from nature mainly depends on the agriculture sector. Crop production is a branch of agriculture dealing with growing crops for use as food and fiber [60]. It is an effective indicator reflecting the total amount of food provisioning ecosystem services. Thus, the P-score in our research was derived from the main crop production. We defined ecosystem service efficiency as the ecosystem service output of a unit area, so the crop yield reflects the food provisioning ecosystem service efficiency. D-score is the trend index of food provisioning ecosystem service, calculated by the annual change in efficiency.
Due to the different value ranges of multi-source datasets used in this study, we standardized the data, i.e., rescaling the data to have a mean of 0 and a standard deviation of 1. This standardization is called a z-score, which can reflect the gap between data points and the overall average. A z-score can be calculated with the following formula: where x * is the z-score, x is the mean of the original data, σ is the standard deviation of the original data. Above-average data points will get a positive z-score, while subaverage data will get a negative z-score. The absolute value of the score indicates the number of standard deviations between the data points and the average. We used this z-score standardization method to calculate the P, D, and Q scores, which reflect the total output, efficiency, and trend of the food provisioning ecosystem services relative to the average level, respectively ( Table 2). Table 2. Formulas of assessment indices.

Index Formulation Introduction
P P i j = P i j −P j σ P j P i j is P-score of county i in the year j, P i j is the crop production of county i in the year j, P j is the average crop production of the state in year j, σ P j is the crop production standard variance of all counties in year j.
Q i j is the Q-score of ecosystem service efficiency of county i in the year j, Q i j is the yield of county i in the year j, Q j is the average crop yield of the state in year j, σ Q j is the yield standard variance of all counties in year j.
D i j is the D-score of the county i in the year j, ∆D j is the annual efficiency change in county i in year j, is the average of annual efficiency change in the entire state in year j, σ ∆D j is the annual efficiency change standard variance of all counties in year j.
The multidimensional assessment framework has the following characteristics: 1.
It is a relative assessment framework based on the average level. P-score, Q-score, and D-score are normalized to (−3,3). A positive index means an above-average dimension; the larger the index, the better the circumstance. Conversely, a negative index indicates a subaverage dimension; the lower the index is, the worse the circumstance is. 2.
P-score, Q-score, and D-score uniquely determine the characteristics of each county. The volume of the assessment cuboid reflects the strength of the characteristics it shows.

Exploratory Spatial Data Analysis
Initially, Tukey [61] proposed a new method of exploratory data analysis (EDA) to solve the problems in traditional statistics: (1) The practical data cannot meet the normal assumption.
(2) The models based on mean and variance tend to have low stability. EDA was able to analyze massive data, made no assumption of the population, and often excluded hypothesis testing. It tried to 'let the data speak' via statistical charts, graphs, and statistical overview methods. Exploratory Spatial Data Analysis (ESDA) is the expansion of EDA in the field of spatial data analysis. ESDA is supported by spatial analysis, emphasizes the spatial correlations of events, and focuses on the nature of spatial data. It is a powerful tool for geospatial data analysis, such as exploring spatial patterns, extracting main characteristics, identifying spatial clusters, evaluating aggregated and discrete patterns, grouping objects, and modeling spatial relationships. We adopted ESDA to explore the geographic spatial distribution and clustering patterns of food provisioning ecosystem services in Minnesota.
Hot-Spot Analysis, one of the main methods of ESDA, is based on the Getis-OrdGi* algorithm Equation (2) proposed by J. Keith Ord and Arthur Getis [62]. It focuses on identifying geospatial clusters of statistically high values (hot spots) and low values (cold spots). In a 'hot spots' zone, the high value is surrounded by high value. Similarly, a low value is located at the low-value gathering area in 'cold spots' zone.
where Gi * is positive for a 'hot spot' and negative for a 'cold spot'. The statistical significance is assessed by Z test in which the z-score and p-value are statistical significance measures and the accordances of whether to reject the null hypothesis or not (Table 3). In fact, they determine whether the observed spatial clusters are more obvious than we expected in a random distribution.

Remote Sensing Image Analysis
The modeling results, based on remote sensing data, could be very different from current conditions. Given that 2018 is the closest time period to the present, we believe that 2018 is more meaningful for recent policy-making and ecosystem management than other time periods (i.e., 1998,2003,2008,2013). Therefore, we selected 2018 as a case in modeling the influencing factors of food provisioning ecosystem services. We used remote sensing images acquired in 2018 to extract the influencing factors of the multidimensional assessment results. The typical crop growing season in Minnesota ranges from May to September [63][64][65]. To supplement some unusable satellite images acquired in May due to high cloud coverage, we extended the time range from April to September. Vegetation conditions and climatic conditions are the key factors affecting crop yield [45,47,66]. We chose the Perpendicular Drought Index (PDI) [47], Normalized Difference Vegetation Index (NDVI), Rainfall (RF), Daytime Temperature (DT), and Diurnal Temperature Difference (DIF) as the main influencing factor indicators ( Figure 3).
We derived the five influencing factors from remote sensing images at the pixel scale and then aggregated them to the county level. We used various remote sensing images in this paper and combined remote sensing data with vector data, resampling the pixels of remote sensing images from a micro perspective to a wider range of county boundaries to avoid the problems of different temporal and spatial resolutions of various remote sensing images. We used the USGS Landsat 8 Surface Reflectance Tier 1 (30-m resolution) data, which have been atmospherically corrected, to calculate the PDI and NDVI. The formulas are as follows: We derived the five influencing factors from remote sensing images at the pixel scale and then aggregated them to the county level. We used various remote sensing images in this paper and combined remote sensing data with vector data, resampling the pixels of remote sensing images from a micro perspective to a wider range of county boundaries to avoid the problems of different temporal and spatial resolutions of various remote sensing images.
We used the USGS Landsat 8 Surface Reflectance Tier 1 (30-m resolution) data, which have been atmospherically corrected, to calculate the PDI and NDVI. The formulas are as follows: where R red and R nir are the red band and near-infrared band of Landsat images, respectively. M is the slope of the soil line (Formula (5)), β is a constant.
DT and DIF data are from the dataset MOD11A2 Terra Land Surface Temperature and Emissivity 8-Day Global, which provides an average 8-day land surface temperature (LST). Daily Surface Weather and Climatological Summaries provide the RF data. These datasets are readily available in the Earth Engine public data catalog.

Statistical Analysis (1) Principal Component Analysis
Principal component analysis (PCA) is one of the most widely used statistical methods for dimensionality reduction. It was first introduced by K. Pearson [67]. It is the process of computing the principal components and using them to perform a change of basis on the data, sometimes using only Remote Sens. 2020, 12, 3955 8 of 19 the first few principal components and ignoring the rest [68]. We used the five influencing factors as the original dataset and extracted their principal components.
(2) Multiple Linear Regression Multiple linear regression attempts to model the relationship between two or more explanatory variables x i and a response variable y i by fitting a linear equation to observed data. The model for multiple linear regression, given N observations, is: In the least-squares model, the best-fitting line for the observed data is calculated by minimizing the sum of the squares of the vertical deviations from each data point to the line (if a point lies on the fitted line exactly, then its vertical deviation is 0). The restriction formula is: In this paper, x i indicates the five influencing factors; y i represents the multidimensional assessment indicators; the values fitted by Equation (6) are donatedŷ i .

Multidimensional Assessment Area Results
As shown in Figure 4, there are four dominant spaces in Minnesota (see Figure 2 for the eight spaces), including: progression zones with above-average output and efficiency (Space I), degradation zones with above-average output and efficiency (Space II), progression zones with below-average output and efficiency (Space VII), and degradation zones with below-average output and efficiency (Space VIII). These four dominant spaces accounted for more than 77% of the total area. Thus, our subsequent analysis focused on these four spaces (i.e., I, II, VII, and VIII).  The value of score indicates the level of research target compared to the average. A positive score means an above-average level and a negative score indicates a below-average level. We divided the results into six categories with an equal interval: (1) Scores greater than 2, indicating they are more than twice the standard deviation above the average level. (2) Scores in [1,2) are one to two times of standard deviation above the average level. P-score in Minnesota is mainly in (−2, 1), accounting for more than 80% of the total area, in which the negative P-score zones account for around 60%. It shows in 80% of the area in Minnesota that the output ranges from two standard deviations below average to one standard deviation above average, and the below-average part is dominating. Figure 5 shows that the area pattern of the P-score has been stable from 1998 to 2018. The negative P-score area has always been larger than the positive, The value of score indicates the level of research target compared to the average. A positive score means an above-average level and a negative score indicates a below-average level. We divided the results into six categories with an equal interval: (1) Scores greater than 2, indicating they are more than twice the standard deviation above the average level. (2) Scores in [1,2) are one to two times of standard deviation above the average level. P-score in Minnesota is mainly in (−2, 1), accounting for more than 80% of the total area, in which the negative P-score zones account for around 60%. It shows in 80% of the area in Minnesota that the output ranges from two standard deviations below average to one standard deviation above average, and the below-average part is dominating. Figure 5 shows that the area pattern of the P-score has been stable from 1998 to 2018. The negative P-score area has always been larger than the positive, indicating that the below-average output area has been larger than the above-average output area. The value of score indicates the level of research target compared to the average. A positive score means an above-average level and a negative score indicates a below-average level. We divided the results into six categories with an equal interval: (1) Scores greater than 2, indicating they are more than twice the standard deviation above the average level. (2) Scores in [1,2) are one to two times of standard deviation above the average level. (3) Scores in [0, 1) are within one standard deviation above the average level. (4) Scores in [−1, 0) are within one standard deviation below the average level. (5) Scores in [−2, −1) are one to two standard deviations below the average level. (6) Scores less than −2 are more than twice the standard deviation below the average level. P-score in Minnesota is mainly in (−2, 1), accounting for more than 80% of the total area, in which the negative P-score zones account for around 60%. It shows in 80% of the area in Minnesota that the output ranges from two standard deviations below average to one standard deviation above average, and the below-average part is dominating. Figure 5 shows that the area pattern of the P-score has been stable from 1998 to 2018. The negative P-score area has always been larger than the positive, indicating that the below-average output area has been larger than the above-average output area.
Similar to the P-score, more than 70% of the area in Minnesota has a Q-score within (−2, 1). That is to say, in 70% of the area, the food provisioning ecosystem service efficiency is from two standard deviations below average to one standard deviation above average, and this section increased to 91% in 2018. The above-average efficiency area (positive Q-score) is similar to the below-average efficiency area (negative Q-score) and the largest section is one standard deviation above-average area (0 < Qscore < 1), accounting for around 40% of the total ( Figure 5). More than 70% of the area had a D-score between (−1, 1), indicating that the trend index was very close to the average level. Figure 5 shows that the above-average trend index area (positive D- Similar to the P-score, more than 70% of the area in Minnesota has a Q-score within (−2, 1). That is to say, in 70% of the area, the food provisioning ecosystem service efficiency is from two standard deviations below average to one standard deviation above average, and this section increased to 91% in 2018. The above-average efficiency area (positive Q-score) is similar to the below-average efficiency area (negative Q-score) and the largest section is one standard deviation above-average area (0 < Q-score < 1), accounting for around 40% of the total ( Figure 5).
More than 70% of the area had a D-score between (−1, 1), indicating that the trend index was very close to the average level. Figure 5 shows that the above-average trend index area (positive D-score), the progression area, gradually overtook the below-average part (negative D-score), the degradation area. Thus, it is reasonable to believe that the food provisioning trend in Minnesota has improved. The areas of (1, 2) and (−2, −1) were gradually shrinking, which means more and more areas have been approaching the average trend level. In other words, the regional gap in the trend index has been decreasing.

Multidimensional Assessment Geography Spatial Results
The multidimensional assessment results were analyzed based on the nine agricultural districts as defined by the United States Department of Agriculture (USDA) ( Table 4):

Multidimensional Assessment of Spatial Analysis
We identified Minnesota in the inferior zones, mixed zones, and inferior zones according to the distribution of assessment spaces and found that they located as a 'sandwich geo-configuration' pattern. The superior zones were covered by progression zones with below-average output and efficiency (Space VII) and degradation zones with below-average output and efficiency (Space VIII), including North Central, Northeast, and East Central. The superior zones were dominated by progression zones with above-average output and efficiency (Space I) and degradation zones with above-average output and efficiency (Space II), spreading over West Central, Southwest, and South Central. The others were identified as mixed zones since they were mixed by the inferior zones and the superior ones. Figure 6 shows that the inferior zones and superior zones have stable temporal and spatial distribution but the mixed zones have changed a lot.

Total Output Geospatial Analysis
A 'sandwich geo-configuration' located in P-score distribution and in the geospatial cluster ( Figure 7): The inferior zones covered by the much below average P-score, [−3, −1), with the cold spots. The hot spots are located in the superior zones with the much above average P-score, [1,3]. The mixed zones did not have significant clusters and they were dominated by the around average

Total Output Geospatial Analysis
A 'sandwich geo-configuration' located in P-score distribution and in the geospatial cluster ( Figure 7): The inferior zones covered by the much below average P-score, [−3, −1), with the cold spots. The hot spots are located in the superior zones with the much above average P-score, [1,3]. The mixed zones did not have significant clusters and they were dominated by the around average P-score, [−1, 1).

Total Output Geospatial Analysis
A 'sandwich geo-configuration' located in P-score distribution and in the geospatial cluster ( Figure 7): The inferior zones covered by the much below average P-score, [−3, −1), with the cold spots. The hot spots are located in the superior zones with the much above average P-score, [1,3]. The mixed zones did not have significant clusters and they were dominated by the around average P-score, [−1, 1). Figure 7. Spatial distribution of P-score (total output of food provisioning ecosystem services). Figure 7. Spatial distribution of P-score (total output of food provisioning ecosystem services).

Efficiency Geospatial Analysis
The distributions of Q-score did not follow the 'sandwich geo-configuration' (Figure 8). The much below average Q-score was located in the inferior zones, while others were randomly spread over the superior zones and the mixed zones. However, the geospatial cluster in Q-score was similar to the P-score; cold spots were in the inferior zones, and hot spots were located in the superior zones. The hot spots were transferred to the southeast from the southwest. The distributions of Q-score did not follow the 'sandwich geo-configuration' (Figure 8). The much below average Q-score was located in the inferior zones, while others were randomly spread over the superior zones and the mixed zones. However, the geospatial cluster in Q-score was similar to the P-score; cold spots were in the inferior zones, and hot spots were located in the superior zones. The hot spots were transferred to the southeast from the southwest.

Trend Geospatial Analysis
The D score does not have a stable spatial distribution pattern, and it changes frequently from 1998 to 2018 (Figure 9). The hot spot analysis results showed that the geospatial cluster in Minnesota disappeared gradually, indicating that the trend index was geospatially randomized.

Trend Geospatial Analysis
The D score does not have a stable spatial distribution pattern, and it changes frequently from 1998 to 2018 (Figure 9). The hot spot analysis results showed that the geospatial cluster in Minnesota disappeared gradually, indicating that the trend index was geospatially randomized.

Trend Geospatial Analysis
The D score does not have a stable spatial distribution pattern, and it changes frequently from 1998 to 2018 (Figure 9). The hot spot analysis results showed that the geospatial cluster in Minnesota disappeared gradually, indicating that the trend index was geospatially randomized.

Principal Component Analysis of Influencing Factors
Based on existing studies [45][46][47] and the data availability, five indicators were selected in this paper, including the Perpendicular Drought Index (PDI), Normalized Difference Vegetation Index (NDVI), Rainfall (RF), Daytime Temperature (DT), and Diurnal Temperature Difference (DIF). We Figure 9. Spatial distribution of D-score (trend of food provisioning ecosystem services).

Principal Component Analysis of Influencing Factors
Based on existing studies [45][46][47] and the data availability, five indicators were selected in this paper, including the Perpendicular Drought Index (PDI), Normalized Difference Vegetation Index (NDVI), Rainfall (RF), Daytime Temperature (DT), and Diurnal Temperature Difference (DIF). We used the average value of remote sensing images during the crop growing season (from May to September). The principal component analysis divided these variables into three comprehensive indices and the principal component loads are shown in Table 5. The five factors influencing the food provisioning ecosystem service are classified into three types, i.e., soil moisture, diurnal temperature difference, and crop growth. It can be seen from Table 5 that the first principal component Z 1 has a large positive correlation with the variables related to soil moisture, the PDI, RF, and DT. Thus, the first principal component is regarded as a representative of soil moisture. The second principal component Z 2 has a large positive correlation with DNT, which is related to the temperature difference between day and night, so the second principal component is regarded as a representative of the day and night temperature difference. The third principal component Z 3 has a large positive correlation with NDVI, which is related to crop growth, so the third principal component is regarded as a representative of crop growth.

Multiple Linear Regression Analysis of Influencing Factors and Assessment Results
To explore the influence of each factor on different assessment dimensions, a multiple linear regression analysis based on assessment results was adopted: Compared with the F-test critical value table, all the F-test scores ( Table 6) we got are greater than the standard values at their corresponding confidence levels. Therefore, we believe that the linear regression calculation results are credible. The positive and negative D-scores have different regression equations, indicating that there are different ways the influencing factors affect the degrading food provisioning ecosystem service and the evolving one.

Multidimensional Assessment Results Discussion
Our results indicate that the food provisioning ecosystem service in Minnesota is significantly polarized, which is dominated by superior zones (Space I, Space II) and inferior zones (Space VII and Space VIII). The results also show that it has come a long way since the area of extra inferior zones (Space VIII) has shrunk considerably. While there is an exigent circumstance with degrading zones in Minnesota, which reached a portion of 50.21% in 2018. However, some improvements have been witnessed by the shrinking area of bad ecosystem service zones and the geospatially steady distributions of total output and efficiency, as well as the spatial homogenization of trend index.

Spatiotemporal Patterns of Multidimensional Assessment Results
The various dimensions have different development patterns. The P-score has developed steadily with a pretty good overall condition. It is witnessed by its range, (−2, 3), without extra-low value. However, the Q-score fluctuated sharply, especially in the last two five-year periods. The D-score has been always located in (−1,1), which means the trend index is relatively evenly distributed in the studied region. It is known that government agricultural disaster payments have irreplaceable positions in fighting weather extremes and ensuring food provisioning [69,70], whose distributions are affected by disastrous weather events. Thus, government agricultural disaster payments reflect the strength of agricultural disasters in the year. We noticed that the most agricultural disaster payments in Minnesota happened in the second period (2003)(2004)(2005)(2006)(2007)(2008)2008 is excluded), 303.48 million dollars, followed by the next five years, 199.60 million dollars [71]. Interestingly, the two sharp fluctuations in assessment results followed these two disaster periods, respectively. Thus, we believed that the disaster financial support stabilized the total output of food provisioning ecosystem service, but seemed unable to eliminate the impact on the efficiency and trend. This is consistent with the argument that agricultural disaster payments always subsidized market failure and focused on market returns [69]. At the same time, our findings bear a resemblance to the reasoning that subsidized disaster assistance and insurance may not prevent negative environmental effects from arising [72].
As the spatial patterns of multidimensional assessment results might be the consequences of ecological land distribution [73][74][75][76], it is reasonable that there is a significant boundary that separates Northeast, East Central, and North Central from the others in our results. This coordinates the reality that forests in Minnesota are mainly distributed in the northeast of the state, while cropland is mainly distributed in the southwest. The total output of food provisioning ecosystem service in Minnesota has steady temporally geospatial distribution and clustering patterns. However, there is a big difference between the amelioration in the northeast side of the state and the deterioration on the other side in ecosystem service efficiency. Its geospatial hotspot located in the south of the state has been moving toward the east, which is in line with the eastward population shifting in south Minnesota [77]. An obvious randomization in D-score reflects the spatially homogeneous development of food provisioning in the state of Minnesota.

Relationship between Multidimensional Assessment Results and Impacting Factors
Remote sensing vegetation indices (e.g., surface water indices, soil moisture, temperature condition, rainfall) have been proved to affect the crop yield [47,78,79]. In our research, we found that soil moisture, diurnal temperature difference, and crop growth are the three main influencing aspects of food provisioning ecosystem service. The influencing factors are ranked according to their strength as follows: PDI, NDVI, RF, DT, and DNT. It is worth noting that the negative coefficients of NDVI in the linear equations are the reflections of the extremely high NDVI of the forest zones with the weaker food provisioning ecosystem service, which is consistent with the results of Prasad's study of Iowa [47].

Limitations and Future Research Direction
Our analysis constituted the first step in multiperspective assessment in ecosystem services evaluation. Nevertheless, there are some limitations in our research: (1) We took food provisioning as the case study, and subsequent research could address other ecosystem services. (2) Referring to the existing studies [45][46][47], we only chose five remote sensing indicators as the influencing factors. More remote sensing products and indicators that may have an impact on food provisioning ecosystem services need to be further explored. (3) We adopted linear models to explore the relationship between impact factors and assessment results. With the rise of machine learning and artificial intelligence application in the field of ecosystem service evaluation [80][81][82][83][84], we may explore smarter multidimensional assessment models with higher precision in subsequent research. (4) We only selected remote sensing images acquired in the growing season in 2018 as a case to explore the ecosystem services influencing factors. Future research can further explore how the quantitative relationship varies according to the dates and seasons of remote sensing images.

Conclusions
Our multidimensional ecosystem services assessment framework evaluates the total output, efficiency, and trend of the ecosystem service simultaneously. It could identify the hidden and degrading ecosystems which are disguised by their good status quo as well as the low-quality ecosystems whose outputs are highly dependent on the area scales, thereby helping the environment management and sustainable development of the region. The serious polarization of food provisioning ecosystem service in Minnesota should be paid more attention, and relevant departments should take some targeted measures to reduce the area of Space VII and Space VIII. The degradation of some zones cannot be ignored either, some superior zones (Space II) are degrading, and a considerable area of inferior zones (Space VII) is getting worse. To ensure the sustainable development of food provisioning ecosystem services, these degradations must be slowed down. We analyzed the temporal development of every assessment dimension and found that the agricultural financial support only stabilized the total output of food provisioning and it worked a little on the efficiency and trend when facing natural disasters such as tornadoes and extreme weather. However, the determination of the quantitative influence relationships of the main influencing factors on each assessment dimension may help us think about how to fix the weak sides of ecosystem services.
To understand and assess the ecosystem services more comprehensively and explore their hidden aspects that are largely ignored, we tried to build a quantitative assessment framework for evaluating ecosystem services from multiple perspectives. We only selected food provisioning ecosystem service as a case in this paper. There is no doubt that more multidimensional assessment methods of other ecosystem services can be explored in future work. Funding: This research was supported by the Technology Innovation Center for Land Spatial Eco-restoration in Metropolitan Area, MNR, Shanghai, 200003. We would also like to thank the China Scholarship Council for sponsoring D.S. to study at the University of Tennessee as a visiting scholar.

Conflicts of Interest:
The authors declare no conflict of interest.