Abstract

According to the correlation of tree ring and solar activity, the cycle analysis method based on variational mode decomposition (VMD) and Hilbert transform is proposed. Firstly, the tree ring width of cypress during 1700 to 1955 beside the Huangdi Tomb and the long-term sunspot number during 1700 to 1955, respectively, are decomposed by VMD into a series of intrinsic mode functions (IMFs). Secondly, Hilbert transformation on the decomposed IMF component is performed. Then, the marginal spectra are given and analyzed. Finally, their quasiperiodic properties are obtained as follows: the tree ring width has the quasiperiodicity of 2 to 7a, 10.8a, and 25a; the sunspot number has the quasiperiodicity of 8.3a, 9.9a, 11.1a, 52.2a, and 101.2a. The result obtained by analyzing that quasiperiodicity shows that the main periods of tree ring width and the sunspot number in the same period are basically consistent, and tree ring width has other cycles. This shows that sunspot activity is an important factor affecting tree ring growth, and tree ring width is influenced by other external environments.

1. Introduction

The influence of solar activity on climate has attracted much attention as early as the seventeenth century. Solar activity is an important controlling factor of the earth’s climate. It plays a leading role in the earth’s climate change [1]. Many research studies [24] have shown that solar activity is closely related to natural environment changes such as temperature and various dust storms, floods, and droughts. Therefore, it affects people’s production and life. If the prevention is not proper, it may bring great loss to the normal life of human beings. Sunspot [5] is an important parameter of solar activity. Tree ring [6, 7] has some characteristics such as accurate age determination, strong continuity, high resolution rate, and being easy to acquire the duplicate. It records the processes of climate and environmental change and is an important global climate change research data [8]. Therefore, the correlation between tree ring and sunspot activity has important significance for people’s life.

Many scholars [913] have researched the correlation between solar activity and tree ring by using a lot of data. To research the correlation between solar activity and tree ring growth, the main way is to calculate whether there is a common change cycle. In 1979, Jr et al. [9] validated the phase locking degree between the change in the arid western United States and the Hale sunspot cycle since 1700 and suggested that the drought rhythm is in some manner controlled by long-term solar variability directly or indirectly related to the solar magnetic effect. In 1998, Zhou et al. [10] measured the correlation coefficient between the ring index and the sunspot cycle length with tree ring data over 598 years and confirmed that there was a negative maximum correlation coefficient between sunspot cycle length and climate. In 2002, Rigozo et al. [11] analyzed the sunspot number and tree ring width time series of Concadia in Brazil from 1837 to 1996 and verified the influence of solar activity on the growth of tree ring. In 2004, Miyahara et al. [12] studied the content of radioactive carbon in the tree ring from 1645 to 1715. The results showed that sunspot could maintain periodic polarity reversal during the minimum extension period. In 2017, Yin et al. [13] used EMD to analyze the stalagmite and tree ring time series of Holocene and discussed the possible effects of solar activity on climate change at the millennium scale.

Climate signal is a nonlinear and nonstationary signal. Since both sunspot number and tree ring width belong to the climate signal [14], it is not possible to directly analyze the periodic characteristics by using the traditional power spectrum estimation method or the wavelet analysis method based on the prior basis function. The application of convolution in the calculation process of wavelet will lead to mutation at the end point, and the wavelet does not have self-adaptability [15]. If power spectrum analysis and wavelet analysis are used to calculate the periodic characteristics of the sunspot number and tree ring, the calculation results are not ideal [16]. Empirical mode decomposition (EMD) [17, 18] provides a new idea for processing the nonlinear and nonstationary signal. When EMD is applied to other fields [19], it is found that mode decomposition may not be unique in physics, and sometimes there are mode mixing and end-point flying wings. The introduction of ensemble empirical mode decomposition (EEMD) [20] and variational mode decomposition (VMD) [21] has improved the mode mixing to some extent. Especially, VMD transforms signal decomposition into the nonrecursive variational mode. Its number of components is also less than that of EMD and EEMD, and it shows better noise robustness [22]. The resolution of the Hilbert spectrum in the time-frequency domain is much higher than that of the wavelet spectrum, so the result can accurately reflect the original physical characteristics of the system. The marginal spectrum obtained by integral transformation can truly reflect whether the frequency exists in the signal and avoid false periodic peaks. Zhong et al. [23] verified the superiority of the marginal spectrum in application through the measured vibration signal.

In view of the characteristics of VMD and Hilbert transform, the cyclical analysis method based on VMD and Hilbert transform is proposed. It is adopted to research the cyclical change in the tree ring in the Huangdi Tomb and sunspot number. This will provide further evidence for the research of the relationship between the solar activity and climate change.

2. Data and Research Method

2.1. Tree Ring Data

Tree ring width data are derived from the study in [24]. The sample sequence is composed of tree ring width of the cypress in the Xuanyuan Huangdi Tomb (altitude 980 m, slope aspect northwest, and slope 40°) and the cypress in Liu’s Tomb (altitude 930 m, slope aspect south, and slope 15°) 7 km away from the Liyuan Park sampled by the Meteorological Bureau of Shaanxi Province in Huangling Country (35°31′N, 109°12′E) in 1975. The cypress sample is sampled on the upper part of the forest belt, and the section is taken at 30 m above the ground. The sequence length is 256 years (from 1700 to 1955) as shown in Figure 1.

2.2. Sunspot Data

In this paper, the annual average of sunspot [25] is from the Solar Influences Data Analysis Center (http://www.sidc.be/silso/datafiles) during 1700 to 1955, as shown in Figure 2. This sequence only has a length of 256 years, so it is difficult to accurately research the periods of more than 100 years. Therefore, the periods of less than 100 years are only researched.

2.3. Research Method
2.3.1. Variational Mode Decomposition

VMD is an effective signal decomposition method. Its overall framework is a variational problem [2628]. VMD is different from the EMD of circular filtering. Each IMF is assumed to be a finite bandwidth with a different center frequency. Its goal is to minimize the sum of the estimated bandwidths for each IMF. The algorithm can be divided into the structure and solution of the variational problem. The detailed description is as follows [29]:

(1) The Structure of Variational Problem. The VMD algorithm converts the signal decomposition into an iterative solution process of the variational problem. The original signal is decomposed into mode functions, , such that the sum of the estimated bandwidths for each mode function is minimized. Its constrained variational model is constructed as follows:

(2) The Solution of Variational Problem. (a) The above constrained variational problem can be changed into a nonbinding variational problem by introducing a quadratic penalty factor and Lagrange multipliers, where guarantees the reconstruction accuracy of the signal and maintains the rigor of the constraint. The augmented Lagrange is denoted as follows:

(b) The alternate direction multiplier method (ADMM) is used to solve the saddle points of the above variational problem, and then , , and are updated alternately (n represents the number of iterations), which are given bywhere is equal to the Wiener filtering results of the current remaining amount . is the center of gravity of the power spectrum of the current IMF. Inverse Fourier transform on is performed to obtain the real part .

(c) Given a discriminant accuracy > 0, the convergence condition of the stop iteration is as follows:

The specific process of the VMD algorithm is summarized as follows:

Step 1. Initialize , , , and .

Step 2. Update the value of , , and according to equations (3)–(5).

Step 3. Judge whether or not it meets the convergence condition (6). Repeat the above steps to update parameters until the convergence stop condition is satisfied.

Step 4. The corresponding mode subsequences are obtained according to the given mode number.

2.3.2. Hilbert Transform

Hilbert transform has the advantages of good time-frequency resolution and self-adaptation. In the analysis of the nonlinear and nonstationary signal, it can effectively avoid the high-frequency interference result. The marginal spectrum [30] is obtained by integral transformation of the Hilbert spectrum. The existence of energy at a certain frequency in the marginal spectrum indicates that there must be a vibration wave at that frequency. The marginal spectrum is the essential difference between the marginal spectrum, the power spectrum, and the Hilbert spectrum. The marginal spectrum is a major breakthrough in signal analysis, and its calculation process is detailed in [23].

2.3.3. The Proposed Periodic Analysis Method

The cycle analysis method of tree ring and solar activity based on the VMD and Hilbert transform is proposed.(1)The original signal is decomposed into a series of IMFs by VMD(2)The Hilbert transform and variance contribution rate are calculated for the IMF components, and the marginal spectrum is given(3)According to the marginal spectrum, the quasiperiodicity of each modality is extracted

3. Data Analysis

3.1. Multiscale Analysis of the Tree Ring Width

The multiscale decomposition of tree ring width is carried out by VMD, and the result is shown in Figure 3. It can be seen from Figure 3 that the IMF components of 9 different time scales and trend components have been obtained by the multiscale decomposition. Each IMF component has a regular fluctuation around the zero mean, so it is a stationary signal. It can be seen from the decomposition results that tree rings have multiple time scales. The IMF1 component has the minimum frequency, namely, the maximum cycle. The IMF2 component has the largest amplitude and the most regular fluctuation, in which basically the maximum or minimum value keeps appearing once every 10 years. IMF3 and IMF2 are similar in fluctuation, and the latter 4 components have smaller amplitude. Generally, IMF1–IMF9 frequencies are increasingly larger, while the cycles are decreasing in a proper sequence. The amplitude of vibration remained basically unchanged from 1940s to 1950s, which indicates that the ten-year period is the dormant period of local tree growth.

Hilbert transform is performed on the IMF components decomposed by VMD. The marginal spectrum of tree ring width is obtained. In order to highlight the sharpness of the periodic peak, the marginal spectrum shown in Figure 4 is normalized. The peak of the marginal spectrum corresponds to the average period of each IMF component. The reciprocal of abscissa frequency at the peak point is the cycle. The normalized marginal spectrum peaks are all over 0.2. It shows that the periodic components of each mode are real. The corresponding frequencies of the main period peak are 0.039, 0.102, 0.137, 0.174, 0.205, 0.283, 0.346, 0.410, and 0.447 from small to large. Their cycles are 25.6a, 9.8a, 7.3a, 5.8a, 4.9a, 3.5a, 2.9a, 2.4a, and 2.2a, respectively.

The variance contribution rate, correlation coefficient calculation, and sorting condition of IMF components of tree ring width are shown in Table 1. According to Table 1, 9.8a and 7.3a have the largest variance contribution. They are real correlation. Next is the quasiperiodicity of 25.6a, and its variance contribution rate is more than 10%. The remaining components are low-frequency quasiperiodic components lower than 7a. It shows the accuracy of the VMD decomposition. Although the correlation coefficient between these low-frequency components and the original signal is less than 0.3, it is also an indispensable part of the multiscale information analysis of the original tree ring width.

3.2. Multiscale Analysis of the Sunspot Number

The multiscale decomposition of sunspot number is carried out by VMD, and the result is shown in Figure 5. It can seen from Figure 5 that the IMF components of 5 different time scales and trend components have been obtained by the multiscale decomposition. Each IMF component has a regular fluctuation around the zero mean, so it is a stationary signal. It can be seen from the decomposition results that sunspot number has multiple time scales. The IMF5 component has the largest amplitude and the most regular fluctuation, in which basically the maximum or minimum value keeps appearing once every 10 years. The amplitudes of the other components are sequentially weakened. Their frequency increases, and the period decreases gradually.

Hilbert transform is performed on the IMF components decomposed by VMD. The marginal spectrum of sunspot number is obtained and normalized, as shown in Figure 6. The normalized marginal spectrum peaks of the sunspot number are all over 0.2. The corresponding frequencies of the main period peak are 0.0096, 0.0192, 0.0923, 0.1006, and 0.1198 from small to large. Their cycles are 102.1a, 52.2a, 10.8a, 9.9a, and 8.6a, respectively.

The variance contribution rate, correlation coefficient calculation, and sorting condition of IMF components of the sunspot number are shown in Table 2. According to Table 2, 10.8a has the largest variance contribution rate up to 54.53%, and the correlation with the original data is also as high as 0.74. It is significant correlation. The periodic components of 8.8a and 9.9a might belong to the same cycle of 10.8a. Although 52.2a and 102.1a are only slightly related to the sunspot number original signal, the variance contribution rates are 6.54% and 10.15%. It indicates that the low-frequency components are also very important and cannot be ignored.

4. Discussion

The tree ring width of Huangdi Mausoleum in Shaanxi Province and sunspot number have nonlinear and nonstationary characteristics, so the selection of the method is very important for the periodic analysis. In this paper, the VMD and Hilbert transform are adopted to make analysis of tree ring width and sunspot number in the same period, which obtained the multiscale time-series change rule. The tree ring width of the Huangdi Tomb in Shaanxi Province is decomposed into 9 IMF components, and there exists quasiperiodicity of about 25.6a, 10a, and 2–7a. The sunspot number in the same period is decomposed into 5 IMF components, and there exists quasiperiodicity of about 102.2a, 52.2a, and 11a.

The 11a cycle obtained by multiscale analysis of sunspot number by VMD and Hilbert transform is the famous Schwabe [31] cycle, and the variance contribution rate is as high as 57.91%. It is the most significant. The average cycle of IMF5 is 102.1a and close to the Gleissberg cycle [32]. The average cycle of IMF4 is 52.2a. This is called the double Hale cycle. The 9.9a and 8.3a cycle components might belong to the same cycle of 11a. Le and Wang [33] obtained cycles of 11a, 53a, and 101a by analyzing the sunspot number through the wavelet. The research results are more accurate than the previous research results [33].

The 25a cycle of tree ring width in the Huangdi Tomb is calculated by VMD, and the 2–7a cycles of the tree ring are generally regarded as related to the El Niño-Southern Oscillation (ENSO), so there is no further discussion here. For the highly correlated 10a quasiperiodic-mode component, it is the focus of research. It is basically consistent with the Schwabe cycle. They are the main scale components of the original sequence, and the variance contribution rate and correlation are the highest.

Moreover, Liu et al. [34] obtained an average period of precipitation with 10a in Shaanxi Province by wavelet. Tang [35] showed the correlation between precipitation in Shaanxi Province and sunspot by EEMD. Generally, the tree ring reflects the variations of local precipitation [36]. Precipitation is affected by sunspot. Therefore, the tree growth is indirectly influenced by sunspot.

5. Conclusions

In this paper, the cycle analysis method based on VMD and Hilbert transform is proposed. This method is used in periodic analysis for tree ring width and sunspot number in the same period, which provides a new method for characteristics analysis of the nonlinear climate signal. Through analytical comparisons, the proposed method demonstrated its superiority and the following contributions:(1)Tree ring width of the Huangdi Tomb in Shaanxi Province from 1700 to 1955 has the quasiperiodicity of 25.6a, 10a, and 2–7a. The sunspot number from 1700 to 1955 has the quasiperiodicity of 102.1a, 52.2a, and 11a.(2)The major cycles of tree ring width of the Huangdi Tomb in Shaanxi Province and sunspot number have the consistency for 10a, which has provided further evidence for the tree growth being driven by solar activity.(3)Besides the 10a cycle of the tree ring width of the Huangdi Tomb, there are other cycles different from the sunspot number. It indicates that, in addition to the solar activity, the tree ring growth is also influenced by other factors.(4)The proposed method is a cycle analysis method which is suitable for nonlinear and nonstationary signals. It can calculate the real and effective signal period. The proposed method is more accurate than the traditional spectral analysis method.

Data Availability

The annual average of the sunspot number time series in this paper comes from the Solar Influences Data Analysis Center of the Royal Observatory of Belgium (http://sidc.oma.be/silso/datafiles). It records the sunspot data from 1700 to the present. In particular, the data analysis center significantly revised the number of sunspots on July 1, 2015. The data in this article are derived from the data set before revision. The tree ring width data used to support the findings of this study have not been made available because they are part of an ongoing research project and its data have not yet been made public.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (No. 51709228).