Applicability Study of Hydrological Period Identiﬁcation Methods: Application to Huayuankou and Lijin in the Yellow River Basin, China

: Identifying implicit periodicities in hydrological data is signiﬁcant for managing river– basin water resources and establishing ﬂood forecasting systems. However, the complexity and randomness of hydrological systems make it difﬁcult to detect hidden oscillatory characteristics. This study discusses the performance and applicability of ﬁve period identiﬁcation methods, namely periodograms, autocorrelation analysis (AA), maximum entropy spectral analysis (MESA), wavelet analysis (WA), and the Hilbert–Huang transform (HHT). The annual and monthly runoff data are sampled from two stations (Huayuankou and Lijin on the Yellow River in China) in the years 1949–2015. The conclusions are as follows: (i) All methods identify the signiﬁcant periods of 6 months, 12 months, and 18–19 months, which have relatively high energy of peaks; (ii) WA and HHT perform best when dealing with nonstationary time series, but they are ineffective for identifying large-scale periods; (iii) MESA has high resolution and stability but is prone to oscillate at small-scale periods when applied to monthly series; and (iv) periodograms and AA are relatively simple, but their results lack stability and are signiﬁcantly affected by the data length—the resolution of AA is too low when applied to annual data, and periodograms can easily produce “false peaks”. Generally, it is better to apply multiple methods comprehensively than each method singularly, and this can be effective in reducing subjective inﬂuences.


Introduction
The hydrological cycle is a particularly complex global system driven by solar radiation and gravity [1,2]. Hydrological time series (HTS) contain the influences of numerous processes involved in the transfer of water in the hydrological cycle [3,4], and they are influenced by many physical factors that are often interrelated, especially large-scale fluctuations in atmospheric circulation, the Earth's rotation and revolution, and sunspots [5][6][7][8]. Therefore, to reveal the variability of hydrological processes, it is essential to be able to analyze HTS that have complicated stochastic characteristics [9,10]. The identification of hydrological periods is a key issue in the analysis of HTS and is of great significance for hydrological simulation and prediction, hydrological design, and the planning and management of water resources [11][12][13]. However, because strict periodicities are difficult to obtain, we are restricted to approximate implied periods, also known as quasi-periodicities or potential cycles [14].
The methods currently used to identify periodicities in HTS have different theoretical bases and levels of performance. The main methods are three types of spectral analysis, namely Fourier methods (FMs), wavelet analysis (WA), and the Hilbert-Huang transform (HHT) [15,16]. The FMs referred to herein are forms of spectral analysis based on the principle of Fourier analysis, including periodograms, autocorrelation analysis (AA), and maximum entropy spectral analysis (MESA).
Traditional spectral analysis methods, mainly periodograms and AA, are widely used and have high computational efficiency when a fast Fourier transform is applied (periodogram) or a few lags are needed (AA) [17]. However, they also have some inevitable shortcomings. For instance, the frequency resolution of a periodogram is approximately equal to the reciprocal of the observation data length, so it is limited by the duration of the available record. The rough spectrum of AA is smoothed by a window, meaning that the statistical stability and resolution depend on the subjective choices of the maximum autocorrelation lag and the type of window function [18,19]. Since Fourier algorithms require integration over infinite time but the actual observed data are limited, periodograms and AA make certain artificial assumptions about data extension. Due to the inherent window effect, the spectrum will suffer "leakage" that distorts it and may cause a low amplitude sinusoid to be submerged completely by the side lobes of a sinusoid of greater amplitude [20]. To improve the performance of traditional spectral analysis, MESA based on the principle of maximum entropy (POME) has been proposed. This assumes that unknown information has maximum uncertainty and randomness, so the entropy of the probability density function (PDF) should be kept at a maximum value [21]. For onedimensional analysis of wide sense stationary, Gaussian processes, using a filter of order m to estimate the maximum entropy spectrum is identical to fitting an autoregressive (AR) model of order m to the data and estimating its spectrum by the fitted parameters [17,18]. MESA methods have relatively high resolution for low noise levels and good spectral fidelity for short data records. However, traditional approaches perform better for processes whose signal-to-noise ratio (SNR) is relatively low. Spectral line splitting and biased frequency estimates may occur when using the Burg algorithm for MESA [22], so the results of period identification must be discriminated with care whenever either an HTS contains an appreciable level of noise or spectra with multiple peaks are encountered [23,24].
While FMs are powerful tools for detecting average periodicities in hydrological data, they cannot show the distribution of periodic variabilities with time [25]. Moreover, Fourier analysis is only suitable for stationary and ergodic time series [26,27]. Wavelet analysis (WA) is derived from Fourier analysis, however, it has further development. Wavelet functions can show the localized and transient phenomena occurring in the time domain through the window with a fixed area but a changeable shape [28,29]. WA is suitable for studying HTS with multitemporal scales and nonstationary characteristics [30,31]. However, it should be noted that the length of the measured runoff sequence is limited by the observation conditions. Boundary effects may occur when applying WA to period identification, especially in relatively short data.
While Fourier analysis uses trigonometric functions as its basis, which may obscure the real characteristics of the time series, the selection of prototype functions also limits the application of WA [32,33]. To alleviate this drawback, Huang et al. [34] developed a new time-series analysis technique, namely the HHT, for nonlinear and nonstationary hydrological data. The HHT is based on empirical mode decomposition (EMD), which is an adaptive decomposition method based on local temporal scales. By assuming that the signal consists of different intrinsic modes of oscillations, EMD separates it into a series of modes known as intrinsic mode functions (IMFs), each of which is associated with a welldefined timescale. To identify a set of IMFs without further decomposition of the residual, this procedure emphasizes completeness, orthogonality, locality, and adaptiveness in the time domain [35]. Applying a Hilbert transform to each IMF leads to the Hilbert spectrum, which reflects the relationship among time, amplitude, and frequency. Compared to the harmonic functions of FMs, whose amplitude and frequency are constant a priori, the IMFs of HHT represent an a posteriori-defined basis, which is nonlinear and nonstationary in amplitude and frequency [36]. EMD is widely used because of its capacity to characterize the multitemporal scales directly. However, signal intermittency and noise may cause mode mixing that obscures the physical meaning of individual IMFs and renders the EMD algorithm unstable [37]. To improve the effectiveness of EMD, ensemble EMD (EEMD) has been proposed, which involves adding white noise of a particular amplitude to the original data [38]. Additionally, some other issues should also be considered when decomposing the original sequence, such as the selection of the fitting curve for extreme points, the treatment of boundary points, and the criteria for ending the sifting process of IMF [14].
In addition to applying the aforementioned single methods, many scholars have explored multimethod ways to identify the multitemporal scales in HTS. For instance, Sang et al. [27] studied the relationship between period identification and noise and proposed a new method, namely main series spectral analysis (MSSA), which is based on two conventional methods of HTS analysis, namely WA and MESA, to improve period identification by reducing noise interference. Sang et al. [14] used empirical mode decomposition (EMD) to decompose an HTS into a set of intrinsic mode functions (IMFs), whereupon the noise could be distinguished so that the true IMFs could be analyzed by MESA. This EMD-MESA method can improve overall period identification by distinguishing among noise, periods, and trends. Yu et al. [39] combined continuous wavelet transform (CWT) and HHT, effectively eliminating the critical drawback in CWT regarding fine-scale mode mixing and revealing multiscale periodicities in precipitation signals.
To summarize, although there are various approaches to identifying periods in HTS, previous studies have lacked systemic analysis of: (i) The conditions under which such methods are applicable and; (ii) How the consequences of such methods vary.
Considering certain limitations, the method to be applied should be chosen according to the application environment. As such, the present study analyzes the observed annual and monthly runoff data (series on two different time scales) of the Huayuankou hydrologic station (HHS) and Lijin hydrologic station (LHS) along the Yellow River in China to explore the applicability of the aforementioned methods in period identification and to improve their efficiency and accuracy.

Study Area and Data
The Yellow River basin lies between 95 • 53 -119 • 05 E and 32 • 10 -41 • 50 N. For this study, we selected HHS and LHS on the Yellow River as study sites, as shown in Figure 1. HHS is the demarcation point between the middle and lower reaches of the Yellow River basin, and plays a key role in preventing and mitigating downstream flooding. The control area covers 97% of the total catchment area, which is about 730,000 km 2 , and the water flow accounts for 96.6% of the total water of the Yellow River basin. LHS is the final hydrological control station on the Yellow River before it reaches the sea. The exchange of material and energy between land and sea is most frequent at the estuary and has the most impact on runoff. LHS is 103.6 km from the sea entrance and controls an area of roughly 751,900 km 2 , accounting for 99.92% of the total basin area. Therefore, it is very important to study the characteristics of the spatial and temporal variation of runoff at these two stations. In this study, we used the annual and monthly mean runoff series, which were calculated from the arithmetic average of the observed daily runoff data. Sixteen sets of runoff data adopted to period identification are listed in Table 1. The specific records in this study were provided by the Yellow River Conservancy Commission (YRCC) and the Yellow River Water Resources Bulletin (YRWRB).
13, x FOR PEER REVIEW

The Mann-Kendall Test
The underlying trend existed in HTS can affect the stability of the data and sho eliminated before the period analysis [40]. Methods for analyzing trends in HTS a ally divided into parametric tests and non-parametric tests. Non-parametric tests, ing Spearman's rho test and the Mann-Kendall (MK) test, are often preferred becau are robust with respect to non-normality, nonlinearity, missing values, serial depen censored data, and outliers (extremes) [41]. Herein, we used the MK test to ident nificant trends in HTS, this being a rank-based distribution-free method; see Ham for specific details. In this study, all series show a decreasing trend at the significanc of 0.01. We used the least-squares method to fit and then remove the trend in e quence before analyzing periods. It should be noted that with the HHT, the origin  Note: the unit of series I is 'year', which is represented by 'a'; the unit of series II is 'month', which is represented by 'm'.

The Mann-Kendall Test
The underlying trend existed in HTS can affect the stability of the data and should be eliminated before the period analysis [40]. Methods for analyzing trends in HTS are usually divided into parametric tests and non-parametric tests. Non-parametric tests, including Spearman's rho test and the Mann-Kendall (MK) test, are often preferred because they are robust with respect to non-normality, nonlinearity, missing values, serial dependency, censored data, and outliers (extremes) [41]. Herein, we used the MK test to identify significant trends in HTS, this being a rank-based distribution-free method; see Hamed [42] for specific details. In this study, all series show a decreasing trend at the significance level of 0.01. We used the least-squares method to fit and then remove the trend in each sequence before analyzing periods. It should be noted that with the HHT, the original data will be used directly for period identification.

Periodogram
Schuster [43] proposed the concept of the periodogram, which is applicable to any time sequence that has natural periodicity or that consists of harmonic signals imbedded in noise [20,44]. Periodogram directly estimates the power spectrum of the original sequence. It is widely used and has high computational efficiency when a fast Fourier transform is applied [17]. The specific calculation steps can be found in Ding and Deng [15].

Autocorrelation Analysis (AA)
Compared to the periodogram, AA indirectly obtains the power spectral density (PSD) by Fourier transform of the autocorrelation function [45]. The spectral density function, S f j , is calculated by Equation (1): where k is the step and m is the maximum lag of autocorrelation; r k is the autocorrelation coefficient of the sample; f j = j 2m with j varying from 1 to m; and D k is the window. This study adopted hanning window as the form of D k , as shown in Equation (2): More information about AA can also be found in Ding and Deng [15].

Maximum Entropy Spectral Analysis (MESA)
Limited by the duration of available records, periodograms and AA both make artificial assumptions about data extension. To overcome this shortcoming, MESA based on the principle of maximum entropy (POME) has been proposed [22]. For one-dimensional analysis of wide sense stationary, Gaussian processes, using a filter of order m to estimate the maximum entropy spectrum is identical to fitting an autoregressive (AR) model of order m to the data and estimating its spectrum by the fitted parameters [17,18]. The AIC criterion and trial-and-error testing were selected to determine the optimal model order. Additionally, m is the optimal order when the AIC(m) reaches minimum value, which is shown in Equation (3): In this study, both HHS and LHS take the AR model order of 9 in the annual series, while for the monthly series, HHS takes the order of 38 and LHS takes the order of 49. More details about MESA can be found in Padmanabhan and Ramachandra Rao [18].

Wavelet Analysis
FMs are powerful tools for detecting average periodicities in hydrological data, but they cannot show the distribution of periodic variabilities with time [25]. Moreover, Fourier analysis is only suitable for stationary and ergodic time series [26,27]. In the early 1980s, Morlet proposed WA based on the short Fourier transform. WA is suitable for studying HTS with multitemporal scales and nonstationary characteristics [30,31]. Wavelet functions are more diverse than the basis functions of Fourier analysis [46]. They can show the localized and transient phenomena occurring in time domain through the window with a fixed area but a changeable shape [28,29]. In this study, Morlet wavelet consisting of a plane wave modulated by a Gaussian is adopted as expressed in Equation (4): where η is dimensionless time and ω 0 is dimensionless frequency. To satisfy the admissibility condition, the value of ω 0 is usually taken to be 6. More information about WA can be seen in Torrence and Compo [47].

The Hilbert-Huang Transform
FMs use trigonometric functions as basis, which may obscure the real characteristics of the time series, and the selection of prototype functions also limits the application of WA [32,33]. To alleviate this drawback, Huang et al. [34] developed the HHT. It consists of two important parts, namely EMD and the Hilbert transform. According to the internal features of an HTS, EMD decomposes the signal adaptively and efficiently into several IMFs with different characteristic time scales [14,48]. The IMF represents an a posteriori-defined basis, which is nonlinear and nonstationary in amplitude and frequency [36]. When the decomposition process ends, the original sequence can be expressed as Equation (5): where c i (t) with i varying from 1 to p represents a series of IMFs and r P (t) is the residual. Applying a Hilbert transform to each IMF leads to the Hilbert spectrum, which reflects the relationship among time, amplitude, and frequency. The specific process of HHT and more details are described as Huang et al. [34].

Period Identification Results
In this study, we calculated the aforementioned 16 sets of data and analyzed them by periodogram, AA, MESA, WA, and HHT. During spectral analysis, there were no stable peaks if the period exceeded a certain range. Therefore, only partial spectra were given below to ensure the integrity and clarity of the recognition results. In addition, if the period identified in a monthly series exceeded 24 months, the results were not sufficiently stable at different series lengths and it was difficult to find regularities directly when viewed on a monthly scale. Therefore, turning these periods into an annual scale facilitated not only analysis but also comparison with the identification results for the annual series (minimum two years). Therefore, we divided the periods identified from the monthly series into two categories: (i) Periods of less than 24 months (which are recorded directly) and; (ii) Periods of 24 months or more.
We measured the second category of period in years to analyze whether the periods in annual and monthly mean runoff series were consistent.
The periodogram spectra are shown in Figures 2 and 3. When adopting the AA method, it is necessary to determine the maximum delay m in different series lengths. Since all sample lengths n are greater than 50, we require m < n/4 and usually take m ≈ n/10 [15]. Herein, we took three values of m for each data set: m = 8, 10, and 12 for annual series and m = 80, 100, and 120 for monthly series. Figures 4 and 5 show the power spectra for the different values of m. We approximated the refined discrete frequency to the continuous frequency, whereupon we obtained the continuous MESA spectrum. However, if the frequency was insufficiently accurate, the spectrum was insufficiently smooth. Herein, the initial frequency was 0.001 (1/T, T is period) and the frequency interval (accuracy) was 0.001 (1/T). When using MESA, both HHS and LHS took the AR model order of 9 in the annual series, while for the monthly series, HHS took the order of 38 and LHS took the order of 49. Finally, the maximum entropy spectra are shown in Figures 6 and 7. In the CWT, the wavelet scale (a) cannot represent the period directly, which need be converted into an appropriate scale sequence according to the arithmetic sequence of the corresponding frequency [26]. For WA, scale (a) should range from 2f to infinity, but herein the scale maximum was 1500 (the unit of scale (a) was equal to period 'T'), which corresponded to a frequency interval of 0.001 (1/T); realistically, this scale was large enough for practical applications. We then used "symmetry extension" to reduce the boundary effect, that is, we extended the length of both ends of the series symmetrically once. Figures 8 and 9 show the wavelet variance diagrams. In the application of HHT analysis, EMD was used to decompose a runoff sequence into a residual and various IMFs. The results regarding the decomposition of sequences 67a, 67m, 66a, and 66m are shown in Figures 10-13, respectively. The residual of each sequence shows a decreasing trend, which was in keeping with the results of previous MK tests. A Hilbert transform was then applied to each IMF to obtain the Hilbert spectra, as shown in Figures 14-17, respectively. The frequency corresponding to each IMF was not constant but fluctuated around a central frequency; the higher the frequency, the greater the fluctuation energy (amplitude squared). Additionally, the fluctuation range was limited and there were relatively few overlaps between different components, so each IMF component had a relatively clear distribution. The central frequencies and the corresponding average periods of the 16 sequences are given in Tables 2-5. Finally, during the five identification processes, we excluded some periods that were not stable at different sequence lengths. In addition, we regarded some low-energy peaks near the main period as oscillations and did not list them separately. The specific period identification results are given in Tables 6-9.                                          10-11 6-7 9-10 9-10 4-5 4-5 20-21 6-7 3

Intrinsic mode functions
Note: the scale of the period is 'year'; the minor periods of low energy and some oscillation components are not listed. Note: 'I' is the first category of the periods, of which the scale is 'month'; 'II' is the second category of the periods, of which the scale is 'year'; in addition to the main periods listed in this table, some periods of low energy and some oscillation components are not listed. Table 8. The periods identified from annual series at Lijin station.  From Table 6, the period recognition results for series 67a, 62a, 57a, and 52a were 3 years, 4-5 years (periodogram, MESA, and WA), 6-7 years (periodogram and WA), 9-11 years (all methods except AA), and 20-21 years (HHT). The second-category periods for series 67m, 62m, 57m, and 52m series are given in Table 7 as 2-3 years (all methods), 3-5 years (periodogram, WA, and HHT), 6-7 years (periodogram and HHT), 8-10 years (periodogram and WA), and 21 years (WA). Comparing these two sets of results shows good correspondence, that is, the periods identified in the annual series can also be found in the monthly series. The first-category periods for series 67m, 62m, 57m, and 52m are also given in Table 7. All methods identify the significant periods of 6 months, 12 months, and 18-19 months, which have relatively high energy of peaks, especially the period of 12 months. Additionally, periodograms and the HHT identify periods of 9-10 months and 7-8 months, respectively. Additionally, the three FMs all identify periods of 4 months.

Periodogram
From Table 8, the period recognition results of series 66a, 61a, 56a, and 51a were 3 years (all methods except HHT), 4-5 years and 9-11 years (all methods except AA), 6-7 years (periodogram and WA), and 20-22 years (HHT). The second-category periods for series 66m, 61m, 56m, and 51m are given in Table 9 as 3 years (all methods), 4-5 years and 6-7 years (periodogram, WA, and HHT), 9-12 years (all methods except AA), and 20-21 years (WA). Comparing these two sets of results also shows good correspondence. The first-category periods for series 66m, 61m, 56m, and 51m are also given in Table 9. All methods identify the significant periods of 6 months, 12 months, and 17-19 months, AA and MESA both identify the period of 4 months, and HHT and periodograms identify the periods of 7-8 months and 9-11 months, respectively.
These five methods identify different periods, some of which (the obvious ones) were identified by most of the methods while others were identified by only some of the methods. In general, overall comparison shows that (i) the results of periods identification obtained from HHS were relatively consistent to those obtained from LHS and (ii) the second-category periods in the monthly series corresponded reasonably well with those in the annual series.
The specific period identification results are given in Tables 6-9. In general, overall comparison shows that: (i) The results of periods identification obtained from HHS are relatively consistent to those obtained from LHS and (ii) The second-category periods (24 months or more) in the monthly series correspond reasonably well with those in the annual series.

Period Validation
This study shows the implicit runoff periodicity, mainly the periods of 3-5 years, 6-7 years, 9-11 years, and 20-21 years at HHS and 2-5 years, 6-7 years, 9-11 years, and 20-22 years at LHS. Sang et al. [49] conducted cross-correlation WA of two runoff series from HHS and LHS and found them to be well cross-correlated on four time scales, namely 3 years, 7 years, 11 years, and 20 years, which correspond to the four evident periods of each runoff series. Wang and Zhu [50] used the POME-based MEM1 spectrum to analyze the monthly HHS runoff series from 1952 to 1990 and obtained the main periods of 12 months, 6 months, 4 months, and 3 months. Those periods were basically consistent with the results in the present study.
Potential physical causes about possible periods were also analyzed. Deterministic periods of 24 h and 12 months (1 year) were due mostly to the Earth's rotation and revolution, respectively. Runoff periods of less than 12 months may be related to seasonal changes, and other periods may be due to the air-sea interaction or sunspot activity [51]. For example, the Western Pacific subtropical high has the 3-4 years quasi-period [52]; the polar-motion amplitude variation has the 6.5 years period [53]; the sunspot-cycle period is about 11 years [54] and solar activity has the 22 years quasi-period [51]. In fact, the complex connections among the air-sea interaction and solar activity can also impact runoff either directly or indirectly [55,56]. Most are reflected intuitively in variations of climate factors such as precipitation, evaporation, and temperature, which have significant influences on runoff [57]. In addition, runoff variations can also be affected by topography, human activities, and other factors. The specific mechanisms of these influences require further exploration.

Performance of Five Methods
The periodogram could identify most periods, however, an obvious "false peak" problem makes the identification process difficult. Meanwhile, its results lacked stability and consistency at different series lengths, especially in the annual series. MESA performed well with both annual and monthly series. AA performed well with monthly series, but its resolution was particularly low when applied to annual series. Additionally, both MESA and AA produced obvious short-period oscillations. WA had good stability at short periods, but the spectral lines began to shift as the period increased and this phenomenon was more obvious with shorter series. Compared to FMs, WA did not oscillate at short periods and were affected less by the series length. The HHT obtained the central frequencies by modulating the IMFs directly and then transformed them into periods, meaning that the identification results were no longer influenced by the spectral shape. The periods were stable and consistent when the IMF frequency was relatively high (i.e., the period was short), whereas the results became unstable as the IMF frequency decreased (i.e., the period increased).
With the increasing demand for industrial water, agricultural irrigation, and domestic water in the Yellow River Basin, the contradiction between water supply and demand in the Yellow River Basin has become increasingly severe. Water shortage has become an important factor restricting regional socioeconomic development and agricultural production. Studying the characteristics of the Yellow River runoff change was directly related to the efficient use and scientific allocation of water resources in the basin and areas along the Yellow River. Related research should pay attention to the periodicity and range of runoff changes. The Yellow River runoff is affected by human activities, climate change, and other factors, which may cause interference in the period identification process. To determine significance levels, it is first necessary to choose an appropriate background spectrum. For many geophysical phenomena, an appropriate background spectrum is either white noise (with a flat Fourier spectrum) or red noise (increasing power with decreasing frequency) [47]. A potential signal depends quite sensitively on our a priori assumptions regarding the nature of the background noise, and the appropriateness of our statistical model for the signal itself [58,59]. It is worth noting that the study recognizes implicit period components and provides a reference standard to choose period identification methods. Therefore, if further analysis or prediction is required, that is, when each period component needs to be analyzed more accurately, different methods should be selected for noise elimination and period significance test according to the applicable conditions.

Theoretical Analysis
The frequency resolution of MESA, WA, and the HHT was relatively high and was not easily affected by the data length, while the recognition abilities of periodograms and AA depended greatly on the data length. A periodogram spectrum is discrete: as the data length n decreases, the frequency interval 1/n increases and the priori frequency points decrease. Thus, its spectral resolution is lower with fewer data points, furthermore, the recognized periods basing different sequence lengths become unstable. The frequency interval of AA is 1/2m (where m is the maximum lag and its value is related to n). Herein, m values of different series were equal, so the recognized periods based on different data lengths could be stable. However, AA assumes that the autocorrelation coefficient approaches zero rapidly [18,20], thereby requiring a sufficiently long data record to achieve the required frequency resolution. For different sequence lengths, the frequency intervals of MESA and WA were identical; herein, both were 0.001 (1/T, T is period). Besides, sufficiently refined discrete frequencies were selected to approximate the continuous spectra in MESA. CWT was adopted in WA and proper scales could be selected to meet the required frequency precision. The HHT also had good resolution when the decomposition levels were adequate.
In spectral analysis, waves of multiple frequencies were used to fit the original fluctuation, thereby producing some extra high-or low-frequency components that may not actually exist. The IMFs derived from EMD have certain physical meanings, allowing such "false peak" problems to be avoided effectively. Theoretically, Fourier spectral analysis is the most effective approach if the signal is stable and its spectral characteristics are clearly distinguishable from noise. For nonstationary series, however, the distribution characteristics of random variables cannot be inferred by statistical indicators, so problems such as "false regression" can occur [26,60]. WA and HHT are more suitable for studying nonstationary data with multitemporal scales. WA can reveal the local characteristics of time series from the time and frequency domains [30,31], while HHT can directly give the Hilbert spectrum, which can reflect the instantaneous changes of frequency and amplitude [34].

Artificial Influences
Period identification results can be affected by certain artificial choices when applying these five methods alone. This is reflected mainly in the selection of relevant parameters and various application conditions. Under certain circumstances, the relevant parameters as selected by the researcher may produce differing degrees of deviation, but there are usually no exact standards or criteria for such choices. This also makes it difficult to control artificial influences when applying a single method. The following are some related options that may impact the results.
The subjective choices of the maximum time lag m and the window function both impact the AA calculation results. The MESA results depend greatly on the signal-to-noise ratio and the initial phase. If the initial parameters of the entropy spectra are selected improperly, spectral peak shifting and line splitting may appear [22,23]. The optimal order of the AR model will affect the identification results in practice. AIC can be affected significantly by the sample size when used alone [61]. For example, the model order is clearly small when applying AIC to the annual series in this study. WA is affected mainly by the selection of the wavelet function and the admissibility conditions. Shorter series are affected more easily by the "boundary effect" and the original series should be prolonged to eliminate this influence before using WA. Since WA requires the wavelet radix and the decomposition layer to be chosen in advance, it is not adaptive. By contrast, the HHT is adaptive, having no fixed a priori basis, and the IMFs are extracted based on the temporal characteristics of the data. However, there are incompletely solved problems in HHT, such as the "end effect", the "under-envelope" phenomenon caused by cubic spline fitting, and the "negative frequency" phenomenon that occurs during the demodulation process [14,37].

Conclusions
The main purpose of this study was to assess the period identification abilities of five methods, namely periodograms, AA, MESA, WA, and the HHT. Periodograms and AA are simple and easy to operate but their identification abilities are susceptible to the data length; in particular, the resolution is too low when AA is applied to short series. For periodograms, the spectral lines are prone to shift and the results contain many false peaks. MESA performed relatively well in period identification but oscillated easily at small periods, especially in long series, and this phenomenon also existed in AA. Compared with FMs, WA can give relatively ideal spectra. The HHT can avoid the spectrum shape influencing the recognition results. Theoretically, WA and the HHT are more suitable for nonstationary time series, but those two methods are also not very effective in identifying large periods.
If a single algorithm is used directly, the related parameters must be adjusted according to the results, making the method subject to artificial influences. Furthermore, each method depends greatly on the relevant applicability conditions. These factors may cause bias in the identification process. This study has shown that the results of these five methods under consideration are strongly connected, that is, they are not completely independent. Therefore, selecting appropriate methods based on different applicability conditions and then combining the related results can give more complete and reliable periods. In brief, such a comprehensive method is less susceptible to artificial factors, its operation is intuitive and feasible, and the results are more accurate.