Considering the complexity pattern of the gold price, this paper adopts the decomposition-reconstruction-forecast-mergence scheme to perform the international gold price forecast. The original gold price data are decomposed into 12 intrinsic mode functions and a residual by the complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) method, and then the 13 sequences are reconstructed into a high-frequency subsequence (IMFH), a low-frequency subsequence (IMFL), and the residual (Res). According to the different characteristics of the subsequences, the IMFL and Res are forecasted by the support vector regression (SVR) model. Besides, in order to further improve the prediction accuracy of IMFH, we have developed a novel hybrid method based on the support vector regression (SVR) model and the grey wolf optimizer (GWO) algorithm with SVR for predicting the IMFH of gold prices, i.e., the CEEMDAN-GWO-SVR model. This hybrid model combines the methodology of complex systems with machine learning techniques, making it more appropriate for analyzing relationships such as high-frequency dependences and solving complex nonlinear problems. Finally, the final result is obtained by combining the forecasting results of the three subsequences. The empirical results show that the proposed model demonstrates the highest prediction ability among all of the investigated models in a comparison of prediction errors with other individual models.

1. Introduction

As a precious metal, gold also acts as a reserve. The central banks of many countries hold an amount of gold, for the purpose of mitigating potential risks, resisting inflation, and increasing economic growth. Gold plays an important role in reserve assets. Especially in the last few years, the international situation has undergone complex and deep changes due to the continuous escalation of the China-US trade war. In this context, it is of great significance to mitigate international financial risks and maintain economic stability. Central banks’ demand for gold has been increasing. In 2021, central banks added a total of 463 tons of gold, 82% more than in 2020, bringing global central bank gold reserves to the highest level in nearly 30 years.

The trend of the gold price is different from other minerals [1], which is impacted by various factors, including demand and supply, inflation, and political issues, among others. In the financial crisis of 2008, most commodities and equities dropped about 40%, while the gold price remained at a high point and increased about 6% [2]. With the continuous adjustment of international patterns in recent years, gold is facing a more complex market environment [3], and its price has been subject to severe fluctuations. Figure 1 shows the movement of the gold price for 21 years from 2000 to 2021. Obviously, the price of gold fluctuates up and down, showing a high degree of nonlinearity and uncertainty. Therefore, an accurate gold price forecast is generally considered to be challenging [4, 5].

The main purpose of this research is to contribute to the accurate prediction of gold price and movement. To this end, we adopt the decomposition-reconstruction-forecast-mergence scheme to perform the international gold price forecast, and propose the hybrid CEEMDAD-GWO-SVR model. Our experiments demonstrate that this scheme and the hybrid model can significantly improve the forecasting performance. The rest of this paper is organized as follows: in Section 2, a series of recent works on gold price forecasting are presented, in Section 3, the CEEMDAN method, the GWO algorithm, and the SVR model are recalled for completeness, and the main steps of the hybrid CEEMDAD-GWO-SVR model are introduced, in Section 4, the experimental results are presented and forecasting accuracies are compared across different models, and in Section 5,the conclusions are discussed.

2. Literature Review

2.1. Forecasting Models

Research on gold price forecasting was performed for decades, and numerous methods have been proposed. The commonly used methods mainly include multiple regression, time-series analysis, neural network methods, and hybridization or combination models.

The multiple regression model makes prediction of gold prices with multiple influencing factors. The time-series methods forecast the future gold price based on the past time series, including ARIMA (auto-regressive integrated moving average), the ARCH (auto-regressive Conditional Heteroscedasticity) class models, and the GARCH (Generalized Auto-Regressive Conditional Heteroscedasticity) class models among others. For example, Wang and Xia [6] analyzed the marginal and joint impact of economic variables on the gold price, and provided a short-term forecast of the gold price by multiple regression, achieving a satisfactory result. Yang [7] used the ARIMA model, and the results showed that the model can reflect the true gold price trend to a certain extent. However, the above two methods have certain limitations: the explanatory variables of the multiple regression are difficult to choose in empirical research, the time-series methods have too many parameters to be estimated, and the optimal lag order is hard to choose. In practice, these two methods may fail to characterize the nonlinear and complex behavior of the gold price time series and decrease the forecasting accuracy.

The neural network methods are well suited to address nonlinear and nonstationary problem, and have been widely applied in time-series forecast. Parisi et al. [8] analyzed recursive and rolling neural network models to forecast one-step-ahead sign variations in gold price, and they found that the rolling ward networks perform well in forecasting gold price sign variation. Zhang and Ci [9] proposed a deep belief network (DBN) model for gold price forecast, and the empirical results showed that the proposed DBN model outperforms the BP (back propagation) neural network model, genetic algorithm optimization for BP model (GA-BP), and the ARIMA model. However, the neural network models have transparency problem and it is difficult to understand the economic implications and impact mechanisms [10].

In recent years, research on gold price forecast has focused on two types of hybrid models and combined models, such as (1) the hybridization or combination of different models and (2) the hybridization or combination of models with metaheuristic algorithms or advanced intelligent approaches. Since not every model is suitable for hybridization or combination with others, the second type of hybridization has more applications than the first one. Most importantly, the hybridization of machine learning approaches achieved many successes and provided more useful findings about the gold price forecast. Weng et al. [11] proposed a novel online extreme learning machine (OS-ELM) model based on genetic algorithm (GA) and regularization to forecast gold price. The proposed hybrid model solves the singularity problem of OS-ELM model and outperforms the ARIMA, BP, ELM, and OS-ELM models. Alameer et al. [12] presented a hybrid model for accurately forecasting long-term monthly gold price fluctuations. The model uses whale optimization algorithm (WOA) as a trainer to learn the multilayer perceptron neural network (NN), and their empirical results indicate the superiority of hybrid WOA-NN model over the classic NN, particle swarm optimization for NN (PSO-NN), genetic algorithm for NN (GA-NN), and ARIMA models.

Owing to its outstanding forecasting performance and interpretability, the support vector regression (SVR) model, especially as part of a hybrid model, has been extensively used in landslide displacement prediction [1315], daily traffic peak flow management [16, 17], financial investment [18, 19], air quality prediction [20, 21], and gold price forecast [22, 23] as well. Compared with other machine learning methodologies, SVR shows rapid convergence [24] and efficiently solves multidimensional function estimation problems [25], but the selection of the SVR penalty parameter and the kernel function parameter has an important influence on its own performance [26]. Especially, SVR would have a large forecasting error in dealing with nonlinear high-frequency time series. Results in [26] showed that the hybrid models of SVR and metaheuristic algorithms can significantly improve the performance of the single SVR model.

Metaheuristic algorithms are widely used by researchers for their remarkable ability to solve optimization problems. Since the commencement of the classical metaheuristics, namely, genetic [27], simulated annealing [28], tabu search [29], and particle swarm optimization [30] in the early 1970s to late 1990s, a new generation of metaheuristics has grown exponentially. In the past three years, nearly 20 metaheuristics has been proposed, including Monarch butterfly optimization [31], butterfly optimization algorithm [32], shuffled shepherd optimization algorithm [33], and water strider algorithm [34]. Currently, there are more than 200 metaheuristic techniques in the literature and several surveys [3537] provided comprehensive overviews on them. Grey Wolf Optimizer (GWO) [38] is one of metaheuristic algorithms, which was proposed by Mirjalili in 2014. It is a swarm intelligent algorithm that mimics the leadership hierarchy and hunting mechanism of grey wolves in nature. It has the advantages of fast convergence speed and strong search ability and can simultaneously optimize the parameters and weights of the involved kernels [39]. Recent studies [38, 40, 41] have demonstrated that GWO provides competitive results compared to some well-known metaheuristics. Therefore, this paper hybridizes SVR with GWO algorithm (i.e., PSO-SVR) to forecast the nonlinear high-frequency part of the gold price time series and the subsequent experimental results verify the satisfactory learning precision and generalization ability with PSO-SVR.

2.2. Data Decomposition

Studies [4245] have shown that data decomposition before forecasting can effectively reduce the nonstationary characteristics of the time series, and it can improve the forecasting performance significantly. Current decomposition algorithms include empirical mode decomposition (EMD) [42], ensemble empirical mode decomposition (EEMD) [43], complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) [44], and variational mode decomposition (VMD) [45]. EMD is an adaptive decomposition method for nonlinear and nonstationary processes. This method can decompose a complicated dataset into a finite and often small number of intrinsic mode functions (IMFs). However, EMD may encounter the mode-mixing problem caused by intermittent signals. EEMD effectively resolves this problem through an ensemble of the signal plus the white Gaussian noise. However, the added white Gaussian noise also creates some of the following new problems: the reconstructed signal includes residual noise and different realizations of signal plus noise may produce different number of IMFs [44]. To overcome these problems, CEEMDAN was proposed. It performs the EEMD over white Gaussian noise added at each stage of the decomposition and makes the decomposition complete. VMD is an entirely nonrecursive variational mode decomposition model, and the experimental results in the study [45] showed that VMD is more robust to sampling and noise than EMD and EEMD. However, the influence parameters of VMD need to be set in advance [4648]. In summary, CEEMDAN not only provides a better spectral separation of the IMFs than EMD and EEMD but also overcomes the parameter selection limitation of VMD. In addition, CEEMDAN reduces the computational cost due to its fewer sifting iterations. Considering the importance of data decomposition and the ease of implementation of the CEEMDAN method, this paper uses CEEMDAN to preprocess the data before forecasting.

2.3. Main Contributions

In this study, we first adopt CEEMDAN to conduct the data preprocessing process and then use GWO-SVR model to forecast the price of gold. The proposed model, named CEEMDAN-GWO-SVR, shows better forecasting ability and performance than the state-of-the-art benchmark models in the numerical experiments. All in all, a list of abbreviations and variables has been added in Table 1 in the appendix. The principal contributions of this paper are as follows:(1)Although numerous methods have been proposed to forecast the gold price, especially the hybrid or combined models were performed for decades, there are few studies based on “decomposing first and then merging.” This paper adopts a decomposition-reconstruction-forecast-mergence scheme to forecast the international gold price, using the CEEMDAN method to decompose the original data into several IMFs and the residual and then using T-test to reconstruct them into three subsequences, namely, high-frequency subsequence (IMFH), low-frequency subsequence (IMFL), and the residual (Res), which are, respectively, forecasted by SVR model and finally the three forecasting results are combined to the final result.(2)Considering the different characteristics of the three subsequences and the large error of SVR when dealing with nonlinear high-frequency time series, different SVR models are used here. Specifically, for IMFL and Res, this paper uses SVR with a linear kernel function to forecast the price of gold, and for IMFH, the GWO algorithm and the SVR model with the Gaussian kernel function are hybridized to generate the GWO-SVR, which overcomes the problem of subjective selection of parameters to improve the accuracy of forecast.(3)A comparison demonstrates that the CEEMDAN method can effectively capture effective information from the original gold price data, and the proposed CEEMDAN-GWO-SVR model provides significant improvements over the benchmark models.

3. Materials and Methods

3.1. Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN)

Based on EMD and EEMD, the CEEMDAN method adds white Gaussian noise at each stage of the decomposition and then computes the unique residue to obtain each IMF of the original dataset. The original signal can be exactly reconstructed by summing the IMFs owing to the complete decomposition of this method. The detailed algorithm is as follows:Step 1. We add white Gaussian noise () to the original dataset and then generate the new datasetwhere is the number of the trials and is the signal-to-noise ratio.Step 2. Each () is decomposed by EMD by getting their first modes and computing the first IMF of Step 3. We calculate the first residual, , asStep 4. We define the operator , which produces the -th mode by EMD and then decomposes () until we get the second IMF,Step 5. We calculate the -th () residue and , in equation (5), where is the total number of IMFStep 6. We decompose () by EMD and compute the -th IMF asStep 7. We repeated Step 4 to 6 until the residue has less than two extrema and cannot be decomposed any longer, so the final residue can be expressed aswhere the original dataset is decomposed

3.2. Grey Wolf Optimizer (GWO)

GWO is metaheuristic algorithm, which mimicks the leadership hierarchy and hunting mechanism of grey wolves in nature. Four types of grey wolves, including , , , and , are employed for simulating the social dominant hierarchy. wolves are the leaders and the decision makers in the pack. wolves are the second group in the hierarchy that help wolves in decision-making and establish a discipline for the pack. wolves obey and wolves, but dominate wolves. They play the roles of scouts, sentinels, elders, hunters, and caretakers. wolves, the lowest in hierarchy, always have to obey all the other types of wolves. The GWO algorithm mathematically models this social hierarchy and hunting process to collect randomly preferred solutions. Specifically, , , , and are considered as the fittest, the second preferred, the third preferred, and the rest solutions, respectively. The GWO algorithm contains the following three steps: tracking, encircling, and attacking the prey.

In the encircling, grey wolves encircle prey during the hunt, and the mathematical formulation is given bywhere is the number of the current iteration, and are the positions of a grey wolf and the prey, respectively, and are coefficient vectors, components of are linearly decreased from 2 to 0, and and are random vectors in [0, 1].

In hunting, it is supposed that the , , and have better knowledge about the potential location of prey. Therefore, the current first three best solutions are saved, and the other search agents (including the ) update their positions according to the position of the best search agent. The formulas are given by equations (10)–(13):

The grey wolf agents diverge in search of the prey and converge for attacking the prey, which can be mathematical modelled by the value of . The grey wolf agents tend to diverge from the prey when and converge towards the prey when . Finally, the GWO algorithm is terminated by minimizing of the fitness function, which would be discussed in Section 3.4.

3.3. Support Vector Regression (SVR)

SVR is a supervised machine learning algorithm. Based on structural risk minimization, SVR can solve nonlinear problems with high dimensionality and small samples. For a set of training data with N observations , where are input vectors, are output vectors, and the SVR function can be described as the following linear regression function:where is the feature mapping function and the parameters and are the weight and bias, respectively. Based on the structural risk minimization, and can be calculated as follows:where is a cost function, and are slack variables, and is an insensitive loss function. By introducing the Lagrange multipliers and , the SVR function is expressed aswhere is the kernel function, which is the inner product of the feature mapping function.

SVR would have a large forecasting error when dealing with nonlinear high-frequency time series. Therefore, unlike the linear kernel function used in forecasting the low-frequency time series, the Gaussian kernel, , is used to forecast the high-frequency time series, owing to its advanced mapping capability and stable forecasting performance. Meanwhile, the selection of the SVR penalty parameter and the kernel function parameter has an important influence on its forecasting accuracy [26]. Thus, this paper tunes these SVR parameters using the GWO algorithm.

3.4. The Hybrid CEEMDAD-GWO-SVR Model

The details of the hybrid CEEMDAD-GWO-SVR model are as follows, and the flowchart is presented in Figure 2. It should be noted that the “decomposition-reconstruction-forecast-mergence” scheme is also shown in Figure 2.Step 1. We conduct CEEMDAN to decompose the original gold price data into several relatively simple IMFs and the residual by using equations (1)–(8). This step corresponds to the “decomposition” of the scheme.Step 2. We reconstruct the IMFs and the residual into three subsequences, namely, IMFH, IMFL, and Res. Specifically, since the high-frequency IMFs are affected by some short-term random factors, they can be assumed to fluctuate around the zero mean. Therefore, each IMF is reconstructed into high-frequency subsequence (IMFH) and low-frequency subsequence (IMFL) by T-test with zero mean. The residual is classified as the third subsequence, i.e., Res. This step corresponds to the “reconstruction” of the scheme.Step 3. Forecasting IMFH by GWO-SVR and forecasting IMFL and Res by SVR. The optimal GWO-SVR model is determined by minimizing the fitness function. In this paper, the mean square error, given by equation (17), is used to calculate the forecasting accuracy. This step is decisive for the prediction performance of the proposed model, corresponding to the “forecast” of the scheme.Step 4. This is the last part of the scheme. By summing up the forecasting results of these three subsequences, the final result is obtained.

4. Empirical Study

This study selects the daily closing prices of international gold from January 3, 2000 to October 21, 2022, getting from the World Gold Council website (https://www.gold.org), with a total of 5,950 data points. Figure 3 displays the fluctuation diagram of the gold price. We take the first 5,655 data points (from January 3, 2000 to September 3, 2021) as the training set, and the last 295 data points (from September 5, 2021 to October 21, 2022) as the test set. Table 2 presents details of each set.

4.1. Data Decomposition and Reconstruction
4.1.1. Data Decomposition

In the CEEMDAN, the parameters of ensemble size and the signal-to-noise ratio are set as 500 and 0.2, respectively, and the maximal number of iterations is 2,000, as proposed in [44]. Decomposition results is shown in Figure 4. The gold price is decomposed into 13 sequences with 12 IMFs and one residual. The frequency of sequences is ranked from high to low. The high-frequency IMFs fluctuate relatively randomly but basically around zero mean, which can be regarded as the interference of short-term random factors (e.g., psychological factors, speculative factors, and short-term policies). The low-frequency IMFs have roughly similar trend to that of the original gold price data, reflecting the influence of important factors such as market policies and macroeconomics. The residual shows the overall upward trend, determining the long-term trend of gold price.

Meanwhile, Pearson correlation coefficient, Kendall rank correlation coefficient, and the variance contribution rate (VCR) are selected to measure the relationship between each sequence and the original gold price data, as depicted in Table 3. From the two correlation coefficients, the high-frequency IMFs have a low association with the original gold price data, while the low-frequency IMFs and the residual have a high correlation with the original. The similar pattern is observed for VCR and the residual covers 66.6% of all variance, indicating that it has the greatest influence on the original data and reflects the internal trajectory of the gold price data. To sum up, the trend of gold price is mainly described by the residual. The low-frequency IMFs have a great influence on gold price in a certain period, while the high-frequency IMFs only have a small influence in a short period.

4.1.2. Data Reconstruction

The high-frequency IMFs have a small influence on gold price in a short period, so they can be combined into one sequence to improve the efficiency of modeling, but the problem is how to determine which sequences are high-frequency IMFs. Since the high-frequency IMFs are affected by some short-term random factors, they can be assumed to fluctuate around zero mean. Therefore, IMFs can be reconstructed by T test with zero mean for each IMF. The results are shown in Table 3. The values of IMF1∼IMF9 are greater than 0.01.

Since the nine IMFs fail to reject the null hypothesis for t test, so they are all divided into the high-frequency subsequence (IMFH) and . The remaining sequences, IMF10, IMF11, and IMF12, reject the null hypothesis for t test at 1% significance level, so we define the low-frequency subsequence (IMFL) as . The residual is classified as the third subsequence of the gold price data, i.e., Res. Table 4 presents the reconstructed results.

Figure 5 shows the reconstructed subsequence trend, and it can be found that: (1) IMFH fluctuates relatively randomly but basically around zero mean, which is mainly interfered by some short-term random factors; (2) IMFL has roughly similar trend to that of the original gold price, which is mainly affected by market policies and macroeconomics; and (3) Res shows the overall upward trend, conducting the trend of gold price.

4.2. Forecasting Results and Discussion
4.2.1. Evaluation Measures

The forecasting accuracies of the proposed model and the benchmark models are evaluated with the following three measures: the mean absolute percentage error (MAPE), the root mean square error (RMSE), and the directional accuracy (DA). Their formulas are given by:where is the number of forecasting points, is the actual gold price at point , and is the gold price forecast at point . Obviously, the smaller the MAPE and RMSE values, the higher the forecasting accuracy of the model, and the larger the DA values, the higher the direction accuracy of the model.

4.2.2. In-Sample Forecast (January 1, 2016–September 3, 2021)

The original gold price data is decomposed and reconstructed into three subsequences: IMFH, IMFL, and Res, which can be expressed aswhere is the actual gold price. The three subsequences are forecasted separately, and then combined to obtain the final forecast of the gold price. Considering the low forecasting accuracy of SVR when dealing with nonlinear high-frequency time series, IMFH is forecasted by GWO-SVR to improve the parameter selection, and IMFL and Res forecasts are conducted by SVR.

In GWO algorithm, the optimization range is set at (0.001 and 200), and the population size and the maximum number of iterations are set to 20 and 50, respectively. For IMFH, the optimal penalty parameter and the kernel function parameter of SVR are −132.7510 and −1.8545, selected by GWO. All calculations are done on the Ryzen R7-4800u CPU, with 16.0 GB DDR4 3200Mhz RAM, and windows 11 system.

Figure 6 shows the forecasting results of the three subsequences and the original gold price data. It can be observed that the blue lines (forecasting results) are very close to the red lines (the reconstructed subsequences), which indicate that the proposed CEEMDAN-GWO-SVR model achieves a relatively good forecasting result. Table 5 presents the measures of the forecasting accuracy. IMFL and Res have relatively small errors and high directional accuracy, achieving good forecasting performance. IMFH has lower accuracy than the other two, and it also reduces the forecasting accuracy of the original gold price data. This is because IMFH is affected by random factors, which increases the difficulty of forecasting.

Table 6 presents the accuracy measures of the forecasting results of the proposed CEEMDAN-GWO-SVR model, SVR-based models, and non-SVR-based models. Among them, the SVR-based models include CEEMDAN-SVR, GWO-SVR, and SVR, and the non-SVR-based models include ARIMA, GARCH, ANN, CART, and BPNN (back propagation neural network). Obviously, the proposed CEEMDAN-GWO-SVR model is superior to the alternatives. For the SVR-based models, the CEEMDAN-SVR, CEEMDAN-PSO-SVR, and CEEMDAN-GA-SVR models can certainly improve the performance of SVR. While the hybrid CEEMDAN-GWO-SVR model achieves the best forecasting result, which shows that the proposed model can critically capture effective information from the original data, and it also demonstrates the superiority of the decomposition-reconstruction-forecast-mergence scheme. Figure 7 shows the convergence plots of GWO, which reflected that the objective can be minimized by the algorithm before the maximum number of iterations.

4.2.3. Out-of-Sample Forecast (September 6, 2021–October 21, 2022)

The above results prove that the CEEMDAN-GWO-SVR model outperforms the benchmark models for in-sample forecast. In order to further test the proposed model for out-of-sample forecasts, this paper estimates the gold price for the next 295 consecutive days (September 6, 2021–October 2021). Parameters of the CEEMDAN-GWO-SVR model here are set in the same way as above.

Out-of-sample forecasting accuracy indexes comparison are respectively shown in Table 7. The superiority of the proposed CEEMDAN-GWO-SVR model can also be demonstrated. Among the three SVR-based models, the CEEMDAN-SVR model obtains the best performance, and at the same time, it also has superior forecast performance than the non-SVR-based models. This shows that the CEEMDAN method could improve the forecasting performance. Moreover, when CEEMDAN-SVR and GWO are hybridized, the forecasting accuracy of gold price is improved further, which shows the superiority of the CEEMDAN-GWO-SVR model proposed in this paper. It should be noted that in short-term forecasts, although the CEEMDAN-GWO-SVR model has relatively lower directional accuracy, it still has smaller forecasting errors (MAPE and RMSE) than alternatives. Therefore, the proposed CEEMDAN-GWO-SVR model can be considered as the best model.

5. Conclusions

This paper proposes a decomposition-reconstruction-forecast-mergence scheme to forecast the international gold price. In this scheme, a useful gold price forecasting model is proposed by combining the CEEMDAN method with GWO-SVR. The experimental results verify the validity of the CEEMDAN-GWO-SVR model for gold price forecast. The following conclusions are drawn:(1)Due to the complex pattern of the gold price, data preprocessing in the early stage of forecasting is necessary. The CEEMDAN method adds white Gaussian noise at each stage of the decomposition and sufficiently reduces the noise of the original data. Also, experimental results show that the CEEMDAN method can improves the forecasting performance of SVR.(2)Although the CEEMDAN-SVR model can certainly improve the performance of SVR, the forecasting accuracy can be further improved by hybridizing CEEMDAN-SVR and GWO. This is because GWO optimizes the selection of SVR parameters for high-frequency time series decomposed by CEEMDAN, thus ultimately improving the forecasting performance.(3)The CEEMDAN-GWO-SVR model significantly outperforms the other six benchmark models, which also demonstrates the superiority of the decomposition-reconstruction-forecast-mergence scheme.

However, this paper still has certain limitations. This paper only uses time series to forecast gold prices. Although the proposed model has achieved good forecasting results, the research on the main influencing factors of gold price is still an important research issue.

Further research can be carried out from the following aspects in the future. The GWO algorithm has certain problems in population diversity and the balance between the exploitation and exploration [49]. Therefore, the algorithm could be improved to further improve the prediction performance of the SVR model. Moreover, except GWO, finding a suitable approach to hybridize more recent metaheuristic algorithms with SVR is also an important research issue.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


This work was supported by the National Science Foundation of China under Grants 71771187 and 72163029, and the Fundamental Research Funds for the Central Universities in China under Grant JBK190602.