Multivariate Grey Prediction Model Application in Civil Aviation Carbon Emission Based on Fractional Order Accumulation and Background Value Optimization

Abstract

The GM(1,N) model, as a classical multivariate grey prediction model, can make a holistic and dynamic analysis of multiple factors and reflect the dynamic change relationship between the variable series and the related factor series. However, numerous works in the literature show that the GM(1,N) model has mechanistic defects, parametric defects, and structural defects. Therefore, the thesis establishes the OGM(1,N) model based on the GM(1,N) model by adding the linear correction term and the amount of grey action. According to the principle of dynamic optimization, the PSO algorithm is used to determine the background value. On this basis, the fractional order idea is introduced to push the model order from the integer field to the real field, and the FOBGM(1,N) model is established to systematically reduce the model error. Second, the literature in the ScienceDirect database for the last ten years is reviewed, and the carbon emission impact factors of civil aviation are selected. The calculated carbon emission values are taken as sample data based on Method 2 of Civil Aviation in Volume 2 of the 2006 IPCC Guide to National Greenhouse Gas Inventories. The results show that the prediction accuracy of the model has an increasing trend after multi-layer and multi-angle optimization. Among them, the MAPE of the OGM model and FOBGM model decreased by 24.40% and 31.86% compared with the GM(1,N) model. The 5-year average prediction accuracy of the FOBGM model reaches 99.996%, which verifies the effectiveness and practicality of the model improvement and has certain practical significance and application prospects.

1. Introduction

The grey prediction method is one of the effective methods for the analysis and prediction of systems with “small samples, poor information, and uncertainty”, and its biggest advantage is that the amount of information required for prediction is relatively small, which is more suitable in the case of difficult data acquisition and high prediction accuracy. In grey prediction theory, the most basic prediction model is the GM(1,1) model, which is based on a small amount of information and has been widely used in industry, agriculture, medicine, and other fields [1]. However, the GM(1,1) model only analyzes and predicts the change pattern of a single variable, ignoring the influence of changes in other influencing factors on the subject of study. Meanwhile, the GM(1,N) model is the basic model of the multivariate grey system modeling method and is compared with the GM(1,1) grey model. This model makes up for the defect of using historical data directly to build a time series prediction model. It can make a holistic and dynamic analysis of multiple factors and reflect the dynamic change relationship between the study variable series and the related factor series.

At present, many scholars have applied multivariate grey prediction models for forecasting in different fields and optimized the GM(1,N) model from different perspectives to improve its prediction accuracy. Jiang S. et al. (2017) combined grey GM(1,N), Markov theory, and metabolic ideas and empirically showed that the relative error of the metabolism GM(1,N)–Markov model was 32.73% smaller than that of the GM(1,N) prediction model [2].

Wu and Zhang (2018) [3] proposed a new GMC(1,N) model for information priority accumulation, which adjusted the weight of the data by adding a parameter and ultimately improved the accuracy of the prediction model.

Xiong et al. (2018) [4] proposed a nonlinear multivariate NGM(1,N) model based on kernel function and grey radius, which was established based on the kernel number and grey radius sequence of an interval grey number sequence. The kernel number and grey radius of the interval grey number sequence were simulated and predicted. According to the formula of kernel function and grey radius, the upper and lower bounds of the interval gray number were derived so as to simulate and predict the model. The NGM(1,N) model based on the interval grey number was used to forecast AQI under haze weather. Experimental results showed that the model had high prediction accuracy.

Wu et al. (2019) [5] conducted a systematic study of the GM(α, n) model using grey modeling techniques and the forward difference method to calculate the simulated and predicted values by transforming fractional order differential equations into fractional order difference equations and proposed a stochastic testing scheme to verify the accuracy of the new GM(α, n) model. The results showed that the GM(α, n) model had a high potential for energy consumption forecasting in China.

Cheng et al. (2020) [6] derived the modified model of the conventional grey differential equation based on the GM(1,N) whitening equation. Using the grey differential equation estimation parameters of the improved GM(1,N) model, the parameter estimation methods under three different background values were given. Compared with the traditional model, the simulation accuracy and prediction accuracy of the improved GM(1,N) model were significantly improved.

According to Zeng et al. (2020) [7], the extreme value of the independent variable was one of the important factors affecting the simulation and prediction results of the GM(1,N) dependent variable, a new multivariate grey prediction model was constructed on the basis of the smooth generation of variable weight independent variable sequences. The performance of NMGM(1,N) was verified with an example

Shen et al. (2021) [8] proposed an optimized discrete GMC(1,N) model (called ODGMC) to further improve the prediction accuracy and stability of GMC(1,N). In particular, a linear correction term was introduced into the new model, and the time response function of the new model was derived. The ODGMC(1,N) model proposed in this paper not only adjusted the relationship between the dependent and independent variables but also exhibited better stability than GMC(1,N) and its discrete form. The algorithm results showed that the proposed ODGMC(1,N) model had better fitting and prediction accuracy than the traditional GM(1,N) model, GMC(1,N) model, and its discrete form regardless of whether the dependent variable series was increasing, decreasing or fluctuating.

By analyzing the relevant literature, it was found that scholars from different countries have improved and optimized the GM(1,N) model in different research fields and obtained relatively better prediction accuracy compared with the traditional GM(1,N) model. Some studies [5,9] combined the kernel function with GM(1,N) through the combination of embedded models, or they cleverly combined the kernel function with GM(1,N). The multi-variable grey prediction model with higher prediction accuracy was established to reduce the prediction error of the model to a certain extent, but the inherent defects of the GM(1,N) model and the optimization of background values were not studied. The authors [2,3,4,7] studied the optimization of background values of GM(1,N) and GM(1,N) models based on the idea of metabolism and the combination of kernel function and GM(1,N), respectively, to achieve the purpose of improving prediction accuracy. However, these studies did not deeply study the mechanism of GM(1,N) model construction and ignored the inherent defects of the GM(1,N) model, and the prediction accuracy of the model was affected. Aiming at the inherent defects of the GM(1,N) model, some studies [7,8] improved the GM(1,N) model from multiple perspectives and achieved good results. However, these studies lacked the optimization of the cumulative order of GM(1,N) and failed to further improve the prediction accuracy of the multivariate grey model.

This paper provided an in-depth analysis of the inherent flaws of the GM(1,N) model and found that there were still some flaws in the construction process of the GM(1,N) model. GM(1,N) solved the whitening equation and derived the approximate time corresponding function in a relatively idealized manner, resulting in its mechanistic flaws [8]. GM(1,N) treated the parameter column estimates as model parameters of the model approximate time response equation, which was its parametric defects [10]. GM(1,N) did not mine enough grey action from itself, lacked the study of the effect of the GM(1,N) model performance by the linear relationship of the number of terms $k$, and when N = 1, the model could not achieve structural equivalence with the GM(1,1) model, indicating its structural defects [8].

In all, multi-stage and multi-angle improvement and optimization of the multi-variable grey model were carried out, including background value optimization, fractional optimization, linear increment optimization, and so on. Specifically speaking, the paper added a linear correction term and grey action quantity to the GM(1,N) model to compensate for its inherent defects and built the OGM(1,N) model. Then the Particle Swarm Optimization (PSO) algorithm was used to find the optimal background value, and the OBGM(1,N) model was constructed. The OBGM(1,N) model was improved and optimized using the PSO algorithm and fractional order algorithm at multiple levels and angles, which was called the FOBGM(1,N) model. Then the literature in the ScienceDirect database for the last ten years was reviewed, and the carbon emission impact factors of civil aviation were selected by research frequency statistics and correlation analysis. The carbon emission calculation method for civil aviation of the 2006 IPCC Guidelines for National Greenhouse Gas Inventories was introduced, relevant data were collected, and the prediction accuracy of the model before and after optimization was compared and analyzed to prove the effectiveness of the model optimization, so as to propose a carbon emission prediction model with high prediction accuracy for civil aviation management to evaluate the effectiveness of implementing current emission reduction measures. It should be noted that the GM(1,N) model based on fractional order accumulation and background value optimization developed in this paper can also be applied to the carbon emission forecasting of other transportation systems according to practical needs.

2. Basic Theory

2.1. Analysis of the Correlation Characteristics of Traditional Multivariate Grey Prediction Model

In terms of the number of input variables to the model, the grey prediction model is divided into a multi-variable grey prediction model and a univariable grey prediction model. Among them, the GM(1,1) model (i.e., a grey prediction model with a first-order equation and only one variable) is a typical representative of univariate grey prediction models, which take a single time-series data as the modeling object. According to the system operation law contained in the time series data, the grey generation method is used to mine the data information so as to predict the future development of the system. As the prediction model with the widest application range and the most research achievements, the univariable grey prediction model has a simple structure and is easy to understand without considering the impact of impact factors on the future development of the system. However, its fitting results often show saturated S-shaped characteristics or exponential forms and do not consider the influence of external environment changes on the development trend of the system, so it cannot predict the “inflection point” in the development process of the system. Therefore, the limitations of this model are relatively large. As a representative of a multivariate grey prediction model, the GM(1,N) model is a causal prediction method, which consists of an N-1 series of related factors and one series of system characteristics. In the process of model construction, the influence of external environment changes on the development trend of the system is included in the research scope. The multiple regression algorithm is similar to this model; however, the former is based on probability statistics, and the latter is based on grey theory, and there are essential differences between them. Compared to univariate grey prediction models, GM(1,N) models no longer have a limited single simulation capability in terms of structure. For a long time, the GM(1,N) model has been used more for systematic analysis, and its prediction ability has not been widely utilized. This is because the GM(1,N) model has certain defects in the model construction mechanism and structure, which results in its prediction accuracy is often lower than the GM(1,1) model in the process of use.

2.2. GM(1,N)

2.2.1. Model Definition

Let $X_1^{\left(0\right)}$ be the sequence of dependent variables (i.e., the sequence of system characteristics), $X_1^{\left(0\right)}=\left(x_1^{\left(0\right)}\left(1\right),x_1^{\left(0\right)}\left(2\right),\cdots ,x_1^{\left(0\right)}\left(m\right)\right)$. Let $X_i^{\left(0\right)}\left(i=2,3,\cdots N\right)$ be the sequence of independent variables with a high correlation with $X_1^{\left(0\right)}$ (i.e., the sequence of explanatory variables), $X_i^{\left(0\right)}=\left(x_i^{\left(0\right)}\left(1\right),x_i^{\left(0\right)}\left(2\right),\cdots ,x_i^{\left(0\right)}\left(m\right)\right)$. Further, $X_j^{\left(1\right)}$ is the 1-AGO sequence of $X_j^{\left(0\right)}$$\left(j=1,2,\cdots ,N\right)$, $X_j^{\left(1\right)}=\left(x_j^{\left(1\right)}\left(1\right),x_j^{\left(1\right)}\left(2\right),\cdots ,x_j^{\left(1\right)}\left(m\right)\right)$. Where $x_j^{\left(1\right)}\left(k\right)={\displaystyle {\sum }_{g=1}^k{x}_j^{\left(0\right)}\left(g\right)},k=1,2,\cdots ,m.$ $Z_1^{\left(1\right)}$ is the immediate mean generating sequence of $X_1^{\left(1\right)}$, which can be calculated by

$$Z_1^{\left(1\right)}=\left(z_1^{\left(1\right)}\left(2\right),z_1^{\left(1\right)}\left(3\right),\cdots ,z_1^{\left(1\right)}\left(m\right)\right).$$

(1)

In Equation (1) $z_1^{\left(1\right)}\left(k\right)=0.5\times \left(x_1^{\left(1\right)}\left(k\right)+x_1^{\left(1\right)}\left(k-1\right)\right),k=2,3,\cdots ,m.$ The GM(1,N) model is expressed as

$$x_1^{\left(0\right)}\left(k\right)+a{z}_1^{\left(1\right)}\left(k\right)={\displaystyle \sum _{i=2}^N{b}_i{x}_i^{\left(1\right)}\left(k\right)}.$$

(2)

2.2.2. Parameter Estimation and Time Response Equation of GM(1,N) Model

In the GM(1,N) model, $-a$ is called the system development coefficient, $b_i{x}_i^{\left(1\right)}\left(k\right)$ is called the driver term, $b_i$ is called the driver coefficient, and $\widehat{a}={\left[a,b_1,b_2,\cdots ,b_N\right]}^T$ is called the parameter column for which the least squares estimate satisfies using

$$P={\left(B^{T}B\right)}^{-1}{B}^{T}Y.$$

(3)

Among them, B and Y can be set

$$B=\left[\begin{array}{cccc}-z_1^{\left(1\right)}\left(2\right)& x_2^{\left(1\right)}\left(2\right)& \cdots & x_N^{\left(1\right)}\left(2\right)\\ -z_1^{\left(1\right)}\left(3\right)& x_2^{\left(1\right)}\left(3\right)& \cdots & x_N^{\left(1\right)}\left(3\right)\\ ⋮& ⋮& & ⋮\\ -z_1^{\left(1\right)}\left(m\right)& x_2^{\left(1\right)}\left(m\right)& \cdots & x_N^{\left(1\right)}\left(m\right)\end{array}\right],Y=\left[\begin{array}{c}{x}_1^{\left(0\right)}\left(2\right)\\ x_1^{\left(0\right)}\left(3\right)\\ ⋮\\ x_1^{\left(0\right)}\left(m\right)\end{array}\right].$$

(4)

Then we can obtain

$$\frac{d{x}_1^{\left(1\right)}}{dt}+a{x}_1^{\left(1\right)}={\displaystyle \sum _{i=2}^N{b}_i{x}_i^{\left(1\right)}}.$$

(5)

Equation (5) is the whitening equation of the GM(1,N) model, also called the shadow equation.

Let the sequence $X_i^{\left(0\right)}$, $X_i^{\left(1\right)}\left(i=1,2,\cdots ,N\right)$, $Z_1^{\left(1\right)}$ and the matrix $B,Y,\widehat{a}$ be as described above, then the solution of the whitening equation $\frac{d{x}_1^{\left(1\right)}}{dt}+a{x}_1^{\left(1\right)}={\displaystyle \sum _{i=2}^N{b}_i{x}_i^{\left(1\right)}}$ is

$$x_1^{\left(1\right)}\left(t\right)=e^{-at}\left[x_1^{\left(0\right)}\left(1\right)-t{\displaystyle \sum _{i=2}^N{b}_i{x}_i^{\left(1\right)}\left(0\right)}+{\displaystyle \sum _{i=2}^N{\displaystyle \int b_i{x}_i^{\left(1\right)}\left(t\right)e^{at}dt}}\right].$$

(6)

When the variation in $X_i^{\left(1\right)}\left(i=1,2,\cdots ,N\right)$ is very small, visualizing ${\displaystyle {\sum }_{i=2}^N{b}_i{x}_i^{\left(1\right)}\left(k\right)}$ as a grey constant, the approximate time response equation of the GM(1,N) model is

$$\begin{array}{c}{x}_1^{\left(1\right)}\left(k+1\right)=\left[x_1^{\left(0\right)}\left(1\right)-\frac{1}{a}{\displaystyle \sum _{i=2}^N{b}_i{x}_i^{\left(1\right)}\left(k+1\right)}\right]e^{-ak}\\ +\frac{1}{a}{\displaystyle \sum _{i=2}^N{b}_i{x}_i^{\left(1\right)}\left(k+1\right)}\end{array}.$$

(7)

Then the cumulative reduction equation is

$${\widehat{x}}_1^{\left(0\right)}\left(k+1\right)={\alpha }^{\left(1\right)}{\widehat{x}}_1^{\left(1\right)}\left(k+1\right)={\widehat{x}}_1^{\left(1\right)}\left(k+1\right)-{\widehat{x}}_1^{\left(1\right)}\left(k\right).$$

(8)

2.2.3. Defects of the GM(1,N) Model

In the construction of the GM(1,N) model, there are still some shortcomings that affect the prediction accuracy of the GM(1,N) model.

(1) The GM(1,N) model solves the whitening equation and derives the approximate time equivalent function in a relatively idealized way. In reality, the “small magnitude of variation” of $X_i^{\left(1\right)}$ is difficult to satisfy. This is because the variables represented by $X_i^{\left(1\right)}\left(i=1,2,\cdots ,N\right)$ are different, so their development trends, dynamic rules, and change characteristics are generally different, and the amplitude of change is difficult to guarantee, which results in the unstable prediction performance of the GM(1,N) model, resulting in the mechanism defects of the GM(1,N) model [10].

(2) GM(1,N) is not a time-response equation derived using $x_1^{\left(0\right)}\left(k\right)+a{z}_1^{\left(1\right)}\left(k\right)={\displaystyle \sum _{i=2}^N{b}_i{x}_i^{\left(1\right)}\left(k\right)}$, but is computed using its shadow equation. Therefore, the parameter column $\widehat{a}={\left[a,b_1,b_2,\cdots ,b_N\right]}^T$ calculates its estimates with $x_1^{\left(0\right)}\left(k\right)+a{z}_1^{\left(1\right)}\left(k\right)={\displaystyle \sum _{i=2}^N{b}_i{x}_i^{\left(1\right)}\left(k\right)}$. The GM(1,N) model, however, regards the estimated value of the parameter column as the model parameter of the approximate time response of the model. The “dislocation” of the GM(1,N) model parameter estimation and application object causes its unstable prediction performance, which is its parameter defect [4].

(3) The GM(1,N) model has a simple structure and is a state model and factor model. The model does not have a sufficient amount of grey action mined from itself and lacks a study on the effect of the GM(1,N) model performance by the linear relationship between the number of terms $k$. In addition, the GM(1,N) model serves as a first-order grey prediction model for $N$ variables, but when N = 1, the model cannot achieve structural equivalence with the GM(1,1) model. This indicates that the GM(1,N) model also has structural flaws, resulting in relatively low prediction accuracy [10].

2.3. Prediction Accuracy Evaluation System

Mean absolute percentage error (MAPE) can be used to compare the prediction accuracy of different models [11], which can be obtained with

$$MAPE=\frac{{\displaystyle \sum _{i=1}^n{q}_i}}{n}\times 100.$$

(9)

where absolute error $e_i$ and relative error $q_i$ are calculated as

$$e_i=\left|x_i-{\widehat{x}}_i\right|,q_i=\frac{e_i}{x_i}.$$

(10)

Further, the mean square error (MSE), root mean square error (RMSE), and mean absolute error (MAE) are calculated as

$$MSE=\frac{{\displaystyle \sum _{i=1}^n{\left(x_i-{\widehat{x}}_i\right)}^2}}{n},$$

(11)

$$RMSE=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n{\left(x_i-{\widehat{x}}_i\right)}^2}},$$

(12)

$$MAE=\frac{1}{n}{\displaystyle \sum _{i=1}^n\left|x_i-{\widehat{x}}_i\right|}.$$

(13)

2.4. Carbon Emission Measurement of Civil Aviation

In order to verify the prediction performance of the model, relevant data from the China civil aviation carbon emission prediction study were used in this paper to verify the effectiveness of the model before and after improvement. The World Meteorological Organization (WMO) and the United Nations Environment Programme (UNEP) jointly established the Intergovernmental Panel on Climate Change (IPCC). This paper calculates carbon emissions based on Method 2 of Civil Aviation in Volume 2 of the 2006 IPCC Guidelines for National Greenhouse Gas Inventories (hereinafter referred to as the Guidelines). The flight of an aircraft has two stages: the LTO stage and the cruise stage. The method separates aviation carbon emissions above and below 914 m (3000 feet) during flight, that is, the LTO phase and cruise phase carbon emissions during flight. The specific calculation process is as

$$E=E_{LTO}+E_X$$

(14)

where $E$ is the carbon emissions from air transport.$E_{LTO}$ indicates the carbon emissions from air transport in the LTO phase and $E_X$ indicates the carbon emissions from air transport in the cruise phase. $E_{LTO}$ and $E_X$ are calculated as

$$E_{LTO}=N_{LTO}\cdot F_{LTO},$$

(15)

$$K_{LTO}=N_{LTO}\cdot U_{LTO},$$

(16)

$$E_X=\left(K-K_{LTO}\right)\cdot F_X,$$

(17)

$$F_X=\frac{H}{(1\times {10}^9)/L}.$$

(18)

$N_{LTO}$ indicates the number of aircraft landings and takeoffs in the national aviation industry; $F_{LTO}$ indicates the CO₂ emission factor in the LTO phase, using the average value of each aircraft type in the Guide, i.e., 4341 kg/LTO; $K_{LTO}$ indicates the fuel consumption in the LTO phase; $U_{LTO}$ indicates the fuel consumption per LTO, using the average value of each aircraft type in the Guide, i.e., 1374 kg/LTO; $K$ indicates the total fuel consumption; $F_X$ indicates the CO₂ emission factor in the cruise phase H; indicates the amount of CO₂ per unit calorific value, using the recommended value of 71,500 kg/TJ in the Guide; and L indicates the low-level heating value of aviation fuel, using the recommended value of 44,100 KJ/kg in the Guide for Accounting Methods and Reporting of Greenhouse Gas Emissions of Chinese Civil Aviation Enterprises prepared by the National Development and Reform Commission.

3. Model Improvement and Optimization

3.1. Model Optimization Ideas

In order to solve the mechanistic defects, parametric defects, and structural defects of GM(1,N), this paper adds linear correction terms and grey action quantities to the GM(1,N) and establishes the OGM(1,N). For the multivariate grey model, the background value coefficient can affect its prediction accuracy to a certain extent, and the value is generally set to 0.5. In this paper, the global search capability of the PSO algorithm is used to search for its background value coefficient, and then the order of OBGGM(1,N) is extended from the integer domain to the real domain using the PSO algorithm and fractional order algorithm to establish the FOBGM model. The improved GM(1,N) model is optimized by multi-level and multi-angle, aiming to further improve the prediction accuracy. The calculation process is shown in Figure 1.

3.2. GM(1,N) Model Structure Improvement

3.2.1. Definition of the OGM(1,N) Model

Let the OGM(1,N) model (optimizing grey model) be

$$x_1^{\left(0\right)}\left(k\right)+a{z}_1^{\left(1\right)}\left(k\right)={\displaystyle \sum _{i=2}^N{b}_i{x}_i^{\left(1\right)}\left(k\right)}+h_1\left(k-1\right)+h_2$$

(19)

where $h_1\left(k-1\right)$ and $h_2$ are the linear correction term and the grey action of the model.

3.2.2. Parameter Estimation and Time Corresponding Equation of OGM(1,N)

The sequences $X_1^{\left(0\right)}$, $Z_1^{\left(1\right)}$, and $X_i^{\left(1\right)}\left(i=1,2,\cdots ,N\right)$ are as described in Section 2.2.1, then the least squares estimates of the OGM(1,N) model parameter columns $\widehat{p}={\left[b_2,b_3,\cdots ,b_N,a,h_1,h_2\right]}^T.$ satisfy

$$\widehat{p}={\left(B^{T}B\right)}^{-1}{B}^{T}Y.$$

(20)

Among them, B and Y can be set

$$\begin{array}{c}B=\left[\begin{array}{ccccccc}{x}_2^{\left(1\right)}\left(2\right)& x_3^{\left(1\right)}\left(2\right)& \cdots & x_N^{\left(1\right)}\left(2\right)& -z_1^{\left(1\right)}\left(2\right)& 1& 1\\ x_2^{\left(1\right)}\left(3\right)& x_3^{\left(1\right)}\left(3\right)& \cdots & x_N^{\left(1\right)}\left(3\right)& -z_1^{\left(1\right)}\left(3\right)& 2& 1\\ ⋮& ⋮& & ⋮& ⋮& ⋮& ⋮\\ x_2^{\left(1\right)}\left(m\right)& x_3^{\left(1\right)}\left(m\right)& \cdots & x_N^{\left(1\right)}\left(m\right)& -z_1^{\left(1\right)}\left(m\right)& m-1& 1\end{array}\right]\\ Y=\left[\begin{array}{c}{x}_1^{\left(0\right)}\left(2\right)\\ x_1^{\left(0\right)}\left(3\right)\\ ⋮\\ x_1^{\left(0\right)}\left(m\right)\end{array}\right]\end{array}.$$

(21)

Then the difference model of OGM(1,N) can be described with

$${\widehat{x}}_1^{\left(0\right)}\left(k\right)={\displaystyle \sum _{i=2}^N{b}_i{\widehat{x}}_i^{\left(1\right)}\left(k\right)}-a{z}_1^{\left(1\right)}\left(k\right)+h_1\left(k-1\right)+h_2.$$

(22)

The time-corresponding equation is

$$\begin{array}{c}{\widehat{x}}_1^{\left(1\right)}\left(k\right)={\displaystyle \sum _{t=1}^{k-1}\left[{\mu }_1{\displaystyle \sum _{i=2}^N{\mu }_2^{t-1}{b}_i{x}_i^{\left(1\right)}\left(k-t+1\right)}\right]}+{\mu }_2^{k-1}{\widehat{x}}_1^{\left(1\right)}\left(1\right)\\ +{\displaystyle \sum _{j=0}^{k-2}{\mu }_2^j\left[\left(k-j\right){\mu }_3+{\mu }_4\right],k=2,3,\cdots }\end{array}.$$

(23)

Among them, ${\mu }_1=\frac{1}{1+0.5a},{\mu }_2=\frac{1-0.5a}{1+0.5a},{\mu }_3=\frac{h_1}{1+0.5a},{\mu }_4=\frac{h_2-h_1}{1+0.5a}$.

The cumulative reduction equation of OGM(1,N) is

$${\widehat{x}}_1^{\left(0\right)}\left(k\right)={\widehat{x}}_1^{\left(1\right)}\left(k\right)-{\widehat{x}}_1^{\left(1\right)}\left(k-1\right),k=2,3,\cdots .$$

(24)

OGM(1,N) does not assume that “${\displaystyle {\sum }_{i=2}^N{b}_i{x}_i^{\left(1\right)}\left(k\right)}$ is regarded as a grey quantity when the variation of $X_i^{\left(1\right)}\left(i=1,2,\cdots ,N\right)$ is small”. There is no “misalignment” between the shadow equation $d{x}_1^{\left(1\right)}/dt+a{x}_1^{\left(1\right)}={\displaystyle {\sum }_{i=2}^N{b}_i{x}_i^{\left(1\right)}}$ and the parameter column $\widehat{p}={\left[b_1,b_2,\cdots ,b_N,a,h_1,h_2\right]}^T$. Furthermore, the linear correction term and the grey action quantity are added to the model to make the model structure more reasonable. OGM(1,N) is an optimization improvement of GM(1,N), which better compensates for the mechanism defects, parameter defects, and structural defects of the latter.

3.3. Optimization of Background Values for OGM(1,N)

The background value coefficient is a parameter in grey theory that can have an impact on the prediction accuracy of multivariate grey prediction models to a certain extent. However, multivariate grey prediction models generally use 0.5 as their background value coefficients with the aim of simplifying the model construction process. Therefore, this section improves the OGM(1,N) model by optimizing its background value coefficients with the help of the global search capability of the PSO algorithm, calls the improved model OBGM(1,N), and derives the time response function and the final reduced equation of OBGM(1,N).

3.3.1. OBGM(1,N)

OBGM(1,N) is defined by Equation (25), which is an OGM(1,N) model with a background value coefficient of $\xi \left(0<\xi <1\right)$.

$$\begin{array}{l}{x}_1^{\left(0\right)}\left(k\right)+a\xi x_1^{\left(1\right)}\left(k\right)+a\left(1-\xi \right)x_1^{\left(1\right)}\left(k-1\right)={\displaystyle \sum _{i=2}^N{b}_i{x}_i^{\left(1\right)}\left(k\right)}+kc+d,\\ k=2,3,\cdots ,m.\end{array}$$

(25)

The parameter column $\widehat{p}={\left[b_2,b_3,\cdots ,b_N,a,c,d\right]}^T$ is computed using the least squares method. Provided that $\widehat{p}={\left[b_2,b_3,\cdots ,b_N,a,c,d\right]}^T$ is known, If a certain $\xi$ minimizes the mean absolute percentage error of sequence $X_1^{\left(0\right)}$, it is the optimal background value coefficient of the OGM(1,N) model. Further, this model is called the OBGM(1,N) model.

The least squares estimate of the parameter column $\widehat{p}={\left[b_2,b_3,\cdots ,b_N,a,c,d\right]}^T$ satisfies

$$\widehat{p}={\left(E^{T}E\right)}^{-1}{E}^{T}S.$$

(26)

Among them, E and S can be set

$$E=\left[\begin{array}{ccccccc}{x}_2^{\left(1\right)}\left(2\right)& x_3^{\left(1\right)}\left(2\right)& \cdots & x_N^{\left(1\right)}\left(2\right)& \left(\xi -1\right)x_1^{\left(1\right)}\left(1\right)-\xi x_1^{\left(1\right)}\left(2\right)& 2& 1\\ x_2^{\left(1\right)}\left(3\right)& x_3^{\left(1\right)}\left(3\right)& \cdots & x_N^{\left(1\right)}\left(3\right)& \left(\xi -1\right)x_1^{\left(1\right)}\left(2\right)-\xi x_1^{\left(1\right)}\left(3\right)& 3& 1\\ ⋮& ⋮& & ⋮& ⋮& ⋮& ⋮\\ x_2^{\left(1\right)}\left(m\right)& x_3^{\left(1\right)}\left(m\right)& \cdots & x_N^{\left(1\right)}\left(m\right)& \left(\xi -1\right)x_1^{\left(1\right)}\left(m-1\right)-\xi x_1^{\left(1\right)}\left(m\right)& m& 1\end{array}\right]$$

(27)

$$S=\left[\begin{array}{c}{x}_1^{\left(0\right)}\left(2\right)\\ x_1^{\left(0\right)}\left(3\right)\\ ⋮\\ x_1^{\left(0\right)}\left(m\right)\end{array}\right]$$

(28)

then ${\widehat{x}}_1^{\left(0\right)}\left(k\right)+a\xi {\widehat{x}}_1^{\left(1\right)}\left(k\right)+a\left(1-\xi \right){\widehat{x}}_1^{\left(1\right)}\left(k-1\right)={\displaystyle \sum _{i=2}^N{b}_i{x}_i^{\left(1\right)}\left(k\right)}+kc+d,k=2,3,\cdots ,m.$ Because ${\widehat{x}}_1^{\left(0\right)}\left(k\right)={\widehat{x}}_1^{\left(1\right)}\left(k\right)-{\widehat{x}}_1^{\left(1\right)}\left(k-1\right)$, then, we can obtain ${\widehat{x}}_1^{\left(1\right)}\left(k\right)-{\widehat{x}}_1^{\left(1\right)}\left(k-1\right)+a\xi {\widehat{x}}_1^{\left(1\right)}\left(k\right)+a\left(1-\xi \right){\widehat{x}}_1^{\left(1\right)}\left(k-1\right)={\displaystyle \sum _{i=2}^N{b}_i{x}_i^{\left(1\right)}\left(k\right)+kc+d}$. The collation is obtained from ${\widehat{x}}_1^{\left(1\right)}\left(k\right)=\frac{1}{1+a\xi }{\displaystyle \sum _{i=2}^N{b}_i}{x}_i^{\left(1\right)}\left(k\right)+\frac{1-a+a\xi }{1+a\xi }{\widehat{x}}_1^{\left(1\right)}\left(k-1\right)+\frac{c}{1+a\xi }k+\frac{d}{1+a\xi }$. Let $v_1=\frac{1}{1+a\xi },v_2=1-\frac{a}{1+a\xi },v_3=\frac{c}{1+a\xi },v_4=\frac{d}{1+a\xi }$, then ${\widehat{x}}_1^{\left(1\right)}\left(k\right)=v_1{\displaystyle \sum _{i=2}^N{b}_i{x}_i^{\left(1\right)}\left(k\right)}+v_2{\widehat{x}}_1^{\left(1\right)}\left(k-1\right)+v_{3}k+v_4,k=2,3,\cdots m.$ When $k=2$, ${\widehat{x}}_1^{\left(1\right)}\left(2\right)=v_1{\displaystyle \sum _{i=2}^N{b}_i{x}_i^{\left(1\right)}\left(2\right)}+v_2{x}_1^{\left(1\right)}\left(1\right)+2v_3+v_4$.

When $k=3$, ${\widehat{x}}_1^{\left(1\right)}\left(3\right)=v_1{\displaystyle \sum _{i=2}^N{b}_i{x}_i^{\left(1\right)}\left(3\right)}+v_2{\widehat{x}}_1^{\left(1\right)}\left(2\right)+3v_3+v_4$. Substituting ${\widehat{x}}_1^{\left(1\right)}\left(2\right)$, and the equation is reduced to ${\widehat{x}}_1^{\left(1\right)}\left(3\right)=v_1{\displaystyle \sum _{i=2}^N{b}_i{x}_i^{\left(1\right)}\left(3\right)}+v_1{v}_2{\displaystyle \sum _{i=2}^N{b}_i{x}_i^{\left(1\right)}\left(2\right)}+v_2^2{x}_1^{\left(1\right)}\left(1\right)+2v_2{v}_3+v_2{v}_4+3v_3+v_4$.

Continuing the derivation, when $k=p$ ${\widehat{x}}_1^{\left(1\right)}\left(p\right)=v_1{\displaystyle \sum _{i=2}^N{b}_i{x}_i^{\left(1\right)}\left(p\right)}+v_1{v}_2{\displaystyle \sum _{i=2}^N{b}_i{x}_i^{\left(1\right)}\left(p-1\right)}+\cdots +v_1{v}_2^{p-2}{\displaystyle \sum _{i=2}^N{b}_i{x}_i^{\left(1\right)}\left(2\right)}+v_2^{p-1}{\widehat{x}}_1^{\left(1\right)}\left(1\right)+2v_2^{p-2}{v}_3+v_2^{p-2}{v}_4+3v_2^{p-3}{v}_3+v_2^{p-3}{v}_4+\cdots +p·v_3+v_4$.

Simplify to get ${\widehat{x}}_1^{\left(1\right)}\left(p\right)={\displaystyle \sum _{u=1}^{p-1}\left[v_1{\displaystyle \sum _{i=2}^N{v}_2^{u-1}{b}_i{x}_i^{\left(1\right)}\left(p-u+1\right)}\right]}+{\displaystyle \sum _{v=0}^{p-2}{v}_2^v\left[\left(p-v\right)v_3+v_4\right]}+v_2^{p-1}{\widehat{x}}_1^{\left(1\right)}\left(1\right)$, which is the time response function of the OBGM(1,N) model, $k=p=2,3,\cdots ,m$.

From ${\widehat{x}}_1^{\left(0\right)}\left(k\right)={\widehat{x}}_1^{\left(1\right)}\left(k\right)-{\widehat{x}}_1^{\left(1\right)}\left(k-1\right)$, we can obtain the final reduced form of the OBGM(1,N) model. Thus, the time response equation is

$${\widehat{x}}_1^{\left(1\right)}\left(k\right)={\displaystyle \sum _{u=1}^{k-1}\left[v_1{\displaystyle \sum _{i=2}^N{v}_2^{u-1}{b}_i{x}_i^{\left(1\right)}\left(k-u+1\right)}\right]}+{\displaystyle \sum _{v=0}^{k-2}{v}_2^v\left[\left(k-v\right)v_3+v_4\right]}+v_2^{k-1}{\widehat{x}}_1^{\left(1\right)}\left(1\right).$$

(29)

The final reduced equation is

$$\begin{array}{r}{\widehat{x}}_1^{\left(0\right)}\left(k\right)=v_1{\displaystyle \sum _{d=2}^k\left[{\displaystyle \sum _{i=2}^N{\left(v_2-1\right)}^{\left[\frac{d-2}{d}\right]}}{v}_2^{\left[\frac{d-2}{d}\right]\left[\frac{d-3}{d}\right]\left(d-3\right)}{b}_i{x}_i^{\left(1\right)}\left(k-d+2\right)\right]}\\ +{\displaystyle \sum _{v=0}^{k-3}\left[\left(v_2^{k-2}{\left(2v_3+v_4+1\right)}^{\left[\frac{v}{k-3}\right]}-v_2^{k-2}\right)+v_2^v{v}_3\right]}\end{array}.$$

(30)

Equation (30) contains two parts. The first half refers to the effect of the independent variable $X_i^{\left(1\right)}(i=2,3,\cdots ,N)$ on the dependent variable $X_1^{\left(0\right)}$. The second half refers to the effect of the constant terms $v_1,v_2,v_3$, and $v_4$ on the dependent variable $X_1^{\left(0\right)}$.

3.3.2. Background Value Coefficient of OBGM(1,N) Model Optimized by PSO Algorithm

In 1995, Kennedy and Eberhart proposed a global optimization algorithm, namely particle swarm optimization (PSO) [12]. The PSO algorithm is simple in concept, easy to program, has fewer parameters, and is widely used in the fields of neural network training and function optimization. In this paper, the PSO optimization algorithm is used to find the optimal background value coefficients of the OBGM(1,N) model and the mean absolute percentage error (i.e., MAPE) of the output results of the OBGM(1,N) model is used as the objective function. When the minimum objective function is found, its corresponding background value is the optimal background value coefficient.

At the minimum MAPE, the optimal background value coefficient of the OBGM(1,N) model is essentially to solve the following optimization problems, which can be seen with Equation (31).

$$\min f\left(\xi \right)=\frac{1}{m-1}{\displaystyle \sum _{k=2}^n\frac{\left|{\widehat{x}}_1^{\left(0\right)}\left(k\right)-x_1^{\left(0\right)}\left(k\right)\right|}{x_1^{\left(0\right)}\left(k\right)}},\xi \in R^{+}.$$

(31)

The computational procedure for finding the optimal background value coefficients for the OBGM(1,N) model using the PSO algorithm $\xi$ is shown below.

Step 1: Set the PSO algorithm parameters to initialize the velocities and positions of the particles in the particle swarm.

Step 2: Set the current position in the particle to $pBest$ and the best particle position in the group to $gBest=1$.

Step 3: Calculate the MAPE of the OBGM(1,N) model when $\xi =pBest$.

Step 4: Update all particles as follows.

Step 4.1: Update of particle velocity and position (shown in Equations (32) and (33))

$$\begin{array}{c}V=\omega \times V+c_1\times rand\times \left(pBest-Present\right)\\ +c_2\times rand\times \left(gBest-Present\right)\end{array},$$

(32)

$$Present=Present+V.$$

(33)

Step 4.2: Update $pBest$ to the new location if the adaptation is better than the adaptation of $pBest$.

Step 4.3: Update $gBest$ to the new position if the adaptation is better than the adaptation of $gBest$.

Step 5: Determine whether the maximum number of iterations is reached or the accuracy requirements are met. If yes, proceed to Step 6; Otherwise, go back to Step 4.

Step 6: Output $gBest$ to get the optimal value of the optimal background value coefficient $\xi$, and output the simulated or predicted value of the OBGM(1,N) model of this optimal background value coefficient $\xi =gBest$, and the run is finished.

3.4. Fractional Order Optimization of OBGM(1,N)

The essence of the OBGM(1,N) model is a multi-variable grey prediction model of first-order equations. By mining the information contained in the data with the help of the first-order cumulative sequence, the fractional order idea is introduced into the OBGM(1,N) model, and the order of the model is extended from the integer domain to the real domain, in order to improve prediction accuracy. The OBGM(1,N) model is further optimized by the optimal cumulative order. Among them, the optimization of cumulative order is determined using the global optimization ability of the PSO optimization algorithm.

3.4.1. Model Definition

Assuming that the original non-negative series sequence is $X_1^{\left(0\right)}=\left(x_1^{\left(0\right)}\left(1\right),x_1^{\left(0\right)}\left(2\right),\cdots x_1^{\left(0\right)}\left(m\right)\right)$, add r-order to the original sequence $X_i^{\left(0\right)}(i=2,3,\cdots ,N)$ to obtain the sequence $X_j^{\left(r\right)}$$(j=1,2,\cdots ,N)$, $X_j^{\left(r\right)}=\left(x_j^{\left(r\right)}\left(1\right),x_j^{\left(r\right)}\left(2\right),\cdots x_j^{\left(r\right)}\left(m\right)\right)$, which can be calculated using Equation (34) [13].

$${\mathrm{x}}_j^{(r)}(\mathrm{k})={\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{k}}\frac{\Gamma \left(r+k-i\right)}{\Gamma \left(k-i+1\right)\Gamma \left(r\right)}}{\mathrm{x}}_j^{(0)}(i),\mathrm{k}=1,2,\cdots ,\mathrm{m}.$$

(34)

The OGM(1,N) model with a background value coefficient of $\xi \left(0<\xi <1\right)$ and fractional order of $r$ can be described by Equation (35), where $\widehat{p}={\left[b_2,b_3,\cdots ,b_N,a,c,d\right]}^T$ is the column of parameters estimated by the least squares method. Under the condition that $\widehat{p}={\left[b_2,b_3,\cdots ,b_N,a,c,d\right]}^T$ is known, if $\xi$ and $r$ exist to minimize the average relative simulation error percentage of sequence $X_1^{\left(0\right)}$, the OGM(1,N) model is said to have the optimal background value coefficient $\xi$ and the optimal order $r$, and the multivariate grey prediction model is denoted as the FOBGM(1,N) model.

$$\begin{array}{l}{x}_1^{\left(0\right)}\left(k\right)+a\xi x_1^{\left(r\right)}\left(k\right)+a\left(1-\xi \right)x_1^{\left(r\right)}\left(k-1\right)={\displaystyle \sum _{i=2}^N{b}_i{x}_i^{\left(r\right)}\left(k\right)}+kc+d,\\ k=2,3,\cdots ,m.\end{array}$$

(35)

3.4.2. Parameter Estimation and Time Response Equation of FOBGM(1,N)

Let $X_j^{(-\mathrm{r})}=({\mathrm{x}}_j^{(-\mathrm{r})}(1),{\mathrm{x}}_j^{(-\mathrm{r})}(2),\cdots ,{\mathrm{x}}_j^{(-\mathrm{r})}(m))$ be the r-order reduction generator of $X_j^{\left(0\right)}$, which can be gained using Equation (36).

$${\mathrm{x}}_j^{(-r)}(\mathrm{k})={\displaystyle \sum _{\mathrm{i}=0}^{\mathrm{k}-1}{\left(-1\right)}^i\frac{\Gamma \left(r+1\right)}{\Gamma \left(r-i+1\right)\Gamma \left(i+1\right)}}{\mathrm{x}}_j^{\left(0\right)}(k-i),\mathrm{k}=1,2,\cdots ,\mathrm{m}.$$

(36)

Let the least squares estimate of the parameter column $\widehat{p}={\left[b_2,b_3,\cdots ,b_N,a,c,d\right]}^T$ satisfy

$$\widehat{p}={\left(E^{T}E\right)}^{-1}{E}^{T}S,$$

(37)

$$E=\left[\begin{array}{ccccccc}{x}_2^{\left(r\right)}\left(2\right)& x_3^{\left(r\right)}\left(2\right)& \cdots & x_N^{\left(r\right)}\left(2\right)& \left(\xi -1\right)x_1^{\left(r\right)}\left(1\right)-\xi x_1^{\left(r\right)}\left(2\right)& 2& 1\\ x_2^{\left(r\right)}\left(3\right)& x_3^{\left(r\right)}\left(3\right)& \cdots & x_N^{\left(r\right)}\left(3\right)& \left(\xi -1\right)x_1^{\left(r\right)}\left(2\right)-\xi x_1^{\left(r\right)}\left(3\right)& 3& 1\\ ⋮& ⋮& & ⋮& ⋮& ⋮& ⋮\\ x_2^{\left(r\right)}\left(m\right)& x_3^{\left(r\right)}\left(m\right)& \cdots & x_N^{\left(r\right)}\left(m\right)& \left(\xi -1\right)x_1^{\left(r\right)}\left(m-1\right)-\xi x_1^{\left(r\right)}\left(m\right)& m& 1\end{array}\right],$$

(38)

$$S=\left[\begin{array}{c}{x}_1^{\left(0\right)}\left(2\right)\\ x_1^{\left(0\right)}\left(3\right)\\ ⋮\\ x_1^{\left(0\right)}\left(m\right)\end{array}\right].$$

(39)

Then the time response equation of the FOBGM(1,N) model is

$$\begin{array}{ll}{\widehat{x}}_1^{\left(r\right)}\left(k\right)& ={\displaystyle \sum _{u=1}^{k-1}\left[v_1{\displaystyle \sum _{i=2}^N{v}_2^{u-1}{b}_i{x}_i^{\left(r\right)}\left(k-u+1\right)}\right]}\\ & +{\displaystyle \sum _{v=0}^{k-2}{v}_2^v\left[\left(k-v\right)v_3+v_4\right]}+v_2^{k-1}{\widehat{x}}_1^{\left(r\right)}\left(1\right)\end{array}.$$

(40)

The final reduced equation is

$$\left\{\begin{array}{l}{\widehat{\mathrm{x}}}_1^{\left(0\right)}(\mathrm{k})={\left({\widehat{x}}_1^{\left(r\right)}\right)}^{\left(-r\right)}\left(k\right)={\displaystyle \sum _{\mathrm{i}=0}^{\mathrm{k}-1}{\left(-1\right)}^i\frac{\Gamma \left(r+1\right)}{\Gamma \left(r-i+1\right)\Gamma \left(i+1\right)}}{\widehat{\mathrm{x}}}_1^{\left(r\right)}(k-i),\\ \mathrm{k}=2,3,\cdots ,\mathrm{m}.\\ {\widehat{\mathrm{x}}}_1^{\left(0\right)}(1)={\mathrm{x}}_1^{\left(0\right)}(1)\end{array}\right.$$

(41)

4. Empirical Study

4.1. Influencing Factor Screening

Although there is relatively little in the literature on China’s civil aviation carbon emission prediction, some scholars have studied the influencing factors of civil aviation carbon emission from different perspectives. In this paper, the literature on carbon emission factors of civil aviation from the most recent ten-year period was reviewed, and some literature on carbon emission of other modes of transport was introduced to calculate the research frequency of carbon emission factors studied so as to preliminarily screen out the impact factors of carbon emission by civil aviation.

As can be seen from

Table 1. Civil aviation carbon emission influencing factors and representative authors.

Impact Factor	Literature Statistics	Representative Authors
Civil aviation fuel consumption	16	[14,15,16,17]
Total civil aviation turnover	9	[18,19,20]
Civil aviation industry-wide operating income	9	[18,19,20]
Civil aviation transport traffic (passenger/cargo volume)	8	[14,17,20,21]
Transportation intensity	5	[16,22]
Number of civil flights	4	[19,20,21,23]
GDP per capital	4	[24]
Population	4	[25]
GDP	2	[25,26]
Energy intensity	2	[22]
Flight Type	1	[23]

, among the main indicators affecting civil aviation carbon emissions in the past decade, the research frequency of several influencing factors such as civil aviation transportation traffic (passenger/cargo volume), total civil aviation turnover, civil aviation fuel consumption, civil aviation industry-wide operating income, and transport intensity are relatively high. Therefore, in this paper, civil aviation transport traffic is divided into passenger traffic and cargo volume. Finally, six influencing factors were determined to predict China’s civil aviation carbon emissions using civil aviation passenger traffic, civil aviation cargo volume, total civil aviation turnover, civil aviation fuel consumption, civil aviation industry-wide operating income, and civil aviation transportation intensity.

4.2. Original Data Collection and Correlation Analysis

According to the above civil aviation carbon emission accounting standards and selected impact factors, the data on civil aviation fuel consumption, passenger traffic, cargo volume, total civil aviation turnover, and civil aviation industry-wide operating income were collected from 2003 to 2017. Among them, the data on civil aviation fuel consumption was from the Statistical Data on Civil Aviation of China. The data on total civil aviation turnover and the data on aircraft takeoffs and landings were from the Statistical Data on Civil Aviation of China and the China Traffic Yearbook. The data of civil aviation industry-wide operating income were from the Statistical Bulletin on the Development of the Civil Aviation Industry, among which the average value of two years before and after was used because the data for 2008 is missing. The data on civil aviation passenger traffic and cargo traffic volume were from the China Statistical Yearbook.

The prediction model in this paper is based on the GM(1,N) model, the construction of which focuses on the structure of matrix B. If the value difference between rows or columns of matrix B is too large or too close, the condition number of the matrix will be too large, resulting in the phenomenon of matrix drift. The order of magnitude and dimension of the original data selected influencing factors are quite different, so it is necessary to preprocess the original data, such as initializing. After initializing the original data, the predicted value is obtained, and then the predicted value is initialized by inverse transformation so as to restore the order of magnitude and dimension of the predicted sequence. Through the above-screened influencing factors, the Pearson correlation coefficient is used to analyze the correlation between each impact factor and civil aviation carbon emissions, and the results are shown in Figure 2.

As shown in Figure 2, the diagonal plot shows the distribution of China’s civil aviation carbon emissions and the six influencing factors themselves; the upper triangle shows the corresponding Pearson correlation coefficients, and the lower triangle shows the corresponding scatter plots. The Pearson correlation coefficient is −0.31583, while the other five factors are strongly correlated with China’s civil aviation carbon emissions, with Pearson correlation coefficients greater than 0.9, and these five factors are also strongly correlated with each other. Therefore, in summary, this paper selects five influencing factors of carbon emissions in China’s civil aviation, which are civil aviation passenger traffic, civil aviation cargo volume, total civil aviation turnover, civil aviation fuel consumption, and civil aviation industry-wide operating income.

4.3. Empirical Results

According to the five influencing factors screened out above, the background values and order optimization of different models were shown in

Table 2. Background values and order of the different models.

	GM(1,5)	OGM(1,5)	OBGM(1,5)	FOBGM(1,5)
Background value	0.5	0.5	1	1
Order	1	1	1	1.0907

, then GM(1,5), OGM(1,5), OBGM(1,5), and FOBGM(1,5) were used to forecast the carbon emissions of civil aviation in China, and the prediction accuracy of the four models was analyzed. In this paper, we collected relevant data from 2003 to 2017 and used the data from 2003 to 2012 as the original data to forecast China’s civil aviation carbon emissions from 2013 to 2017 and then compared it with the actual data to obtain the prediction accuracy of the model. The prediction results are shown in

Table 3. Prediction results of carbon emissions.

Year	Civil Aviation Carbon Emissions Original Data	GM Model Background Value/Order		OGM(1,5)		OBGM(1,5)		FOBGM(1,5)
		Background Value	Cumulative Order	Background Value	Cumulative Order	Background Value	Cumulative Order	Background Value	Cumulative Order
		0.5	1	0.5	1	1	1	1	1.0907
2013	6302.916254	6303.041342		6302.988503		6302.987549		6302.916802
2014	6990.086895	6990.332886		6990.260412		6990.259593		6990.250567
2015	7901.107325	7901.670548		7901.5331		7901.531429		7901.510396
2016	8879.055003	8879.86047		8879.67912		8879.676821		8879.645614
2017	10,009.66542	10,010.64701		10,010.4437		10,010.44095		10,010.40286

.

Based on the prediction results, the absolute and relative errors of each model for each year were obtained, as shown in Figure 3 and Figure 4.

It can be found in Figure 3 and Figure 4 that the prediction accuracy of the model gradually improved during the improvement of the model, and the FOBGM(1,5) model had the highest prediction accuracy among the four models, and the prediction error of the FOBGM(1,5) model gradually decreased by 99.56%, 33.46%, 28.43%, 26.67%, and 24.87%, relative to the GM(1,5) model from 2013 to 2017. In order to further simplify the prediction errors of each model and more directly reflect the prediction accuracy of the model, this paper compares the prediction performance of each model through MAPE, MAE, RMSE, and MAE.

As can be seen from Figure 5, through multi-stage optimization and improvement of the model, the prediction accuracy of the model is gradually improved. The prediction error of the OGM(1,5) model was reduced by 24.40% compared with the GM(1,5) model. The prediction error of the OBGM(1,5) model was reduced by 0.43% compared with OGM(1,5) model. The prediction error of the FOBGM(1,5) model was reduced by 9.48% compared with that of the OBGM(1,5) model, which proves the effectiveness of the model improvement, and the FOBGM(1,5) model had the smallest prediction error and the highest prediction accuracy among the four models.

5. Conclusions

In this paper, based on the GM(1,N) model, first, the mechanism defects, parameter defects, and structural defects of the GM(1,N) model were compensated by adding linear correction term and grey action, and the OGM(1,N) model is established. Then, the background value coefficients of the OGM(1,N) model were optimized using the PSO algorithm, and the OBGM(1,N) model was established. Then, by introducing the fractional order idea, this paper uses the PSO algorithm to optimize the cumulative order of the OBGM(1,N) model and extend the order of the OBGM(1,N) model from an integer field to a real number field to establish FOBGM(1,N) model. Five influencing factors were determined to predict China’s civil aviation carbon emissions using civil aviation passenger traffic, civil aviation cargo volume, total civil aviation turnover, civil aviation fuel consumption, civil aviation industry-wide operating income, and civil aviation transportation intensity. Based on the carbon emission prediction data of civil aviation transportation in China, the improvement effect of each model was empirically studied. It can be seen from the prediction results that the prediction error of the model decreases gradually, and the prediction accuracy shows an increasing trend after multi-level and multi-angle improvement and optimization. Among them, the MAPE of the OGM(1,5) model decreased by 24.40% compared with the GM(1,5) model, the MAPE of the OBGM(1,5) model decreased by 24.72% compared with the GM(1,5) model, the MAPE of the FOBGM(1,5) model decreased by 31.86% compared with the GM(1,5) model. It reflects the effectiveness of model improvement and proves the practicability of the FOBGM(1,5) model. From the perspective of algorithm improvement, the idea of improving the GM(1,N) model based on fractional summation and background value optimization has certain practical significance and application prospects and can be applied to the optimization research of other prediction models of grey system theory. Although the paper has carried out multi-angle and multi-level optimization of the model, such as background value, fractional order, linear increase, etc., the model still has room for improvement, such as from the perspectives of initial value improvement, residual correction, and metabolism, as a subject of further research.