The historical data used in this article are climate and economic data from twenty nine provinces and autonomous regions in China from 1981 to 2018. The climate data include the monthly precipitation and average temperature of each province during the growing period (April to September), which were obtained from the National Meteorological Information Center (http://data.cma.cn). The economic data include the annual data on rural labour, grain sown area, agricultural fertilizer, grain yield and population. The data were obtained from the China National Statistical Center (http://www.bjstats.gov.cn) and the China Statistical Yearbook.

Future climate data include the monthly average of the 2020–2050 growth period (April-September) output by the BCC-CSM2-MR, which is newly developed by the National Climate Center. It is a mode resolution model participating in the sixth phase of the Coupled Model Intercomparison Project (CMIP6). The ability of the BCC-CSM2-MR climate model to simulate the climate mean state and global warming trend, quasi-biennial oscillation (QBO), tropical intraseasonal oscillation (MJO) and diurnal variation in precipitation in the 20th century has significantly improved. (Wu et al., 2019). CMIP6 designed a series of new scenario prediction experiments based on the possible anthropogenic emissions and land use changes caused by the energy structure in different shared socioeconomic paths. This paper mainly focuses on the three radiative forcing scenarios shown in Table 2 (Pu et al., 2020). These three different radiative forcing scenarios represent three future emissions levels, e.g., low, moderate and high, and provide a theoretical basis for studying the impacts of different climate changes on grain yield.

Grain yield is affected not only by the input of socioeconomic elements such as rural labour, grain cultivation area, and agricultural fertilizer but also by meteorological elements, including climatic elements. This article is based on the Cobb-Douglas, which can map the relationship between input and output under certain technical conditions. Chou et al. presented a novel approach by introducing the comprehensive climatic factor into the Cobb-Douglas model. It is called the economy-climate model (C-D-C model) (Chou and Ye, 2006). Then, we use grey system model GM(1,1) to predict the future socioeconomic data. GM(1,1) establishes a grey differential prediction model through a small amount of incomplete information to make long-term fuzzy information. It is of higher adaptability and accuracy (Yin, 2013). We use economy-climate and grey system models to construct a comprehensive model to predict the future grain yield using the following equations: Y it = μ ∗ x 1 β 1 x 2 β 2 x 3 β 3 C γ (1)

Its logarithmic function is as follows: ln ⁡ Y = β 1 ⁡ ln ⁡ x 1 + β 2 ⁡ ln ⁡ x 2 + β 3 ⁡ ln ⁡ x 3 + γ ⁡ ln ⁡ C + ln ⁡ μ (2) where Y_it represents grain yield; X₁, X₂, and X₃ represent the rural labour, the grain cultivation area, and the agricultural fertilizer in each province, respectively; β ₁, β ₂, β ₃, and γ are the coefficients of each element; and µ represents the other elements.

Principal component analysis (PCA) is probably the most commonly used dimensionality reduction method. To ensure that the original information is not missing or is missing as little data as possible, this method uses a few unrelated comprehensive indexes to express the majority of interrelated indexes (Brunetti et al., 2006). Let the original index temperature and precipitation be X ₁ and X ₂ and the main components be F₁ and F₂. { F 1 = a 11 X 1 + a 12 X 2 F 2 = a 21 X 1 + a 22 X 2 (3)

First, a sample matrix was constructed. Suppose a provincial administrative region in China has a sample size of monthly mean temperature precipitation from 1981 to 2018, and the observed climate indicator is temperature and precipitation. X = ( X 1 , X 2 ) T ( x 11 x 12 x 21 x 22 ⋮ ⋮ ⋮ x n 1 x n 2 ) (4)

Second, the principal component was determined. MATLAB was used to normalize the original climate data and calculate the relevant coefficient matrix r, principal component correlation coefficient e i j and eigenvalues of the matrix λ j θ j = λ j ∑ λ (5)

According to Eq. 5, the variance contribution rate θ j was obtained. Then, the correlation coefficient a i j of the i–th index in the j–th principal component was calculated according to Eq. 6. a i j = e i j / λ j (6)

Finally, the index weight was calculated, and the coefficient of each index was calculated according to Eq. 7. F j = ∑ j = 1 p a i j X i (7)

The weight vector was constructed, and Eq. 8 was used to calculate the weight of each index of the original data matrix. q = λ 1 F 1 + λ 2 F 2 + ⋯ + λ j F j ∑ i = 1 p λ i (8)

PCA was used to determine the weight of the comprehensive climatic factor C (Chou et al., 2019). It can represent the relative importance of a certain indicator. If the weight is greater, it means that it has a greater influence overall. C = q 1 X 1 + q 2 X 2 (9) where X ₁ and X ₂ are the temperature and precipitation, respectively, q 1 and q 2 are the weights obtained by PCA calculation. Although the comprehensive climate factors of the different provinces were calculated in the same way, their weights in the principal component analysis were differed due to various factors, such as terrain and climate, that lead to differences in the comprehensive climate factors of each province.

Resilience in physics refers to the ability of an object to resist the deformation of external forces and spring back into shape. Ecologically, resilience is defined as the stability of an ecosystem. In the field of disaster science, Patton defined resilience as the ability of the system to maintain its regular function, as well as its ability to cope with alterations and adjustments under extraneous conflict (Paton and Johnston, 2001). Bruneau believes that resilience focuses on the ability to quickly recover from disasters (Bruneau et al., 2003). Cai believes that resilience is the unique capability of an ecosystem to preserve the operation of its main functions when the system resolves a crisis (Cai et al., 2012). The United Nations defines resilience as the ability of a city to confront risks to withstand, maintain and cope with the influence of risk in a prompt and effective way. Climate change poses significant risks and uncertainties for the sustainable development of grain yield in all countries. As one of the important means of identifying grain yield status in grain yield areas, the assessment of grain yield has attracted the attention of many scholars, and there has been a lot of theoretical, methodological and applied research on this topic (Xu et al., 2013; Ji et al., 2014; Davenport et al., 2015). To quantitatively assess the resilience of future grain yield in China, this paper innovatively proposes climate change resilience indicators for grain yield. "Resilience" is used to indicate the ability of grain-producing areas to recover from the disturbance of climate change and realize stable grain production. Due to the obvious spatial heterogeneity of grain yield affected by climate change, the resilience of grain yield varies across different climate regions. Considering the negative correlation between grain yield resilience and climate change risk, previous research has developed grain yield risk indicators under global climate change (Wu et al., 2017; Sun et al., 2020). The article quantitatively assesses the sustainability of grain yield using two factors: grain output per person and grain unit yield.

We define the resilience index of grain output per person as the ratio of the 2000–2018 average annual grain output per person to the 2020–2050 grain output per person under global climate change. The index is calculated as follows: H m = 1 n ∑ j F j m , F j m = G j m / P j m , j = 2000 − 2018 , n = 19 (10) F w m ′ = G w m / P w m , w = 2020 − 2050 (11) S w m = H m / F w m ′ (12) where G _jm is the grain yield of province m in year j, P _jm is the population of province m in year j, F _jm is the grain output per person of p mrovince in year j, and H _m is the average annual grain output of province m. F’ _wm is the grain output per person of province m in year w, G _wm is the grain yield of province m in year w, P _wm is the population of province m in year w, and S _wm is the grain output per person resilience index of province m in year w.

We define the resilience index of grain unit yield as the ratio of the 2000–2018 average grain unit yield to the 2020–2050 grain unit yield under global climate change. The index is calculated as follows: N m = 1 n ∑ j M j m , M j m = G j m / A j m , j = 2000 − 2018 , n = 19 (13) M w m ′ = G w m / A w m , w = 2020 − 2050 (14) R w m = N m / M w m ′ (15) where G _jm is the grain yield of province m in year j, A _jm is the grain sown area of province m in year j, M _jm is the grain unit yield of province m in year j, and N _m is the average annual grain unit yield of province m. M′ _wm is the grain unit yield of province m in year w, G _yi is the grain yield of province m in year w, A _wm is the planting area of province m in year w, and R _wm is the grain unit yield resilience index of province m in year w.

Based on socioeconomic and climatic data, this paper forecasts the grain yield for SSP1–2.6, SSP2–4.5, and SSP5–8.5. Figure 2 is the flow chart of China’s grain yield resilience under global climate change.

This article analyses the resilience of grain output per person from the perspective of the grain supply balance, i.e., grain yield and grain consumption. Considering the negative correlation between grain yield resilience and climate change risk, previous research has developed grain yield risk indexes (Wu et al., 2017; Sun et al., 2020). This paper defines the resilience index of grain output per person as the ratio of the 2000–2018 average grain output per person to the 2020–2050 grain output per person under different scenarios. According to the principle of dividing by the grain production risk index, the risk index for grain output per person is as follows: if the risk index value is greater than 1.20, there is no risk, and a value under 1.20 is risky. If the risk index value is higher than 1.00, there is a low risk to grain production; if the value is below 0.90, there is high risk; and values between the two indicate moderate risk (Sun et al., 2020). Therefore, "resilience" here is inversely proportional to the risk indicator, i.e., the reciprocal of the risk indicator. In this paper, we define a resilience index of grain yield per capita value of less than 0.83 as a high resilience; other values indicate fragile conditions. A value of less than one indicates moderate resilience; a value of higher than 1.11 indicates poor resilience; and values between them indicate low resilience.

In Figure 5, the regions with poor resilience for grain output per person in China are mainly concentrated in the regions of the subtropical monsoon climate and the tropical monsoon climate. At the same time, under SSP1–2.6, SSP2-4.5, and SSP5-8.5, with increased emissions, the resilience of grain output per person in the regions of the subtropical monsoon climate and the tropical monsoon climate worsens. The gradual increase in temperature exceeds the appropriate temperature required for grain yield, resulting in a reduction in grain yield. Therefore, in subtropical monsoon climate regions and tropical monsoon climate regions, energy-saving measures and emission reductions are important measures for maintaining and stabilizing local grain yield. The grain output per person depends on the grain demand and grain supply. The eastern coastal provinces are located in the tropical monsoon climate region. The cities include Jiangsu, Zhejiang, Fujian, and Guangdong. Due to reform and opening-up policies, the spatial pattern of China's grain-producing areas has changed significantly, the regional differences between grain yield and population flow have gradually increased, and the population has gradually flowed to areas with higher urbanization levels. In the eastern coastal provinces, industrialization, urbanization and internationalization activities are much higher than those in other provinces and cities, which lead to a poor resilience. However, the resilience of grain output per person in the temperate monsoon climate region remained high under the different scenarios; at the same time, with the increase in the emissions intensity under SSP1–2.6, SSP2–4.5, and SSP5–8.5, the resilience of grain output per person improved. The resilience of grain output per person in temperate continental climate regions and plateau alpine climate regions did not change significantly under SSP1–2.6, SSP2–4.5, and SSP5–8.5, and the resilience of grain output per person was maintained in these regions.

Grain security is not only based on grain supply and demand but rather grain yield also reflects the land production capacity, and the agricultural production level is also an important factor in determining grain security. Considering the negative correlation between grain yield resilience and climate change risk, previous research has developed grain yield risk indexes (Wu et al., 2017; Sun et al., 2020). This paper defines the resilience index of the grain unit yield as the ratio of the 2000–2018 average grain unit yield to the 2020–2050 grain unit yield under SSP1–2.6, SSP2–4.5, and SSP5–8.5. According to the principle of dividing by the grain production risk index, the risk index for grain unit yield is as follows: if the risk index value is greater than 1.75, there is no risk; however, a value under 1.75 is risky. If the risk index value is higher than 1.50, there is a low risk to grain production; if the value is below 1.25, there is high risk; and values between the two indicate moderate risk (Sun et al., 2020). Therefore, "resilience" here is inversely proportional to the risk indicator, i.e., the reciprocal of the risk indicator. In this paper, we define a resilience index of grain unit yield value of less than 0.57 as a high resilience; other values indicate fragile conditions. A value of less than 0.67 indicates a moderate resilience; a value of higher than 0.8 indicates a poor resilience; and values in between indicate a low resilience.

Figure 6 shows that, in general, the resilience of grain unit yield in China's temperate monsoon climate regions is higher than that in the other climate zones in China. This is mainly due to the rapid growth of grain yield in temperate monsoon regions; on one hand, it is due to the development and upgrading of the industry; on the other hand, it benefits from global warming, which has led to an increase in the multiple cropping index and an increase in land production capacity. The increase in grain unit yield in the high-latitude area (Heilongjiang) is obvious. Due to the quick industrialization of subtropical and tropical monsoon climate regions, the loss of cultivated land in these regions is severe. In the process of balancing occupied and available cultivated land, the difference in the quality of occupied and available field land causes a decline in terms of grain production, resulting in a lower grain yield sustainability. In addition, the average temperatures in the subtropical and tropical monsoon climate zones are higher during the grain growth period than that in other climate zones. In the high-emissions scenario, SSP5–8.5, the higher the emissions intensity is, the greater the temperature increase. This increase would exceed the appropriate temperature for crop growth, weakening the land production capacity, which is also a key factor causing the low resilience of grain unit yield. In the SSP5–8.5 scenario, the resilience of grain unit yield in the monsoon region will gradually decrease from north to south. In the future, the resilience of grain unit yield in temperate continental climate regions and plateau alpine climate regions will not change substantially. Because the temperate continental climate zone is located inland, rainwater is not sufficient for crop growth; the increase in water discharge leads to temperature increases that further aggravate water shortages (Liu and Lin, 2007). The results show that the increase in temperature in the plateau alpine climate region aggravates the melting of ice and sown, increases the area of ice lakes, and increases the risk of glacier debris flows and ice lake breaks (Cui et al., 2014). However, the change in grain yield in these areas mainly depends on improving production efficiency, increasing investments in science and technology, strengthening farmland infrastructure conditions, increasing the promotion of water-saving irrigation technology, and adjusting the regional layout and variety structure of grains to a achieve sustainable and stable grain unit yield.