# Extracting Modulated Annual Cycle in Climate and Ocean Time Series Using an Enhanced Harmonic Analysis

Authors: Yibo Zhang; Haidong Pan; Shuang Li; Xianqing Lv
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Accurate extraction of the modulated annual cycle (MAC) is important for climatic and oceanic research. A variety of methods are available to extract the annual cycle with inconsistent results. Since actual annual cycles are unknown in the observation series, the reliability and applicability of the results extracted by these methods are difficult to estimate. In this study, three widely used decomposition methods, ensemble empirical mode decomposition (EEMD), nonlinear mode decomposition (NMD), and enhanced harmonic analysis (EHA), are evaluated by idealized numerical experiments for extracting modulated annual cycles from climate series. Idealized numerical experiments are carried out and show that the recently proposed EHA had the most accuracy in extracting the MAC from the constructed data. The optimal independent point (IP) number, which makes the most accurate result for EHA, can be found in each ideal experiment. In the actual experiment, two IP selection criteria are proposed for EHA to extract MAC from observations.


---

## Body

## 1. Introduction

The annual cycle is an important aspect of climatic and oceanic data research. Air temperature, water temperature, water level, and other data have significant annual cycles. For example, the annual cycle can explain 60–90% of monthly sea level variance in the marginal northwest Pacific [1]. In the past, the annual cycle was usually regarded as an exact repeat from year to year [2]. This stationary annual cycle, known as the traditional annual cycle (TAC), is often used in anomalies calculation [1–3]. Anomaly is usually defined as a deviation of a value from its typical position in the annual cycle. Numerous papers on climate science have been published, which use monthly anomalies and propose a physical mechanism to explain it [2]. However, TAC ignored the variability of amplitudes and phases of the annual cycle and mixed the changes in the annual cycle into the anomalies. In order to explain the annual cycle changes which are usually caused by changes within the nonlinear climate system, a description of the annual cycle needed to be improved, including the amplitude changes from one year to the next, or rather a modulated annual cycle (MAC) [2]. How and why the MAC changes have attracted increasingly more attention and various methods have been used to extract the MAC [4].The TAC of climatic and oceanic data is generally considered to be expressed as a trigonometric function, and the amplitude and phase of the TAC can be obtained through fitting models, which is the elementary method to extract TAC. Feng et al. [1] extracted the MAC with moving window fitting models. Although this method can obtain the continuously time-varying amplitude and phase of the annual cycle, it usually underestimates the amplitude change. Wu et al. [2] used ensemble empirical mode decomposition (EEMD) [5] to extract MAC. Qian et al. [3] used EEMD to analyze the Niño-3 SST index from 1880 to 2008 and found that its MAC has strong interannual modulation and long-term changes in both amplitude and phase. Deng et al. [6] used nonlinear mode decomposition (NMD) to extract the change of MAC and found that the interannual change of the annual cycle ignored by TAC would affect the quantification of long-term persistence (LTP).To test the effectiveness of various methods, Deng and Fu [7] compared the extraction results of five methods, including the fitting models, complex demodulation (CD), EEMD, NMD, and seasonal trend decomposition procedure based on loess (STL), through idealized numerical experiments involving MAC with three different time-varying amplitudes, and they concluded that NMD has the optimal performance. In addition, NMD is a more complicated method, which requires more resources on computation.Jin et al. [8] proposed an enhanced harmonic analysis (EHA) to analyze temporal variations of internal tidal currents in the South China Sea. Based on the independent point (IP) scheme [9], EHA can generate time-varying amplitudes and phases for specific frequencies. Pan et al. [10] developed a MATLAB toolbox of S_TIDE based on the EHA and applied it to explore tide-river interaction in the Columbia River Estuary. Recently, EHA has been applied to investigate the seasonal variation of the principal tidal constituents in the Bohai Sea [11] and the changes of nodal cycles in the Gulf of Maine [12]. In this paper, we firstly apply EHA to the extraction of the MAC and compare it with EEMD and NMD.The paper is organized as follows. In Section2, we introduce the extraction methods: EEMD, NMD, and EHA. Section 3 describes the design of idealized numerical experiments. Experimental results are discussed inSection 4. Finally, the conclusions are drawn in Section 5.


## 2. Methods

### 2.1. EEMD

Empirical mode decomposition (EMD) is an adaptive signal decomposition algorithm for nonlinear, nonstationary signals proposed by Huang et al. [13]. Through EMD, a complex signal in the time domain can be decomposed into a finite number of components, which are called intrinsic mode functions (IMFs). As an intrinsic mode of a signal, each IMF can be modulated by both amplitude and frequency, meaning that its amplitude and frequency change with time. A time series of climate element data can be decomposed with the EMD method in the following form:(1)St=∑i=1nηit+rt,where n is the number of IMFs, which is determined by factors like the length and variation of the data, end effects, and stoppage criteria, ηit is the i-th IMF, and rt represents the residual trend. Compared with the traditional Fourier transform-based signal analysis method, such as wavelet transform, EMD breaks the limitation of the preselection of wavelet function and has good temporal resolution and adaptability.While the EMD has these wonderful properties, it also has a serious “mode mixing” problem, which is defined as any IMF consisting of oscillations of dramatically disparate scales. To solve the mode mixing problem, the Ensemble EMD (EEMD) was developed [5]. In this method, multiple noises are added to an observed time series to simulate the scenario of multiple trials for a time series, in order to carry out an ensemble average approach for corresponding IMFs and extract scale-consistent signals. The main steps of the EEMD method are as follows [2]:(1)
Add a white noise series to the targeted data and decompose the data with added white noise into IMFs.(2)
Repeat step 1 several times, but the white noise series is different in each time.(3)
Obtain the (ensemble) means of corresponding IMFs of the decompositions as the final result.The effect of EEMD is that the added white noise series cancel each other, and the mean IMFs stay within the natural dyadic filter windows, significantly reducing the chance of mode mixing. EEMD is widely used in many fields, such as MAC [2–4], seismic signal [14], and machinery diagnosis [15].


### 2.2. NMD

Nonlinear mode decomposition (NMD) was proposed by Iatsenko et al. [16] for nonlinear and nonstationary data. This method can decompose a time series into a set of physically meaningful oscillations in the frequency space, based on a combination of time-frequency analysis techniques and surrogate data tests. A given signal St can be decomposed by the NMD method into nonlinear modes (NMs) cit and some additive noise ζt as follows:(2)St=∑i=1ncit+ζt.The key of the NMD method is to obtain the time-frequency representation (TFR) of the signal. Theoretically, all traditional TFR techniques can be used, like wavelet transform, windowed Fourier transform, and so on. Based on the TFR, all possible harmonics of each NM can be reconstructed. Then, using surrogate data sets’ tests to eliminate false harmonics, all remaining harmonics are summed to obtain a complete NM. As a new decomposition method, NMD has been applied to many fields, such as city engineering modeling [17], analysis of surface air temperature [6], and photoplethysmography [18].


### 2.3. EHA

In the traditional fitting models, the amplitude and phase of the annual cycle can be estimated by the least-squares methods, and this method is expressed as follows:(3)St=S0+Hcos2π365t−g,where St is the observations, S0 is a constant, H is the amplitude, and g is the phase. Equation (3) can be rewritten as follows:(4)St=S0+acos2π365t+bsin2π365t,where(5)H=a2+b2,g=Arctanba.EHA assumes the parametersa and b to be time-varying functions. Then, equation (4) can be modified as(6)St=S0+atcos2π365t+btsin2π365t.The time-varying parameters in equation (6) are solved with the help of an independent point (IP) scheme. The IP scheme assumes that the coefficients a and b are changing and selects k representative time moments as IP, which are evenly distributed in the time domain. The coefficients a and b at IPs are called independent parameters (ai and bi), and the values at the other times are obtained by cubic spline interpolation between two adjacent IPs. k is called the IP number. Therefore, the coefficients a and b can be expressed as a function of ai and bi:(7)at=∑i=1kft,iai,bt=∑i=1kft,ibi,where ft,i is the interpolation weight coefficient of the i-th independent point at a time t and the coefficient is determined by the interpolation method. Both ai and bi are unknown at this time. Substituting equation (7) into equation (6), we get(8)St=S0+∑i=1kft,iaicos2π365t+∑i=1kft,ibisin2π365t.AtN times t=t1,t2,t3,…,tN, we can get S=St1,St2,St3,…,StN, so N equations can be obtained:(9)S0+∑i=1kft1,iaicos2π365t1+∑i=1kft1,ibisin2π365t1=St1,S0+∑i=1kft2,iaicos2π365t2+∑i=1kft2,ibisin2π365t2=St2,⋮S0+∑i=1kftN,iaicos2π365tN+∑i=1kftN,ibisin2π365tN=StN.There are2k+1 unknowns in the system of equation (9). When N is much larger than 2k+1, S0, ai, and bi can be solved by the least-square fitting. Then, the time-varying amplitude and phase can be calculated from ai and bi after interpolation. The selection of IP numbers is discussed in detail in Section 4.3. It should be noted that when IP number is one, the amplitude and phase obtained by EHA are constant and the same as the results of the traditional fitting model. This means that the traditional fitting model can be seen as a special case of EHA.


## 2.1. EEMD

Empirical mode decomposition (EMD) is an adaptive signal decomposition algorithm for nonlinear, nonstationary signals proposed by Huang et al. [13]. Through EMD, a complex signal in the time domain can be decomposed into a finite number of components, which are called intrinsic mode functions (IMFs). As an intrinsic mode of a signal, each IMF can be modulated by both amplitude and frequency, meaning that its amplitude and frequency change with time. A time series of climate element data can be decomposed with the EMD method in the following form:(1)St=∑i=1nηit+rt,where n is the number of IMFs, which is determined by factors like the length and variation of the data, end effects, and stoppage criteria, ηit is the i-th IMF, and rt represents the residual trend. Compared with the traditional Fourier transform-based signal analysis method, such as wavelet transform, EMD breaks the limitation of the preselection of wavelet function and has good temporal resolution and adaptability.While the EMD has these wonderful properties, it also has a serious “mode mixing” problem, which is defined as any IMF consisting of oscillations of dramatically disparate scales. To solve the mode mixing problem, the Ensemble EMD (EEMD) was developed [5]. In this method, multiple noises are added to an observed time series to simulate the scenario of multiple trials for a time series, in order to carry out an ensemble average approach for corresponding IMFs and extract scale-consistent signals. The main steps of the EEMD method are as follows [2]:(1)
Add a white noise series to the targeted data and decompose the data with added white noise into IMFs.(2)
Repeat step 1 several times, but the white noise series is different in each time.(3)
Obtain the (ensemble) means of corresponding IMFs of the decompositions as the final result.The effect of EEMD is that the added white noise series cancel each other, and the mean IMFs stay within the natural dyadic filter windows, significantly reducing the chance of mode mixing. EEMD is widely used in many fields, such as MAC [2–4], seismic signal [14], and machinery diagnosis [15].


## 2.2. NMD

Nonlinear mode decomposition (NMD) was proposed by Iatsenko et al. [16] for nonlinear and nonstationary data. This method can decompose a time series into a set of physically meaningful oscillations in the frequency space, based on a combination of time-frequency analysis techniques and surrogate data tests. A given signal St can be decomposed by the NMD method into nonlinear modes (NMs) cit and some additive noise ζt as follows:(2)St=∑i=1ncit+ζt.The key of the NMD method is to obtain the time-frequency representation (TFR) of the signal. Theoretically, all traditional TFR techniques can be used, like wavelet transform, windowed Fourier transform, and so on. Based on the TFR, all possible harmonics of each NM can be reconstructed. Then, using surrogate data sets’ tests to eliminate false harmonics, all remaining harmonics are summed to obtain a complete NM. As a new decomposition method, NMD has been applied to many fields, such as city engineering modeling [17], analysis of surface air temperature [6], and photoplethysmography [18].


## 2.3. EHA

In the traditional fitting models, the amplitude and phase of the annual cycle can be estimated by the least-squares methods, and this method is expressed as follows:(3)St=S0+Hcos2π365t−g,where St is the observations, S0 is a constant, H is the amplitude, and g is the phase. Equation (3) can be rewritten as follows:(4)St=S0+acos2π365t+bsin2π365t,where(5)H=a2+b2,g=Arctanba.EHA assumes the parametersa and b to be time-varying functions. Then, equation (4) can be modified as(6)St=S0+atcos2π365t+btsin2π365t.The time-varying parameters in equation (6) are solved with the help of an independent point (IP) scheme. The IP scheme assumes that the coefficients a and b are changing and selects k representative time moments as IP, which are evenly distributed in the time domain. The coefficients a and b at IPs are called independent parameters (ai and bi), and the values at the other times are obtained by cubic spline interpolation between two adjacent IPs. k is called the IP number. Therefore, the coefficients a and b can be expressed as a function of ai and bi:(7)at=∑i=1kft,iai,bt=∑i=1kft,ibi,where ft,i is the interpolation weight coefficient of the i-th independent point at a time t and the coefficient is determined by the interpolation method. Both ai and bi are unknown at this time. Substituting equation (7) into equation (6), we get(8)St=S0+∑i=1kft,iaicos2π365t+∑i=1kft,ibisin2π365t.AtN times t=t1,t2,t3,…,tN, we can get S=St1,St2,St3,…,StN, so N equations can be obtained:(9)S0+∑i=1kft1,iaicos2π365t1+∑i=1kft1,ibisin2π365t1=St1,S0+∑i=1kft2,iaicos2π365t2+∑i=1kft2,ibisin2π365t2=St2,⋮S0+∑i=1kftN,iaicos2π365tN+∑i=1kftN,ibisin2π365tN=StN.There are2k+1 unknowns in the system of equation (9). When N is much larger than 2k+1, S0, ai, and bi can be solved by the least-square fitting. Then, the time-varying amplitude and phase can be calculated from ai and bi after interpolation. The selection of IP numbers is discussed in detail in Section 4.3. It should be noted that when IP number is one, the amplitude and phase obtained by EHA are constant and the same as the results of the traditional fitting model. This means that the traditional fitting model can be seen as a special case of EHA.


## 3. Ideal Experiments

In this work, we design seven experiments, each of which has a different amplitude change in the annual cycle. Since the amplitude change is more apparent than the phase change, we mainly consider amplitude change with fixed frequency and phase in this work. Nevertheless, in the section of the discussion, we also conduct some experiments where both amplitudes and phases are time-varying. The constructed data are obtained by adding noise to the preset annual cycles. And seven kinds of amplitude changes are set in this work: constant, linear increase, periodic change, and four kinds of complex nonlinear changes. In each experiment, noise with 1%, 2%, and 4% noise-power ratio (NPR) whose Hurst index isα=0.65 is added, respectively, to test the influence of noise on the extraction results. The noise is generated by using a Fourier-filtering method (FFM) [19]. Then, three methods, EMD, NMD, and EHA, are used to extract the annual cycle. By comparing the extraction results with the preset annual cycle, we examine the effectiveness of EHA and provide some suggestions in the actual application.Climatic data can be decomposed into long-term trends, annual cycles, and daily high-frequency varying fluctuations [2]. Like Deng and Fu [7], in order to examine the accuracy of the extraction results of various methods, the constructed data was used instead of the observed data, and the influence of the long-term trend was not considered. The original data consisted of a MAC and high-frequency fluctuations on a time scale of 50 years (one number per day).In this paper, the trigonometric function was used to construct the annual cycleSt with fixed frequency f and phase and changing amplitude At. St can be expressed as the following form:(10)St=Atcos2πft,t=1,2,…,N,where the frequency is f=2.74×10−3d−1. In this paper, 7 kinds of amplitude changes were designed: constant, linear increase, periodic change, and 4 kinds of nonlinear changes (Figure 1). The 7 kinds of amplitude changes At all ranged from 9 to 11. Similar to Deng and Fu [7], noise was considered in this paper to represent high-frequency fluctuations. Noise with a noise-power ratio of 1%, 2%, and 4% is added in each experiment, respectively.Figure 1
Preset annual cycle amplitude changes.


## 4. Result and Discussion

### 4.1. Result

The previously constructed data has been processed by three methods, EEMD, NMD, and EHA, to extract different MACs. This section mainly discusses the deviation between the extraction results and the preset amplitude-modulated annual cycle. In order to avoid the edge error, only 46 years results were selected, excluding the first two years and the last two years. The deviation can be expressed by the root mean square error (RMSE) as follows:(11)RMSE=∑i=1NSi−Si′2N,where Si represents the extracted annual cycle and Si′ represents the preset annual cycle. N is the length of Si. The RMSE in each experiment is shown in Figure 2. The RMSE of the three methods increases with noise power. The RMSEs of EEMD are the largest and an order of magnitude higher than the RMSEs of other methods. When the changes of the preset annual cycle are relatively simple (the first three types), the RMSE of EHA is an order of magnitude smaller than that of NMD. When the preset annual cycle changes are complex, RMSEs of EHA and NMD are in the same order of magnitude. In each group of experiments, the RMSE of EHA is the smallest, EEMD is the largest, and NMD is between them. The RMSEs of NMD and EHA are also affected by the complexity of amplitude variation. At the same noise-power-ratio level, the more complex the amplitude variation is, the larger the RMSE will be. However, the RMSE of EEMD is almost only affected by the noise-power ratio and is not sensitive to the complexity of amplitude changes. EHA can directly obtain the instantaneous time-varying amplitude. For the extraction results of NMD and EEMD, the time-varying amplitudes of the MAC are obtained through a cubic spline fitting to the local maxima of MAC2 series [4]. The result of EHA is closest to the preset amplitude (Figure 3).Figure 2
The RMSE of extraction results with variable amplitude annual cycle.Figure 3
The amplitude changes obtained by the three methods in experiment 6 with 1% white noise.Among the three methods, both EEMD and NMD are adaptive, while EHA, which is subjective, needs users to tune the IP number. When the IP number is too small, the result obtained by EHA is relatively simple and can only reflect low-frequency changes of MAC. The more the IP is, the more the MAC changes can be reflected. When the IP number is too large, the result of EHA is likely to be severely abnormal. This means that we not only need a large enough IP number to simulate the change of MAC but also need certain restrictions to prevent it from becoming a pure data fitting result. In the idealized numerical experiments, the IP number with the smallest RMSE is the optimal IP number. However, there is no preset MAC in the observed data, so the subjectivity of EHA is a weakness that needs to be overcome. For this reason, we try to propose two criteria to select the optimal IP number in Section4.4.


### 4.2. Extended to Time-Varying Phases

In fact, not only the amplitudes but also the phases of the annual cycles are modulated according to the observations. Thus, we also use EHA, EEMD, and NMD to extract the annual cycle with both amplitudes and phases modulated from constructed data. The new MAC can be expressed as the following form:(12)St=Atcos2πft+π180φt.The frequencyf (once a year) was fixed. The time-varying amplitude At and noise of the annual cycle were set as in Section 3, and the newly added phase change was designed as follows:(13)φt=−10+2Nt+10sin2πNt,t=1,2,…,N.The phaseφt varies between positive and negative 20°. Then, we obtained 7 more complicated MACs. The data used in the experiments 8–14 was constructed with the new MACs and noise. And the new MACs are extracted by the EEMD, NMD, and EHA. The RMSEs of experiments 8–14 are shown in Figure 4. It can be seen that, in experiments 8–10, the RMSEs of EHA are significantly larger than those without phase changes, but they are still the smallest among the three methods. And in experiments 11–14, EHA still shows the highest accuracy. As for the amplitude, the error of the NMD result becomes larger, and the EHA result is still closest to the preset amplitude (Figure 5), indicating the superiority of EHA.Figure 4
The RMSE of extraction results of annual cycle with variable amplitude and phase.Figure 5
The amplitude changes obtained by the three methods in experiment 11 with 1% white noise.


### 4.3. Impact of IP Number on the EHA Extraction Results

The extraction result of EHA depends on the IP number, which is a feature of the EHA method. According to Pan et al. [10], amplitudes and phases obtained by EHA using different IP numbers represent oscillations on different time scales. If the amplitudes and phases of the given MAC are very complicated, it is obvious that more IP numbers should be used in EHA. However, too many IP numbers will induce a seriously abnormal extraction result. In fact, for a given MAC, there always exists an optimal IP number that can describe it the best.For example, when the amplitude and phase are constant, the RMSE increased with IP number, and the RMSE was the smallest when IP number is one. If the amplitude and phase were not constant, the RMSE decreased first and then increased with the increasing IP number. In other words, in each experiment, an optimal IP number can always be found to minimize the deviation of the EHA result from the preset annual cycle (Tables1 and 2). The optimal IP number is determined by the complexity of the annual cycle and the noise power. In our idealized experiments, the optimal IP number increases with the complexity of MAC.Table 1
The IP number making the smallest RMSE in the annual cycle extraction of amplitude changes. NPR is the noise-power ratio.
NPR%EX 1EX 2EX 3EX 4EX 5EX 6EX 711248991421241010917412489814Table 2
The IP number making the smallest RMSE the smallest in the annual cycle extraction of simultaneous changes in amplitude and phase. NPR is the noise-power ratio.
NPR (%)EX 8EX 9EX 10EX 11EX 12EX 13EX 1416678991424461010917444489814


### 4.4. Annual Cycle Extraction of Observation Data

In order to test the performance of these methods (EEMD, NMD, and EHA) in the real observations, they are applied to the daily water temperature at the Rhône Sion Station from 1984 to 2013, for a total of 30 years. Sion Station is located at the bottom of a densely populated Alpine Valley through the Rhône River. We use EEMD, NMD, and EHA with 80 IPs (the reason for the selection of 80 is given later) to extract the annual cycle of this data (Figure6(a)). The lack of true annual cycle data made it impossible to judge the accuracy of their results. The trend of annual cycle amplitude changes obtained by the three methods was generally consistent, but the EEMD results changed more drastically, and the difference between NMD and EHA is relatively small (Figure 6(b)). It shows that EHA can also effectively extract the annual cycle from the real observed climate data and can decompose annual cycles into time-varying amplitudes and phases.Figure 6
(a) Annual cycles extracted from real observed data using three methods; (b) the same as (a), but from 16th year to 19th year.
(a)(b)Compared with the constructed data, the observation data of the temperature time series are more complicated. The way of selecting the optimal IP number has changed since we do not know the real MAC. As indicated before, the more the IP numbers, the more complicated the results obtained by EHA. As many IP numbers as possible should be used in the analysis of observed data. Too many IP numbers, however, will lead to serious abnormal extraction results. Then, the resulted time-varying amplitudes and phases will lose physical meaning and become only artifacts. Here, we attempt to propose two criteria for selecting optimal IP number when extracting MAC from observations using EHA. Firstly, if the IP number is appropriate, the time-varying amplitudeHt obtained by EHA usually fluctuates around the constant amplitude H′ obtained by the traditional fitting model [10]. Thus, the relative size of the time-averaged amplitude Ht¯ can be used to verify the reliability of the results. Ht¯ cannot exceed 5% of the amplitude H′ extracted by traditional fitting model, namely,(14)0.95H′≤Ht¯≤1.05H′.Secondly, S_TIDE will also provide a confidence interval of 95% for the obtained amplitude. We can compare the amplitude and interval length (i.e., the signal-to-noise ratio of amplitude) to determine if the result is reasonable:(15)SNRH=Ht¯Lt¯,where Ht¯ is time-averaged amplitude and Lt¯ is the corresponding time-averaged confidence interval. It is generally believed that SNRH greater than 2 is credible, and SNRH usually decreases with increasing IP number (Figure 7). In the observed series analyzed in this paper, the IP number less than 80 can meet the above criteria, so we treat 80 as the optimal IP number.Figure 7
Changes of SNR and mean amplitude with increasing IP number.In the observed data, there may be many cycles, such as the annual cycle and half-year annual cycle. The previous numerical experiments were based on a monofrequency MAC. Therefore, applying the EHA described in Section2.3 directly to the observed data may not perform as well as the experimental results. In fact, the description of EHA in Section 2.3 is a simplified version for the monofrequency MAC experiment. The complete EHA can process the data which has multiple cycles. When the data has multiple frequency cycles, the form of EHA changes from equation (5) to(16)St=S0+∑i=1naitcosfit+bitsinfit,where fi is the frequency of the i-th cycle and n is the number of cycles. After calculating equation (15), multiple modulated cycles with different frequencies can be obtained. For different types of data, which cycles should be used in EHA decomposition is also different, and its effect needs further verification.


## 4.1. Result

The previously constructed data has been processed by three methods, EEMD, NMD, and EHA, to extract different MACs. This section mainly discusses the deviation between the extraction results and the preset amplitude-modulated annual cycle. In order to avoid the edge error, only 46 years results were selected, excluding the first two years and the last two years. The deviation can be expressed by the root mean square error (RMSE) as follows:(11)RMSE=∑i=1NSi−Si′2N,where Si represents the extracted annual cycle and Si′ represents the preset annual cycle. N is the length of Si. The RMSE in each experiment is shown in Figure 2. The RMSE of the three methods increases with noise power. The RMSEs of EEMD are the largest and an order of magnitude higher than the RMSEs of other methods. When the changes of the preset annual cycle are relatively simple (the first three types), the RMSE of EHA is an order of magnitude smaller than that of NMD. When the preset annual cycle changes are complex, RMSEs of EHA and NMD are in the same order of magnitude. In each group of experiments, the RMSE of EHA is the smallest, EEMD is the largest, and NMD is between them. The RMSEs of NMD and EHA are also affected by the complexity of amplitude variation. At the same noise-power-ratio level, the more complex the amplitude variation is, the larger the RMSE will be. However, the RMSE of EEMD is almost only affected by the noise-power ratio and is not sensitive to the complexity of amplitude changes. EHA can directly obtain the instantaneous time-varying amplitude. For the extraction results of NMD and EEMD, the time-varying amplitudes of the MAC are obtained through a cubic spline fitting to the local maxima of MAC2 series [4]. The result of EHA is closest to the preset amplitude (Figure 3).Figure 2
The RMSE of extraction results with variable amplitude annual cycle.Figure 3
The amplitude changes obtained by the three methods in experiment 6 with 1% white noise.Among the three methods, both EEMD and NMD are adaptive, while EHA, which is subjective, needs users to tune the IP number. When the IP number is too small, the result obtained by EHA is relatively simple and can only reflect low-frequency changes of MAC. The more the IP is, the more the MAC changes can be reflected. When the IP number is too large, the result of EHA is likely to be severely abnormal. This means that we not only need a large enough IP number to simulate the change of MAC but also need certain restrictions to prevent it from becoming a pure data fitting result. In the idealized numerical experiments, the IP number with the smallest RMSE is the optimal IP number. However, there is no preset MAC in the observed data, so the subjectivity of EHA is a weakness that needs to be overcome. For this reason, we try to propose two criteria to select the optimal IP number in Section4.4.


## 4.2. Extended to Time-Varying Phases

In fact, not only the amplitudes but also the phases of the annual cycles are modulated according to the observations. Thus, we also use EHA, EEMD, and NMD to extract the annual cycle with both amplitudes and phases modulated from constructed data. The new MAC can be expressed as the following form:(12)St=Atcos2πft+π180φt.The frequencyf (once a year) was fixed. The time-varying amplitude At and noise of the annual cycle were set as in Section 3, and the newly added phase change was designed as follows:(13)φt=−10+2Nt+10sin2πNt,t=1,2,…,N.The phaseφt varies between positive and negative 20°. Then, we obtained 7 more complicated MACs. The data used in the experiments 8–14 was constructed with the new MACs and noise. And the new MACs are extracted by the EEMD, NMD, and EHA. The RMSEs of experiments 8–14 are shown in Figure 4. It can be seen that, in experiments 8–10, the RMSEs of EHA are significantly larger than those without phase changes, but they are still the smallest among the three methods. And in experiments 11–14, EHA still shows the highest accuracy. As for the amplitude, the error of the NMD result becomes larger, and the EHA result is still closest to the preset amplitude (Figure 5), indicating the superiority of EHA.Figure 4
The RMSE of extraction results of annual cycle with variable amplitude and phase.Figure 5
The amplitude changes obtained by the three methods in experiment 11 with 1% white noise.


## 4.3. Impact of IP Number on the EHA Extraction Results

The extraction result of EHA depends on the IP number, which is a feature of the EHA method. According to Pan et al. [10], amplitudes and phases obtained by EHA using different IP numbers represent oscillations on different time scales. If the amplitudes and phases of the given MAC are very complicated, it is obvious that more IP numbers should be used in EHA. However, too many IP numbers will induce a seriously abnormal extraction result. In fact, for a given MAC, there always exists an optimal IP number that can describe it the best.For example, when the amplitude and phase are constant, the RMSE increased with IP number, and the RMSE was the smallest when IP number is one. If the amplitude and phase were not constant, the RMSE decreased first and then increased with the increasing IP number. In other words, in each experiment, an optimal IP number can always be found to minimize the deviation of the EHA result from the preset annual cycle (Tables1 and 2). The optimal IP number is determined by the complexity of the annual cycle and the noise power. In our idealized experiments, the optimal IP number increases with the complexity of MAC.Table 1
The IP number making the smallest RMSE in the annual cycle extraction of amplitude changes. NPR is the noise-power ratio.
NPR%EX 1EX 2EX 3EX 4EX 5EX 6EX 711248991421241010917412489814Table 2
The IP number making the smallest RMSE the smallest in the annual cycle extraction of simultaneous changes in amplitude and phase. NPR is the noise-power ratio.
NPR (%)EX 8EX 9EX 10EX 11EX 12EX 13EX 1416678991424461010917444489814


## 4.4. Annual Cycle Extraction of Observation Data

In order to test the performance of these methods (EEMD, NMD, and EHA) in the real observations, they are applied to the daily water temperature at the Rhône Sion Station from 1984 to 2013, for a total of 30 years. Sion Station is located at the bottom of a densely populated Alpine Valley through the Rhône River. We use EEMD, NMD, and EHA with 80 IPs (the reason for the selection of 80 is given later) to extract the annual cycle of this data (Figure6(a)). The lack of true annual cycle data made it impossible to judge the accuracy of their results. The trend of annual cycle amplitude changes obtained by the three methods was generally consistent, but the EEMD results changed more drastically, and the difference between NMD and EHA is relatively small (Figure 6(b)). It shows that EHA can also effectively extract the annual cycle from the real observed climate data and can decompose annual cycles into time-varying amplitudes and phases.Figure 6
(a) Annual cycles extracted from real observed data using three methods; (b) the same as (a), but from 16th year to 19th year.
(a)(b)Compared with the constructed data, the observation data of the temperature time series are more complicated. The way of selecting the optimal IP number has changed since we do not know the real MAC. As indicated before, the more the IP numbers, the more complicated the results obtained by EHA. As many IP numbers as possible should be used in the analysis of observed data. Too many IP numbers, however, will lead to serious abnormal extraction results. Then, the resulted time-varying amplitudes and phases will lose physical meaning and become only artifacts. Here, we attempt to propose two criteria for selecting optimal IP number when extracting MAC from observations using EHA. Firstly, if the IP number is appropriate, the time-varying amplitudeHt obtained by EHA usually fluctuates around the constant amplitude H′ obtained by the traditional fitting model [10]. Thus, the relative size of the time-averaged amplitude Ht¯ can be used to verify the reliability of the results. Ht¯ cannot exceed 5% of the amplitude H′ extracted by traditional fitting model, namely,(14)0.95H′≤Ht¯≤1.05H′.Secondly, S_TIDE will also provide a confidence interval of 95% for the obtained amplitude. We can compare the amplitude and interval length (i.e., the signal-to-noise ratio of amplitude) to determine if the result is reasonable:(15)SNRH=Ht¯Lt¯,where Ht¯ is time-averaged amplitude and Lt¯ is the corresponding time-averaged confidence interval. It is generally believed that SNRH greater than 2 is credible, and SNRH usually decreases with increasing IP number (Figure 7). In the observed series analyzed in this paper, the IP number less than 80 can meet the above criteria, so we treat 80 as the optimal IP number.Figure 7
Changes of SNR and mean amplitude with increasing IP number.In the observed data, there may be many cycles, such as the annual cycle and half-year annual cycle. The previous numerical experiments were based on a monofrequency MAC. Therefore, applying the EHA described in Section2.3 directly to the observed data may not perform as well as the experimental results. In fact, the description of EHA in Section 2.3 is a simplified version for the monofrequency MAC experiment. The complete EHA can process the data which has multiple cycles. When the data has multiple frequency cycles, the form of EHA changes from equation (5) to(16)St=S0+∑i=1naitcosfit+bitsinfit,where fi is the frequency of the i-th cycle and n is the number of cycles. After calculating equation (15), multiple modulated cycles with different frequencies can be obtained. For different types of data, which cycles should be used in EHA decomposition is also different, and its effect needs further verification.


## 5. Conclusion

The extraction of annual cycles has great significance in the study of climatic and oceanic data. We design a set of idealized experiments to test the effect of EHA on extracting the annual cycle. In the idealized experiments involving the annual cycle with only amplitude changed, the result of EHA has the smallest deviation from the preset annual cycle. The deviations of EHA and NMD extraction results are affected by the complexity of the annual cycle and the noise power. In the experiments extended to time-varying phases, the RMSEs of EHA with modulated phases are significantly larger than the RMSEs with a fixed phase, but EHA can still generate the most accurate result among the three methods.The performance of EHA is affected by the IP number. For each group of experiments, the optimal IP number minimizing the deviation of the extraction results from the preset annual cycle can always be found. For the actual observation data, due to the lack of the reference annual cycle, the definition of the optimal IP number has changed. Thus, we proposed two criteria for selecting the optimal IP number. In future work, the new method explored in this study is expected to process a variety of climatic and oceanic data.



---

*Source: 9625795-2021-09-30.xml*