Cycle Analysis Method of Tree Ring and Solar Activity Based on Variational Mode Decomposition and Hilbert Transform

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.*en, 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.*e 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. *is shows that sunspot activity is an important factor affecting tree ring growth, and tree ring width is influenced by other external environments.


Introduction
e 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 [2][3][4] have shown that solar activity is closely related to natural environment changes such as temperature and various dust storms, floods, and droughts.erefore, 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].erefore, the correlation between tree ring and sunspot activity has important significance for people's life.
Many scholars [9][10][11][12][13] 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.e 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.e application of convolution in the calculation process of wavelet will lead to mutation at the end point, and the wavelet does not have selfadaptability [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 endpoint flying wings.e 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].e 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.e 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.is will provide further evidence for the research of the relationship between the solar activity and climate change.

Data and Research Method
2.1.Tree Ring Data.Tree ring width data are derived from the study in [24].e 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.
e cypress sample is sampled on the upper part of the forest belt, and the section is taken at 30 m above the ground.e sequence length is 256 years (from 1700 to 1955) as shown in Figure 1.

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. is sequence only has a length of 256 years, so it is difficult to accurately research the periods of more than 100 years.erefore, the periods of less than 100 years are only researched.

Variational Mode Decomposition.
VMD is an effective signal decomposition method.Its overall framework is a variational problem [26][27][28].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.e algorithm can be divided into the structure and solution of the variational problem.e detailed description is as follows [29]: (1) e Structure of Variational Problem.e VMD algorithm converts the signal decomposition into an iterative solution process of the variational problem.
e original signal is decomposed into K mode functions, u k (t) (k � 1, 2, . . ., K), such that the sum of the estimated bandwidths for each mode function is minimized.Its constrained variational model is constructed as follows: (2) e Solution of Variational Problem.(a) e above constrained variational problem can be changed into a nonbinding variational problem by introducing a quadratic penalty factor C and Lagrange multipliers, where C guarantees the reconstruction accuracy of the signal and maintains the rigor of the constraint.
e augmented Lagrange is denoted as follows: (b) e alternate direction multiplier method (ADMM) is used to solve the saddle points of the above variational problem, and then u n+1 k , ω n+1 k , and θ n+1 are updated alternately (n represents the number of iterations), which are given by where  u n+1 k (ω) 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  u k (ω) is performed to obtain the real part u k (t)  .(c) Given a discriminant accuracy e > 0, the convergence condition of the stop iteration is as follows: e specific process of the VMD algorithm is summarized as follows: 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. e corresponding mode subsequences are obtained according to the given mode number.

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.e marginal spectrum [30] is obtained by integral transformation of the Hilbert spectrum.e existence of energy at a certain frequency f in the marginal spectrum indicates that there must be a vibration wave at that frequency.e marginal spectrum is the essential difference between the marginal spectrum, the power spectrum, and the Hilbert spectrum.e marginal spectrum is a major breakthrough in signal analysis, and its calculation process is detailed in [23].

e Proposed Periodic Analysis Method.
e cycle analysis method of tree ring and solar activity based on the VMD and Hilbert transform is proposed.
(1) e original signal is decomposed into a series of IMFs by VMD (2) e 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

Multiscale Analysis of the Tree Ring
Width. e 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.e IMF1 component has the minimum frequency, namely, the maximum cycle.
e 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.e 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. e 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.e peak of the marginal spectrum corresponds to the average period of each IMF component.
e reciprocal of abscissa frequency at the peak point is the cycle.e normalized marginal spectrum peaks are all over 0.2.It shows that the periodic components of each mode are real.
e 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.eir cycles are 25.6a, 9.8a, 7.3a, 5.8a, 4.9a, 3.5a, 2.9a, 2.4a, and 2.2a, respectively.e 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.ey are real correlation.Next is the quasiperiodicity of 25.6a, and its variance contribution rate is more than 10%.e 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.

Multiscale Analysis of the Sunspot Number.
e 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.e 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.e amplitudes of the other components are sequentially weakened.
eir frequency increases, and the period decreases gradually.
Hilbert transform is performed on the IMF components decomposed by VMD. e marginal spectrum of sunspot number is obtained and normalized, as shown in Figure 6.
e normalized marginal spectrum peaks of the sunspot number are all over 0.2.e corresponding frequencies of the main period peak are 0.0096, 0.0192, 0.0923, 0.1006, and 0.1198 from small to large.eir cycles are 102.1a,52.2a, 10.8a, 9.9a, and 8.6a, respectively.e 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.e 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 lowfrequency components are also very important and cannot be ignored.

Discussion
e 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.
e 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.e sunspot number in the same period is decomposed into 5 IMF components, and there exists quasiperiodicity of about 102.2a, 52.2a, and 11a.e 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.e average cycle of IMF5 is 102.1a and close to the Gleissberg cycle [32].e average cycle of IMF4 is 52.2a. is is called the double Hale cycle.e 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 4 Advances in Meteorology the wavelet.e research results are more accurate than the previous research results [33].e 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.ey 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.erefore, the tree growth is indirectly influenced by sunspot.Advances in Meteorology

Conclusions
In this paper, the cycle analysis method based on VMD and Hilbert transform is proposed.is 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.rough analytical comparisons, the proposed method demonstrated its superiority and the following contributions: (    Data Availability e 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.e data in this article are derived from the data set before revision.e 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.

Figure 3 :
Figure 3: IMF components of tree ring width.

Figure 4 :
Figure 4: Normalized marginal spectrum of the tree ring.

Figure 6 :
Figure 6: Normalized marginal spectrum of the sunspot number.

Table 1 :
Variance contribution rate, correlation coefficient calculation, and sorting condition of IMF components of tree ring width.

Table 2 :
Variance contribution rate, correlation coefficient calculation, and sorting condition of IMF components of the sunspot number.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) e proposed method is a cycle analysis method which is suitable for nonlinear and nonstationary signals.It can calculate the real and effective signal period.e proposed method is more accurate than the traditional spectral analysis method.