Global Hopf bifurcation of a delayed phytoplankton-zooplankton system considering toxin producing effect and delay dependent coefficient

1.   Introduction

1.1.   Research background

A typhoon is a kind of tropical cyclone. When the sustained wind speed in the center of a tropical cyclone reaches category 12 to 13 (32.7–41.4 meters per second), it is called a typhoon ^[1]. Typhoons usually form on the sea surface in tropical areas 3–5 latitudes away from the equator, such as the North and South Pacific, the North Atlantic, and the Indian Ocean. Its movement is mainly affected by large-scale weather systems, and eventually dissipates at sea, becomes an extratropical cyclone, or dissipates after arriving on land. A typhoon is sudden and destructive; it will bring strong winds and rainstorms wherever it goes. The storm surge disaster caused by typhoons is one of the most serious natural disasters in the world. Sun et al. ^[2] found that typhoons will become stronger, larger, and more destructive in the context of global warming. Kossin ^[3] mentioned in his study that the translation speeds of tropical-cyclones have decreased globally by 10% and decreased by 30% in Asia during 1949–2016. The lower the velocity of a cyclone, the more rain falls on the location, and the higher the risk of flooding and other disasters. Therefore, slower translation speeds will lead to more serious disasters. In order to reduce the economic losses and casualties caused by typhoons, it is urgent to establish accurate and timely prediction methods.

The motion of a typhoon is highly complex. Unlike many time series problems, this cannot be solved by using the characteristics of periodicity ^[4]. A meteorological department mainly uses numerical simulation to achieve typhoon track prediction. This requires a lot of typhoon observational data. In order to establish complex atmospheric equations to simulate the motion and variation of the atmosphere, researchers need to fully understand the internal structural changes, key physical processes, and the mutation mechanism of the atmosphere ^[5]. South Korea, for example, conducts numerical simulations on a Cray XC40 supercomputer with 139, 329 CPUs. Due to its difficulty, numerical simulation requires a lot of calculation time. Even so, numerical simulation cannot obtain the analytical solution, only the numerical solution. Also, expensive hardware means high maintenance costs.

A deep learning method can use a single graphics processing unit (GPU). The prediction results can be generated within seconds by updating the input data of the prediction model. Compared with using hundreds of thousands of CPUs in a simulation, this method is obviously quicker and cheaper. There have been many investigations into the use of deep learning methods for typhoon track prediction. Song et al. ^[6] used the spatial location information and meteorological factors of typhoons to train a neural network which is constructed based on convolutional neural networks (CNNs) and the gated recurrent neural network (GRU). They state that the proposed model is better than the traditional prediction methods for typhoon track prediction. Song et al. ^[7] used BiGRU with an attention mechanism to predict typhoon tracks. Gao et al. ^[8] trained a LSTM with typhoon track data. In addition to describing the relevant physical quantities of the typhoon process, satellite images can also effectively reflect the spatial information.

Satellite images can provide the distribution of clouds, and the changes between satellite images can reflect the evolution of weather systems. With the help of satellite images, more accurate weather forecasts can be made. Lee et al. ^[9] were probably the first to consider using satellite images of tropical cyclones for typhoon track prediction with the help of neural networks. However, the images did not function as input data for the LSTM. Therefore, they extract information such as Dvorak number, maximum wind speed, or spatial location of the cyclone from satellite images and input these into the network for time series prediction. Kovordanyi et al. ^[10] used satellite images as neural network input data for the first time. They predicted the shape of cloud clusters and the direction of movement of the typhoon via a neural network. Rüttgers et al. ^[11] marked the typhoon center on the satellite images and used generative adversarial networks (GANs) ^[12] to generate the satellite images to predict the shape of the cloud cluster and the typhoon track.

The prediction of typhoon tracks is affected by many factors such as the typhoon model structure, terrain, initial field, and boundary conditions of the model ^[13]. For example, the movement of a typhoon is affected by the steering flow, which is determined by the environment around the storm. When a typhoon is affected by multiple weather systems at the same time, it is difficult to predict the steering flow ^[14]. Therefore, a single neural network cannot accurately predict the typhoon track. Therefore, the fusion of multiple neural networks is more suitable for the task of typhoon track prediction rather than a single neural network ^[7]. Xu et al. ^[15] modeled the spatiotemporal relationship. They attached the trained generator to the output of the LSTM to generate the cloud images. Hong et al. ^[16] applied a convolutional LSTM network (ConvLSTM) to typhoon track prediction.

1.2.   Motivation and contributions

At present, research based on deep learning methods combined with satellite images mostly focuses on typhoon eye tracking ^[17] or typhoon intensity prediction ^[18], and less on typhoon track prediction. We believe that it is necessary to extract the temporal and spatial features of time series satellite images for typhoon path prediction. Therefore, we propose a new typhoon track prediction model based on a CNN and LSTM: DeepTyphoon. The specific contributions of this paper are as follows: We mark the satellite images obtained from the Japan Meteorological Agency to create the data set. According to the characteristics of the images in the data set, we establish a prediction network. The features of the clouds and marked squares in the image are extracted by dilated convolution. In addition, a convolutional block attention module (CBAM) is used to enhance the feature extraction ability of the network from two dimensions: Channel and spatial.

2.   Materials and method

2.1.   Related work

2.1.1.   Convolutional block attention module

The attention module can focus on the local information of the network and suppress invalid information, so as to enhance the expression of effective features. The attention mechanism was initially applied to natural language processing and achieved good results. However, there are few studies on the application of attention mechanisms to time series prediction in meteorology. The convolutional block attention module is a new convolutional attention module proposed by Sanghyun et al. in 2018 ^[19]. The convolution operation extracts information features by mixing cross channel information and spatial information. Therefore, they innovatively proposed the attention mechanism of integrating channel attention and spatial attention and applied a channel attention module and a spatial attention module for the two respective dimensions. The principle of CBAM is shown in Figure 1.

In this paper, CBAM is applied to the encoder of the neural network to enhance the feature extraction ability of the convolution network, and extracts the features of marked squares and clouds from the two dimensions.

2.1.2.   Long short-term memory network

The long short-term memory network is a special recurrent neural network. It was first proposed by Hochreiter and Schmidhuber ^[20] and later improved by Grave ^[21], which can address the problems of gradient disappearance and gradient explosion in the process of long sequence training. Like RNNs, LSTMs can also connect previous information to the current task, such as using the previous video frame information to understand the current frame. Compared with RNNs, LSTMs perform better in learning long-term dependencies.

The operation of LSTM is mainly realized by the input gate i_t, forget gate f_t, and output gate o_t. First, the LSTM uses the forget gate f_t to decide which information in c_t-1 to delete, and the activation state of the forget gate is determined by the activation function σ(·):

where W_f is the weight matrix of the forget gate; b_f is the offset term of the forget gate; f_t represents the forget gate at time t; x_t is the input of network at time t. Second, the LSTM uses the input gate i_t to determine the information stored in the new memory unit c_t:

where W_i is the weight matrix of the input gate; W_c is the weight moment of the cell state; b_i is the offset term of the input gate; b_c is the bias term of the cell state; i_t is the input gate at time t; tanh is the hyperbolic tangent function; $\widetilde {{c_t}}$ indicates the candidate value c_t added to the new memory unit status. Then, the LSTM updates the status of the previous memory unit c_t-1 to the status c_t:

where c_t is the input status of the memory unit. Finally, the output gate o_t is used to calculate the hidden state h_t. The specific formula is as follows:

where W_o is the weight matrix of the output gate; b_o is the output gate offset term; σ is the sigmoid function; $\odot$ represents the multiplication of matrix elements.

2.1.3.   Hybrid dilated convolution

The dilated residual network ^[22] was proposed to improve the receptive field of convolutions. It maintains the width and height of input features by setting dilation rates. However, the setting of dilation rates leads to the "gridding issue"; when the dilation rate increases, the input sampling will become very sparse, which will be detrimental to the learning of features. Therefore, Wang et al. ^[23] proposed hybrid dilated convolution (HDC) which alleviates the gridding issue by setting the dilation rates for a series of convolution layers and connecting them. In this paper, the size of the marked square on the satellite images is 3 × 3 pixels. As shown in Figure 2, the three-layer HDC layer with a dilation rate of 1/2/3 can ensure that the input feature of the square is not lost.

2.2.   Establishment of the data set

2.2.1.   Data acquisition and data cleaning

The self-built typhoon image data set used in this paper is provided by the National Institute of Information (NII). The satellite image data are taken from Himawari-1~8 series and GOES meteorological satellites. Of these, Himawari-8 can achieve a 10-minute time resolution and 500 m spatial resolution. It has an improved imager observing five bands (one for visible (0.64 μm) and four for infrared (10.4, 12.4, 6.2 and 3.9 μm)). GOES is the geostationary orbit operational satellite series of NOAA in the United States. It also has an improved imager observing five bands (one for visible (0.55–0.75 μm) and four for infrared (10.20–11.20, 11.50–12.50, 6.50–7.00 and 3.80–400 μm)). Satellite images in the recent 40 years are selected as data samples, which are all infrared images (10.4 μm) with a size of 128 × 128 pixels, representing an actual range of 2600 × 2600 km. The middle of the image is the typhoon center, and NII provides the longitude and latitude of the typhoon position corresponding to each image.

Due to the different life cycles of typhoons, the typhoon processes with typhoon cycle less than 60 hours are screened out through data cleaning in order to ensure that the neural network can obtain sufficient information in the training process. We selected satellite images located in the Northwest Pacific and Southwest Pacific to demonstrate the effectiveness of the model. Of these, 513 typhoon processes located in the Northwest Pacific from 1980 to 2020, that is, 513 typhoon sequences, and a total of 13400 satellite images are selected to construct the training set, verification set, and test set, which are 338, 103 and 72 typhoon sequences, respectively, corresponding to 9110, 2680 and 1610 satellite images. Also, 198 typhoon processes located in the Southwest Pacific from 2000 to 2020, a total of 8420 satellite images, are selected to construct the training set, verification set, and test set, which are 108, 60 and 30 typhoon sequences, respectively, corresponding to 5470, 2110 and 840 satellite images.

The original infrared image is a single channel gray image, and the range of digital number of satellite image is 0–255. Typhoons are sequential from start to finish. Moreover, the same typhoon has different forms in different lifecycle stages, and the spiral cloud system shown on the satellite cloud map is also different ^[24]. Figure 3 shows some images in the lifecycle of typhoon Damrey, which lasted for 11 days and 18 hours.

2.2.2.   Marking the satellite images

As described in Section 2.2.1, the typhoon center in the satellite image provided by NII is located in the middle of the image. However, the typhoon center on the satellite images is an absolute position. We need the relative positional information between satellite image sequences to judge the movement of clouds. In order to determine the direction of movement and distance of the typhoon through the satellite image, for each satellite image, we mark the typhoon position at the previous time (for example, six hours ago) on the satellite image at the current time. The marking method is as follows: The spatial coordinate of the typhoon at the current time is (φ, λ). The horizontal offset D_φ and vertical offset D_λ (in km) between the coordinates of the previous time and the coordinates of the current time is calculated according to the Havers formula ^[25]; the earth radius R is taken from the real coordinate system:

Then, according to D_φ and D_λ, we calculate the horizontal and vertical offset P_φ and P_λ at the pixel level on the image using the reference system. This determines the position of the typhoon at the previous time on the satellite image at the current time.

We use black squares to mark satellite images (gray image) for the experiments. Compared with using RGB image (image data marked with red squares) as the input of neural network, the use of gray images as the input of neural network means a shorter training time which is convenient for the experiments to adjust the parameters of the neural network. In the experimentation, it is found that using a red square marked image can eliminate the interference of the gray background, so as to determine the position of the typhoon more accurately. As shown in Figure 4: Figure 4(a) is the original satellite image, Figure 4(b) is the satellite image marked with a black square, and Figure 4(c) is the satellite image marked with a red square. In Figure 4(b), (c), the center of the image is the typhoon center at the current time, and the marked square represents the typhoon center at the previous time. The size of the marked square is 3 × 3 pixels. There are 13400 single channel gray images in the typhoon image data set with black square marks, and 13, 400 RGB three channel images in the typhoon image data set with red square marks. The range of digital number of image data marked with square (black or red) are both 0–255.The size of the dataset does not change as a result of marking. The data set size of the Northwest Pacific is 9110 images for the training set, 2680 images for the verification set, and 1610 images for the test set. For the Southwest Pacific data set, this is 5470, 2110 and 840, respectively.

2.3.   DeepTyphoon prediction model

2.3.1.   Design of DeepTyphoon

In the research of using deep learning to predict typhoon tracks, most researchers use reanalysis data as the input of neural network. For example, the work in ^[26] converts the reanalysis data into a density map and uses it as the input of a neural network to determine the spatial position of the typhoon according to the predicted density map. However, compared with the thermodynamic diagram transformed according to physical quantity data, satellite images contain more intuitive information. In order to consider the temporal and spatial characteristics of the typhoon process at the same time, we propose the DeepTyphoon prediction model for satellite image training and prediction. DeepTyphoon combines the ability of a CNN to extract advanced semantic features and the utility of LSTM for time series prediction. It takes the marked satellite image as the input, learns the relative position between the square mark and the cloud cluster and the direction of movement of the cloud cluster in the training process, and finally predicts the position of the typhoon according to the square mark on the predicted image.

2.3.2.   Prediction process of DeepTyphoon

The prediction process of DeepTyphoon is shown in Figure 5. The whole neural network model includes a convolution neural network (as encoder), an LSTM, and a deconvolution neural network (as decoder). Taking the satellite images marked with black squares as an example, the prediction process is as follows: First, the marked satellite images (denoted by m) are input into the encoder. We use the hybrid dilated convolution layer to extract the features of the satellite image m. At the same time, we add the channel attention module and the spatial attention module to further improve the feature extraction ability of the coding network. After each image of a satellite image sequence is input into the encoder, the encoder outputs the feature vector of m (denoted by t) as the input of LSTM. The output of the LSTM network t' is the predicted value of t. Finally, we deconvolute t' into an image m' in order to determine the typhoon center according to the prediction result t' of LSTM. Due to the black square mark in the input image m, the final output of the model, the predicted image m', will also have a predicted black square mark.

The black square in the predicted image is detected by applying a filter on the predicted image and a convolution with the size of the black square on the filtered image. The pixel with the highest value in the convolved image gives the location of the predicted typhoon center. After identifying the pixel coordinates of the predicted black square, they are converted back to coordinates of latitude (φ) and longitude (λ) through the geographic reference system.

2.4.   Experiment

2.4.1.   Experimental environment

The experimental environment of this paper is mainly based on the PyTorch framework and built on the Ubuntu operating system. The specific configuration is shown in Table 1 below. In order to fully verify the performance, the experiments in this paper are completed in the same experimental environment. The experimental hyperparameters are set as follows: The batch size is 100, the learning rate is 0.0001, and the number of epochs is 100. The Adam algorithm is used to optimize the network parameters, and the weight attenuation coefficient is set to 0.00002.

2.4.2.   Optimization of DeepTyphoon

In order to test the influence of the number of convolution layers on the experimental results, we set different convolution layers for testing. First, we set the dilation rates of the three continuous convolution layers to 1/2/3, respectively, and take the three convolution layers as one HDC layer. We use the satellite image marked with the red square as the input of the model, and calculate the mean average error (MAE) of the predicted results when one, two, and three HDC layers are used.

As shown in Table 2, we achieved the minimum MAE when using two HDC layers. That is to say, when we use six convolution layers and set the dilation rates to 1/2/3/1/2/3 respectively, the prediction result of the neural network is the best.

An important hyperparameter of the model is the number of input image sequences. If the number of input image sequences (denoted by n) is too large, the memory usage will be increased, and the training time will be increased. Moreover, the lifecycle of a typhoon is not too long, and too large a value of n will make most typhoon processes unable to be used as training data. However, too small a value of n will lead to insufficient input information and affect the prediction performance. We use the satellite images marked with red squares as the input. The time interval of the satellite image is 6 hours. After training and testing, for the three cases where n is 5, 10 and 20, the MAE after 6 hours is calculated.

As shown in Table 3: In the case of n = 5, the change of the typhoon seems to have not been fully captured. In the case of n = 20, due to the reduction of the number of training data sets, the results are not good compared with n = 10. For example, considering typhoon Mitag with 29 sequences, using n = 10 allows to start training from sequence 11, and provides a number of 19 training sequences. However, when n = 20, this only allows to start training from sequence 21, that is, 9 training sequences.

During the experiment, it is found that different convolution kernel sizes will affect the shape of the final generated square. In order to accurately predict the square in the predicted images, the convolution kernel size is tested. We select the satellite image marked with a black square as the input, which can save calculation time.

Figure 6 shows the prediction results corresponding to different convolution kernels. Figure 6(a) shows the prediction results when the convolution kernel size is 3. It is obvious from the prediction image that a 3 × 3 pixel black square is generated. Figure 6(b), (c) correspond to the shape of the predicted square in the generated image when the convolution kernel size is 2 and 4, respectively. Although it can be seen from the predicted image that the marks are generated, they are blurred to varying degrees and the central position cannot be accurately determined by convolution. Finally, the size of the convolution kernel is set to 3 so that the exact position of the square can be obtained.

For the typhoon satellite image at the later stage of the typhoon lifecycle, as shown in Figure 7(a), if the black square is used to mark the satellite image, the square in the predicted image almost coincides with the gray background of the image. It is impossible to accurately detect the mark position in the prediction image through convolution. Therefore, we choose the satellite image marked by the red square as the model input so that we can accurately detect the location of the marked pixels. The channel attention module is added after the convolution layer to enhance the weight of the pixel value of the image R channel (red). In the predicted image, only the pixel value on the R channel in the generated image is detected, and the position of the marked square can be obtained accurately. As shown in Figure 7(b), the red square in the figure is located and redrawn after convolution on the R channel.

The parameters of the neural network are finally determined through experiments, as shown in Table 4. Of these, the number of HDC layers is obtained from Table 2, which is six convolution layers with different dilation rates. In order to decode the prediction vector into an image, the number of deconvolution layers is the same as that of convolution layers. The kernel size is determined according to Figure 6. Based on previous experience, we set the number of hidden layers of LSTM to two.

2.4.3.   Evaluation indicators

In order to evaluate the performance of the prediction model, the root mean square error (RMSE) and the mean absolute error are selected to measure the prediction accuracy of the model.

The absolute error (E) describes the distance between the predicted coordinates (φ_pred, λ_pred) and the real coordinates (φ_real, λ_real). It is calculated in kilometers (km) by applying the haversine formula, and the earth radius R is taken from the real coordinate position:

The mean absolute error formula is as follows, in km, where E represents the absolute error of the typhoon:

The root mean square error formula is as follows, in km, where E represents the absolute error of the typhoon:

3.   Results

3.1.   Experimental results and analysis

As an example of a predicted image, Figure 8 shows the prediction results of typhoon Gray sequence 11, in which Figure 8(a) is the ground truth and Figure 8(b) is a predicted image.

It is obvious that the predicted image does not have the same definition as the real image; in particular, the clouds are blurred compared with the ground truth. However, we focus on accurately predicting the position of the square in the predicted images, so as to accurately predict the spatial position of the typhoon rather than improve the quality of the generated images.

By calculating the relative position between the square and the center of the predicted image, the spatial position of the predicted typhoon is determined. We compare the position of the predicted typhoon with the position of the real typhoon, and finally we can obtain the prediction accuracy. We select an example in the Northwest Pacific and the Southwest Pacific to show the prediction results of DeepTyphoon. Each typhoon is predicted in chronological order, and the sequence is named by the date and time of the corresponding real image UTC. For example, "1993072718" represents 6 p.m. on July 27, 1993 (UTC).

Table 5 shows the coordinate prediction results and prediction errors of Typhoon Mitag. As can be seen from the table, this typhoon has taken 29 steps from occurrence to extinction; each time step is 6 hours, giving a total of 174 hours. Take ten points as a sequence as the input of the prediction model, the predicted values of Mitag sequence 11 to sequence 29 can be obtained in order. After calculation, the average error of Typhoon Mitag from sequence 11 to sequence 29 is 55.037 km.

As shown in Figure 9, in order to display the forecast results more intuitively, they are marked in the map according to the real track (yellow) and predicted track (blue) of Typhoon Mitag.

Table 6 shows the coordinate prediction results and prediction errors of Typhoon Marcia. Typhoon Marcia has taken 22 steps from occurrence to extinction; each time step is 6 hours, giving a total of 132 hours. Take ten points as a sequence as the input of the prediction model, the predicted values of Marcia sequence 11 to sequence 22 can be obtained in order.

As shown in Figure 10, we also display the forecast results for Typhoon Marcia in the map according to Table 6. The real track is represented by yellow points, and the predicted path is represented by blue points.

We tested two test sets located in the Northwest Pacific and Southwest Pacific. Table 7 shows the prediction results of DeepTyphoon for the two test sets. The predicted time range is 6 hours. The average absolute error of the prediction results of DeepTyphoon for the two test sets is less than 70 km. In addition, we recorded the calculation time of the DeepTyphoon prediction model. After 14 hours of training, the DeepTyphoon prediction model can process 1000 satellite images, that is, 100 typhoon points, in one minute.

3.2.   Comparison with other models

Table 8 shows the prediction results for RNN, GAN, AE–GRU, TraijGRU, and DeepTyphoon. The performance indicators include RMSE and MAE.

As shown in Table 8, the evaluation indexes for DeepTyphoon are the best after 6 and 12 hours. Through comparison, it is shown that DeepTyphoon can more accurately predict the location of typhoons using satellite images.

In order to further verify the validity of the proposed model, we calculated the MAEs in longitude and latitude and compared it with SOTA models. As shown in Table 9, our model has the best prediction result at 6 hours and 12 hours. This verifies the validity of the proposed model.

4.   Discussion

The experimental results show that marking the typhoon center at the previous time on the satellite images and using CNN and LSTM for feature extraction can effectively learn cloud motions in the satellite image. Also, we can use the deconvolution network to detect the typhoon center at a future time so as to accurately predict the typhoon path in the future.

However, we expect to get more information from the prediction image obtained by the deconvolution network. Therefore, the focus of future work should be to further optimize the neural network and better extract cloud features, so as to obtain more realistic satellite images through prediction. In this way, we can not only detect the typhoon center in the prediction image, but also determine the intensity of the typhoon. In addition, the prediction accuracy of the typhoon path should be further improved. Then, how to ensure the accuracy of typhoon track prediction and the authenticity of the prediction image at the same time is the difficulty in the future.

5.   Conclusions

In this paper, we use a deep learning method combined with satellite images to predict typhoon tracks, and propose the DeepTyphoon prediction model. The typhoon center is marked on the satellite image. Combined with the characteristics of CNN and LSTM, we extract the features of the satellite image by HDC and CBAM, and use LSTM to predict the features. After the predicted feature vector is transformed into an image, the position of the typhoon is detected on the predicted image. More than 20000 time series satellite images in two data sets released by the Japan Meteorological Agency in recent years are used for training and prediction. Experimental results on the Northwest Pacific and Southwest Pacific data sets show that our prediction method is better than existing deep learning methods in typhoon track prediction.

Acknowledgements

This work was supported by "Shanghai science and technology innovation plan" (Project number: 20dz1203800) and "the Capacity Development for Local College Project" (Project number: 19050502100).

Conflict of interest

The authors declare there is no conflict of interests.