1 Introduction

Remarkable progress has been made in tropical cyclone (TC) track forecasting in recent years, owing to advanced assimilation methods with multi-source observations, improved understanding of atmospheric physical processes and the air–sea interaction, and the application of ensemble forecasting using multiple models. How-ever, due to limited knowledge on the cloud physics associated with TCs and the various uncertainties in the air–sea interaction process (Emanuel, 2000), there have been fewer improvements in TC intensity forecasts than in their associated track forecasts (Qian et al., 2012). Therefore, TC intensity forecasting remains a challenge for the international research and operational forecasting community. After several strong hurricanes (e.g., Katrina, Rita, and Wilma) impacted the United States and caused serious economic losses and numerous casualties in 2005, the Hurricane Forecast Improvement Project (HFIP) was initiated by the National Oceanic and Atmospheric Administration (NOAA) in 2008. HFIP aimed to reduce the average error in TC intensity forecasts by 20% in its first 5 yr and by 50% by the end of its 10-yr duration; it also was designed to improve the efficiency and accuracy of probabilistic forecasts of rapidly intensifying hurricanes (Gall et al., 2013).

Currently, statistics, statistical dynamics, and ensemble techniques are widely applied to TC intensity forecasting by operational forecasting centers worldwide (Chen et al., 2012; Yu et al., 2012). In 2004, based on the Statistical Hurricane Intensity Prediction System (SHIPS), the National Hurricane Center (NHC) implemented ensemble forecasting of hurricane intensity in the operational system for eastern Pacific and Atlantic (DeMaria, 1996; DeMaria and Kaplan, 1999). High skill scores were achieved for long lead-time forecasts by employing the Logistic Growth Equation Model (LGEM) and Decay-SHIPs in ensemble forecasts from the models such as the Hurricane Weather Research and Forecasting (HWRF), Geophysical Fluid Dynamics Laboratory (GFDL), Glo-bal Forecast System (GFS) of NCEP, and ECMWF (added in 2015) models (DeMaria, 1996; DeMaria and Kaplan, 1999; DeMaria et al., 2005; Wang et al., 2015). Yu et al. (2015) also demonstrated encouraging operational applications of ensemble forecasts of TC intensity. They calibrated the multiple model results with a consensus method by selecting major environmental variables with a step-wise regression.

Introducing probabilistic forecasts from traditional deterministic forecasts is a current trend in weather forecasting (Zhu, 2010; Slingo and Palmer, 2011). Simple deterministic forecasts from numerical models do not consider the objective forecast uncertainty that limits the transfer of complete information and proper guidance to forecasters and users (Krzysztofowicz, 1985). Probabilistic forecasts, in contrast, provide quantitative expressions of forecast uncertainty (Kelly and Krzysztofowicz, 1997) and more complete and comprehensive information about future weather conditions. In recent years, the ensemble-based Bayesian method has widely been used to forecast surface temperature, surface pressure, and precipitation (Herr and Krzysztofowicz, 2005; Raftery et al., 2005; Sloughter et al., 2007; Zhi et al., 2014, 2015), among others. However, no advanced studies exist that forecast TC intensity using an ensemble-based Bayesian method. In this study, we examine whether an ensemble-based Bayesian method could improve TC intensity forecasts through multi-model application. These predictions from the Bayesian model averaging (BMA) method could provide probabilistic (and/or uncertainty) forecasts rather than deterministic forecasts. Notably, the forecast uncertainty is quantitatively expressed, which allows for advanced scientific interpretation of the forecasts (Zhu, 2005).

The TC intensity forecasts and observational datasets are described in Section 2. The error analysis and bias correction are discussed in Section 3. The BMA method and its application are summarized in Section 4 and 5, respectively. Conclusions are then made in Section 6.

2 TC intensity forecasts and observational datasets

The locations and intensities of an observed TC were determined from the tropical cyclone best-track dataset for the western North Pacific, available from the Joint Typhoon Warning Center (JTWC, http://www.usno.navy.mil/NOOC/nmfc-ph/RSS/jtwc/best_tracks/wpindex.php).

Forecast data from three models, including two global models and a regional model (Table 1), were obtained from the database created by the NCEP’s Environmental Modeling Center (NCEP/EMC). The global ECMWF Integrated Forecasting System (ECMWF-IFS) (Morcrette et al., 2009) utilizes four-dimensional variational data assimilation and had a resolution of 9 km in 2016. It is truncated at wavenumber 1279 and divided into 91 hybrid layers perpendicular to the vertical direction, with a model top pressure of 0.01 hPa. The model is initialized twice per day at 0000 UTC and 1200 UTC. NCEP-GFS is also a global model that was developed by NCEP, and it uses the Global Data Assimilation System. The model was truncated at wavenumber 574 in 2014 and at 1534 in 2015 and 2016. The GFS consists of 64 hybrid layers perpendicular to the vertical direction, with a model top pressure of 0.27 hPa. Its resolution was 22 km in 2014 and 13 km in 2015 and 2016. It is initialized at 0000, 0600, 1200, and 1800 UTC each day. The NCEP-HWRF model adopts the Gridpoint Statistical Interpolation assimilation system. It has 61 layers, a top pressure of 0.27 hPa, and an inner domain resolution of 2 km. Its forecast is initialized four times per day.

TC intensity is normally represented by the maximum wind speed near the cyclone center (V_max) or the minimum sea-level pressure at the cyclone center (P_min). Research has shown that P_min has more significant linear relationships with initial error and forecast error than does V_max (Knaff and Zehr, 2007; Yu et al., 2013). Therefore, this study used P_min as its TC intensity indicator. This study focused on 24- to 120-h forecasts for all TCs that occurred in the western North Pacific basin during 2014–16, as produced by the three aforementioned models at 6-h intervals. TC intensity forecast error refers to the difference between the analyzed value from the JTWC best-track data and the P_min forecast by the models.

3 Error analysis and bias correction

Error analysis of the forecasts from the aforementioned models for 2014–16 was performed in an attempt to capture the forecast error characteristics. Prior to the BMA adjustment, the systematic errors for the indivi-dual forecast models were corrected to improve forecast reliability.

Forecast error is calculated as the difference between the annual average of the observed and forecast P_min value at each lead time. The results are presented in Fig. 1. Both the GFS and HWRF include the vortex relocation technique, while the HWRF also includes a bogus TC technique. TheP_min values predicted by the HWRF mo-del for 2014 and 2015 are about 10 hPa smaller than the observed values. In 2016, the model’s forecast bias was –2.4 hPa, which was a significant decline that can likely be attributed to major model changes. In particular, the error in the 120-h forecast was weak and changed signs (Fig. 1, third bar). This suggests that the systematic forecast errors in 2014 and 2015 are different from those in 2016. The GFS forecasts for P_min in 2014 were 10–20 hPa higher than the corresponding observed values, while the GFS forecasts for 2015 were 12–20 hPa lower than the observations. Errors in the GFS forecasts for 2016 had absolute values no greater than 3 hPa. There is no consistent sign of model error for these three years. Therefore, no correlation was found between this model’s forecast errors throughout the three years. The European Centre (EC) model was the only one that de-monstrated consistently negative forecast error (or bias) throughout the three years. As the lead time varied, the forecast errors in 2014 and 2015 showed significant correlations with the forecast error in 2016. Therefore, bias correction for the EC’s P_minforecasts in 2014–15 may substantially improve its forecasts for 2016, while the GFS and HWRF are unsuitable for bias correction, due to lack of a systematic pattern of error variation during the three years.

The analysis above reveals that only the EC model’s forecast errors followed a regular pattern of distribution across the years, and its P_minforecasts were relatively high overall (indicating a weak TC). Bias correction was applied to the EC forecasts for 2016 based on the pattern of error distribution across 2014–15 from the following formula: P_min(bc; t) = P_min(t) – Bias(t), where t represents forecast lead time. This was achieved by subtracting 17, 15, 16, 13, and 11 hPa, respectively, from every 24- to 120-h forecast. Figure 2 compares the root-mean-square error (RMSE) of the forecasts from the three models and the bias-corrected EC, denoted as EC-BC. The P_min forecasts from the regional model HWRF were much more accurate than those from the global models EC and GFS, as indicated by the comparatively lower errors at all lead times (Fig. 2). This is possibly due to its higher resolution and the use of a bogus TC. Compared to the RMSE of the EC forecasts, the RMSE of the EC-BC forecasts were significantly reduced: by 4.6 hPa at a lead time of 120 h and by 7–8 hPa at other lead times.

We also plotted the error distribution map for the three original models (not presented in this paper) to identify any regularity in the error distribution across latitudes and regions for further removal of systematic errors. Unfortunately, those error distributions were more random and provided little useful information.

4 Bayesian model averaging 4.1 The theoretical basis for BMA

Bayesian model averaging (BMA) is a method of data post-processing based on Bayesian theory. It has multiple advantages over other statistical prediction models. For example, it considers subjective prior information and unavoidable unpredictability in models, allowing for the effective and accurate extraction of information from model forecasts. It has been successfully applied to temperature, pressure, wind speed, and precipitation forecasts.

There are several key steps in BMA. The first step is to calculate posterior probability. The second step is to derive the posterior probability density function (PDF) of a BMA probabilistic forecast by assuming that the likelihood function and prior PDF follow normal linear distributions (including the Poisson distribution, gamma distribution, Weibull distribution, and similar distributions). These posterior model probabilities are then adopted as the weights for the models, and the expectation of the BMA model is obtained by weighted averaging of individual model forecasts (i.e., weighted average forecasts).

BMA can be used to combine forecasts and inferences from multiple models (or ensembles of forecast models) and determine the predictive variance of the PDF of corresponding forecast variables. According to studies by Madigan and Raftery (1994) and Hoeting et al. (1999), BMA forecast models can be constructed by extending the law of total probability [see Eq. (1) below], where y represents the forecast variable, e.g., P_min; f_k is the kth forecast; w_k is the weight for the kth forecast, i.e., the optimal posterior probability of forecasts from the corresponding model, which indicates each model’s contribution to forecast skill over the training period and satisfies $\mathop \sum \nolimits_{k = 1}^K {w_i} = 1$ ; and $\tilde {{f_k}}$ is the bias-corrected value of ${f_k},{g_k}(y|{\tilde f_k})$ that denotes the conditional PDF for k forecast models about $\tilde {{f_k}}$ :

$p\left( {y{\rm{|}}{f_1}, \ldots ,{f_k}} \right) = \sum\nolimits_{k = 1}^K {{w_k}} \cdot {g_k}(y|{{\tilde f}_k}).$

(1)

In the BMA model, the predicted value of the forecast variable y is the weighted means of the models’ expectations, which were obtained by weighting the expectations with w_k [see Eq. (2)]. The BMA model’s total variance is expressed as a relationship between the observed values of y_{s, t} and ${\tilde f_{k,s,t}}$ at site s (TC-observed) and timet. It consists of two components: inter-model and intra-model variances [see Eq. (3)] (Raftery, 1993):

$E\left( {y{\rm{|}}{f_1}, \ldots ,{f_k}} \right) = \sum\nolimits_{k = 1}^K {{w_k} \cdot {{\tilde f}_k},\;\;\;\;\;\;\;} \quad\quad\quad\quad\quad\quad $

(2)

$\begin{array}{l}Var\left( {{y_{s,t}}{\rm{|}}{{\tilde f}_{1,s,t}}, \ldots , {{\tilde f}_{K,s,t}}} \right) = \\\displaystyle\sum\nolimits_{k = 1}^K {{w_k}} {\left( {{{\tilde f}_{k,s,t}} - \displaystyle\sum\nolimits_{i = 1}^K {{w_i}} \cdot {{\tilde f}_{i,s,t}}} \right)^2} + \sum\nolimits_{k = 1}^K {{w_k}} \cdot \sigma _k^2.\end{array}$

(3)

Next, we need to solve for w_k and $\sigma $ in the BMA mo-del and iteratively solve for maximum expectations using the maximum likelihood expectation. The $\sigma $ represents all weights and variance. The iterative process includes two steps: calculating expectations and solving the maximization problem. The expectation calculation involves estimating the likelihood function and expectations for the full dataset from a known variable and current parameter estimates. In the second step, the parame-ters are maximized in the full dataset from the first step through re-estimation, provided that the estimation in the first step is correct. The likelihood does not decrease at any iteration, indicating that this algorithm is convergent. This guarantees that it converges to the maximum in the dataset. The operation in the second step is iterated until convergence [see Eq. (4)]:

$l\left( \theta \right) = \mathop \sum \nolimits_{s,t} \log \left( {\mathop \sum \nolimits_{k = 1}^K {w_k} \cdot {g_k}({y_{s,t}}|{{\tilde f}_{k,s,t}})} \right).$

(4)

Assume that there is a latent variable Z_k,s,t whose value is 1 if ${\tilde f_{k,s,t}}$ is the best forecast at site s(TC-observed) and timet; and is zero if ${\tilde f_{k,s,t}}$ represents the worst forecast. An initial guess is given to both w_k and $\sigma $ . The equal weight can be used as the initial guess for w_k, and the empirical $\sigma $ history of the forecast variable can be adopted as the initial guess for $\sigma $ . After $\hat Z_{k,s,t}^j$ is obtained by using Eq. (5), w_k and $\sigma $ can be calculated from $\hat Z_{k,s,t}^j$ with Eqs. (6) and (7), respectively:

$\hat Z_{k,s,t}^j = \frac{{w_k^{j - 1}g({y_{s,t}}|{{\tilde f}_{k,s,t}},\sigma _k^{j - 1})}}{{\mathop \sum \nolimits_{i = 1}^K w_i^{j - 1}g({y_{s,t}}|{{\tilde f}_{i,s,t}},\sigma _i^{j - 1})}},$

(5)

$w_k^j = \frac{1}{n}\mathop \sum \limits_{s,t} \hat Z_{k,s,t}^j,\;\;\;\;\;\;\;\;\;\quad\quad\quad\quad\quad$

(6)

${\sigma ^2}_k^j = \frac{{\mathop \sum \nolimits_{s,t} \hat Z_{k,s,t}^j \cdot {{\left( {{y_{s,t}} - {{\tilde f}_{k,s,t}}} \right)}^2}}}{{\mathop \sum \nolimits_{s,t} \hat Z_{k,s,t}^j}},\quad\;$

(7)

where n is the number of cases in the training set, i.e., the number of distinct values of (s, t).

4.2 Assumption and transformation of PDF

Since the BMA method relies on a normal distribution assumption (Raftery et al., 2005), we first need to ensure that the BMA input variables follow normal distributions. In Fig. 3a, we construct a histogram (with a 10-hPa bin width) with 2246 P_min values obtained from the best-track data for the years 2014 and 2015. As shown in the figure, the best-track observed TC intensity is generally skewed toward the high P_min values. Therefore, we must transform it into a normal distribution. Following Chou et al. (1998), we use the logistics transformation method based on the empirical formulas in Eqs. (8)–(10). The coefficients in the empirical equations are obtained through numerous comparisons and tests. The transformed distribution (see Fig. 3b) appears closer to the normal distribution. This study emphasizes the real-time application of forecasting P_min. If we treat the data from 2014–16 as the full sample, then the real-time forecast will only have the data from 2014–15 available as the training sample. Although we acknowledge that including 2016 data in the training sample would likely improve the results, we chose to perform our analysis in the same way as the real-time operational practice.

First, the observations were ranked in ascending order of P_min as in Eq. (8), and the maximum and minimum values of P_min in this sample were determined. We let y(i) (with a value range of 0–1) denote the original variableP_min [see Eq. (9)]. A function f(y(i)) was then constructed by changing the natural exponential [Eq. (10)]; its distribution is closer to normal than is the historical P_min distribution. For example, when y(1) is equal to 0, f(y(1)) is equal to 0; when y(i) is equal to 0.5, f(y(i)) is equal to 0.38; and when y(n) is equal to 1, f(y(n)) is equal to 1:

$p_{\min }^{\rm obs}\left( i \right)\left| {_{i = 1,n}} \right. \sim \left( {p_{\min }^{\rm obs}\left( 1 \right),p_{\min }^{\rm obs}\left( 2 \right),p_{\min }^{\rm obs}\left( 3 \right), \cdots , p_{\min }^{\rm obs}\left( n \right)} \right),$

(8)

where P_min^obs(i) is sorted from low to high, and

$y\left( i \right) = \frac{{P_{\rm min}^{\rm obs}\left( i \right) - P_{\rm min}^{\rm obs}\left( 1 \right)}}{{P_{\rm min}^{\rm obs}\left( n \right) - P_{\rm min}^{\rm obs}\left( 1 \right)}},\;\;\;\;\;\;$

(9)

$f\left( {y\left( i \right)} \right) = \frac{{{e^{y\left( i \right)}} - 1}}{{e - 1}}. \quad\quad\quad\quad\;\;$

(10)

Figure 3b illustrates the distribution of f(y) that was obtained by matching the P_minvalues in the best-track data for 2014–16 to the 24-h forecasts from the three models after the logistics transformation (sample size: 1265). As can be seen in this figure, high frequencies are largely concentrated around the mean of the f(y) distribution.

A normality test is needed to determine whether the f(y) distribution is normal or approximately normal. The standardized skewness and kurtosis of the distribution were calculated by using Eq. (11) and (12), respectively:

${g_1} = \frac{{\sum\nolimits_i^N {{{\left( {{x_i} - \bar x} \right)}^3}/N} }}{{{s^3}}},\quad\;\,$

(11)

${g_2} = \frac{{\mathop \sum \nolimits_i^N {\raise0.7ex\hbox{${{{\left( {{x_i} - \bar x} \right)}^4}}$} \!/\!\lower0.7ex\hbox{$N$}}}}{{{s^4}}} - 3.$

(12)

In the above equations, $\bar x$ and s denote the sample mean and sample standard deviation, respectively. After the logistics transformation, the skewness changes from –0.31 to 0.15, which is closer to zero. This demonstrates that the transformation changes the originally right-skewed distribution of the sample data to a left-skewed, approximately symmetric one (it is perfectly symmetric when the skewness is zero). After a subtraction of 3 from Eq. (12) (Westfall, 2014), the final kurtosis is closer to zero. The transformation changes the kurtosis of the distribution from –0.78 to –0.59, which suggests that the resulting distribution is closer to a standard normal distribution.

The dimensions of the variables were standardized to facilitate subsequent validation and comparison. After a BMA model was constructed and used to predict f(y), Eqs. (9) and (10) were rearranged to create Eqs. (13) and (14), in which f(y) has the same units as P_min. This step will not be detailed in Section 5:

${\rm{y}}( i ) = \ln \left[ {\left( {e - 1} \right) \cdot f\left( {y\left( i \right)} \right) + 1} \right],\quad\quad\quad\quad\quad\quad\;$

(13)

${P_{\min}}\left( i \right) = y\left( i \right) \cdot \left( {P_{\min}^{\rm obs}\left( i \right) - P_{\min}^{\rm obs}\left( 1 \right)} \right) + P_{\min}^{\rm obs}\left( 1 \right).$

(14)

5 Forecasts from the BMA model and analysis 5.1 BMA model construction

To construct a BMA model, we first conducted a full-sample forecast experiment. All the data for 2014–16 were used in the training period. The TC intensity forecasts for the same period were considered in the validation. Figure 4 shows the experimental forecasts for the No. 1623 TC called “MEARI” (WP26), initialized at 1200 UTC 6 November 2016. Of the three forecast models, HWRF had the highest weight, at 0.55, compared to the GFS’s 0.25. The weight of EC for the raw forecast was 0.2, which is lower than that of the other models. The weights of the three models totaled one.

The full-sample forecasts for this study were obtained from all 2014–16 model forecasts with the bias correction for EC, the transformation of the P_min distribution, and dimensional standardization. Additionally, an equally weighted ensemble (EQW) of the three models was used as a benchmark for comparison. As shown in Table 2 and Fig. 5, the BMA forecasts had the lowest average RMSE during the three years, and the average RMSE of the EQW forecasts was lower than those of the three individual forecasts. The proportions of weights for the members of the BMA model show that HWRF was more skillful than the other members at each lead time. As the lead time increased, the weight for EC-BC increased significantly, while the GFS’s weight tended to decline.

The experimental results demonstrate that the constructed BMA model has the best performance. However, this experimental scheme cannot solve or prevent the occurrence of possible overfitting in future applications. To make the experiment closer to the actual operational TC forecasting approach and achieve full usage of the limited sample, it is necessary to design a sample-matching scheme that does not require coincidence between the training period and the forecast period.

5.2 Cross validation

The main purpose of cross validation is to extract as much information as possible from the limited sample data to make up for the deficiency caused by the limited sample size and to avoid overfitting wherever possible. This study adopted an eight-fold cross validation. The dataset from 2014–16 was divided into eight folds: seven folds for training and the remaining one for validation. The results of this validation are presented in Table 3 and Fig. 6. The sample sizes are slightly different from those in Table 2 to allow for even division of the sample size by eight for the eight-fold cross validation in Table 3. As can be seen in the figure, after cross validation, the RMSEs of the BMA forecasts were the lowest at different lead times and decreased about 5%–7% from those of the forecasts from HWRF, which was the best individual model. The predictive variance ( $\sigma $ ) of the normal distribution of the BMA forecasts was slightly lower than the corresponding RMSE. This further confirms the applicability of BMA to TC intensity forecasting.

In addition, BMA can also be used to quantify the uncertainty in TC intensity forecasts and to provide a probabilistic forecast product. In Table 3, the predictive variance from BMA can be seen to be statistically similar (or equal) to the BMA errors. This indicates that BMA may not only provide the best mean forecasts (those with the smallest errors) but also offer reliable forecast uncertainties (probability).

For example, a forecast was initialized at 0000 UTC 12 September 2016. Figure 7 displays the P_min values predicted by the BMA model using a box plot, which visualizes the uncertainty and dispersion in the forecasts. The figure reveals that the degree of dispersion in the forecasts increased with an increasing length of lead time. In the box plot, the upper and lower bounds of each box indicate the forecast values when the cumulative probability was 75% and 25%, respectively. The black line in the middle of each box indicates the forecast value at a cumulative probability of 50%. The upper and lower ends of each dotted line represent the maximum and minimum P_min forecasts at each lead time.

Compared to deterministic forecasting, the probabili-stic forecasting method based on BMA can provide forecasts with higher accuracy and more comprehensive information. It also outperforms the forecasting method of equally weighted ensemble averaging. Moreover, the BMA model can quantify reliable forecast uncertainty, which is a distinct advantage over other models.

6 Conclusions

This study aims to improve the accuracy of TC intensity forecasts, especially forecasts of the minimum pressure at cyclone center (P_min). Based on Bayesian statistics, the posterior PDF of the BMA probabilistic forecasts was derived by normal approximation to the likelihood function and the prior PDF. Weighted averaging of forecasts from multiple models was implemented. The following conclusions can be drawn from the dataset processing and forecast experiments.

(1) The TC intensity (P_min) forecasts from the HWRF, EC, and GFS models were tested and assessed. The results show that, among the three models, the regional model HWRF has the highest accuracy, due to its relatively high resolution and application of a bogus TC. The P_minvalues predicted by the EC were high overall, indicating low TC intensities. Bias corrections can be applied to the EC model forecast to improve its results, because the model demonstrates consistency in its biases.

(2) In the probabilistic forecast experiment using BMA, bias correction for the EC forecasts was performed. The P_min PDF was approximated by an assumed distribution and then transformed. After that, cross validation was im-plemented to process the data for more efficient use of the sample. The validation results show that the BMA model outperformed the HWRF, which was the best individual model. The RMSE of the BMA forecasts was 5%–7% lower than those of other forecasts at various lead times.

(3) The forecasts from the BMA process provided probabilistic forecast guidance (i.e., mean forecasts and uncertainty forecasts) that could provide more useful information for our future research and for forecasters.

This study could be extended to enhance forecast skill by applying BMA to a pressure field and then projecting an associated wind field based on the pressure–wind relationship.

If we could remove systematic error (or bias) for the NCEP-GFS and NCEP-HWRF, the results could be even better; if the three models exhibit similar forecasting skills, the BMA method may have an even greater advantage in improving the final forecasts.

Acknowledgments. The authors thank NCEP/EMC for the instruction on understanding and applying the EMC’s advanced technology in tropical cyclone forecasting and post-processing of ensemble forecasts, in particular for setting up the BMA algorithm. With their help, the first author was able to carry out an operational forecast experiment at the National Marine Environmental Forecasting Center of China and to implement related forecasting techniques applicable to operational marine forecasting in China. The authors also thank Mr. Eric Sinsky for English grammar editing.