Non-Parametric Wavelet Functional Analysis for Horizontal and Vertical Displacements Derived from GPS Stations in Western Alaska during the Year 2012

In order to analyze the dynamic processes of the Earth interior and the effect of the propagation of the seismic waves to the surface, a comprehensive study of the Earth crust kinematics is necessary. Although the Global Positing System (GPS) is a powerful method to measure ground displacements and velocities both horizontally and vertically as well as to infer the tectonic stress regime generated by the subsurface processes (from local fault systems to huge tectonic plate movements and active volcanoes), the complexity of the deformation pattern generated during such movements is not always easy to be interpreted. Therefore, it is necessary to work on new methodologies and modifying the previous approaches in order to improve the current methods and better understand the crustal movements. In this paper, we focus on western Alaska area, where many complex faults and active volcanoes exist. In particular, we analyze the data acquired each 30 seconds by three GPS stations located in western Alaska (AC31, AB09 and AB11) from January 1, 2012 to December 31, 2012 in order to compute their displacements in horizontal and vertical components by vectorial summation of the average daily and annual velocities components. Furthermore, we design non-parametric DMeyer and Haar wavelets for horizontal and vertical velocities directions in order to identify significant and homogenous displacements during the year 2012. Finally, the non-parametric decomposition of total horizontal and vertical normalized velocities based on level 1 and level 2 coefficients have been applied to compute normal and cumulative probability histograms related to the accuracy and statistical evolution of each applied wavelet. The results present a very good agreement between the designed non-parametric wavelets and their decomposition functions for each of the three above mentioned GPS stations displacements and velocities during the year 2012.


Introduction
Alaska southern edge abuts the boundary between the North American Plate and the Pacific Plate.The Pacific Plate is moving northwestward in subduction at the annual average rate of 5-7 millimeters with respect to the North American Plate (Altamimi et al. 2002;Suito et al. 2003;Altamimi et al. 2007;AEIC 2008;Altamimi et al. 2011).
The area of Alaska is formed by several active faults, like Fair-weather fault, Denali fault, Castle mountain fault, Tintina fault and Totschunda fault, whose activity is connected to the generation of many big earthquakes (Freymueller et al. 2008); for example surrounding Denali fault, several events with a magnitude higher than 6 occurred over the last century such as the Mw=7.2 on July 7, 1912, the Mw=6.2 on August 31, 1958, the Mw=6.7 on October 23, 2002and the Mw=7.9 on November 3, 2002(USGS Report 2013).
In recent studies, for example in Freymueller et al. (2008), a complete set of information about Alaska tectonic regime from 1992 to 2007 has been achieved by using ground GPS stations data.Most of the recent studies have been concentrated in southern Alaska where the spatial distribution of GPS sites is quite denser than the other parts, especially in comparison to western Alaska area.Recently, it has been found that the velocity pattern of the stations located across southern and central Alaska is spatially complex because of the influence of several significant contributions to the surface deformation field (e.g.Brocher et al. 1994;Fuis et al. 2008;Ali and Freed 2010).
In western Alaska, most of the continuous ground GPS sites located near the Bering Sea coast and the Bering Sea islands, show southward to south-southwestward velocities of about 5 mm/yr, while the sites placed along the Bering Sea coast show a southeastward motion of ~10 mm/yr, which is in consistent to the trenchward motion of Cook Inle Region (Kogan et al. 2002;Fournier and Freymueller 2007;Freymueller et al. 2008).However, Cross and Freymueller (2008) showed that the displacements of the Bering Sea sites are in a good agreement with the clockwise rotation of the rigid Bering plate located in eastern Asia while the direction of the displacements of the sites located inside St. Paul Island, at the southernmost part of the Bering Sea sites, are in the same directions of the velocities in the Sanak segment of the Alaska Peninsula as well as to the rigid Bering plate.
Volcanic deformation plays also a critical role in driving the kinematics of western Alaska area.Okmok volcano generated up to 50 cm of displacement during two main inflation periods in 2002-2003and 2004(Lu et al. 2000a, 2005;;Miyagi et al. 2004;Fournier and Freymueller 2007;Fournier et al. 2009).At Akutan volcano the deformation was essentially related to the intrusion of a large dike in 1996 and to its after effects (Lu et al. 2000c) as well as the volcanic deformation retrieved in Westdahl volcano, at the western end of Unimak Island (Lu et al. 2000b;Mann and Freymueller 2003).Moreover, the complex seismicity pattern requires accurate methods in order to infer both horizontal (the combination of North-South and East-West directions) and vertical (upward/downward direction) displacements during coseismic periods (Moradi et al. 2014;Moradi et al. 2015;Khamespanah et al. 2016).
In this paper, we use the data of three GPS stations located in western Alaska (AC31, AB09 and AB11) acquired from January 1, 2012 to December 31, 2012 (366 days).These selected GPS stations are equipped with continuous double-frequency GPS receivers and collect data based on 1 epoch per 30 seconds.By using the vectorial summation of the horizontal and vertical components of displacements, we compute the total horizontal and vertical daily displacements and the average daily and annual velocities.Next, we design and apply two non-parametric wavelet functions including DMeyer wavelet for the total horizontal displacement and Haar wavelet for the total vertical displacement.In addition, non-parametric decomposition of the total vertical and horizontal normalized displacements based on level 1 and level 2 coefficients are derived through wavelet reconstructed approximation signals at level 1 and wavelet reconstructed detail at level 1 for the vertical displacement.
Moreover, the histogram of the normal probability distribution and the corresponding cumulative histogram for each applied wavelet signal, such as vertical displacement, reconstructed signal approximation and reconstructed signal detail at level 1, are showed and compared in terms of statistical values and uncertainties.
Finally, the accuracy of the signals coming from concordant shakes (sudden common fluctuations of the displacements of the GPS stations) in western Alaska (especially in level number 1 for some big shocks during the year 2012) is discussed in detail.
According to the above discussion, Section 2 describes theory and methodology; used data and observations are presented in Section 3; Section 4 includes results and discussions; finally, summary and conclusions are given in Section 5.

Theory and Methodology
Figure 1 shows the methodology proposed in this study.The flowchart consists of two main core divisions: displacements computations and wavelet analysis; these two parts are discussed in details in the next paragraphs respectively.

Displacements Computations
To compute East-West, North-South and vertical components of displacement and then the total horizontal and vertical displacements daily, the average conventional regression solution based on average daily epochs is proposed.
It is assumed that the displacements of GPS stations satisfy the normal function distribution characteristics (more information about this are available in Ricker (1953), Antoniadis and Pham (1998), Antoniadis and Gregoire (1999) and Chen et al. (2012)).Moreover, the average annual velocities related to each GPS station can easily be computed by multiplying the average daily velocities by 365.25.
In order to calculate the absolute value of the total velocity, the length of the velocity vector is computed by the quadratic sum of each planar displacement element.These total horizontal and vertical velocities, related to each station, are used as inputs into wavelet decomposition and reconstruction analysis functions.
Figure 1.Flowchart of the proposed methodology to analyze and compare GPS stations displacements for seismic analysis in western Alaska

Wavelet Analysis Approach
The wavelet analysis and transform are mostly based on a main function while translation and dyadic dilation of constitute an orthonormal basis of (where Lp(R)q belongs to the Lebesgue spaces; see more in Dunford and Schwartz (1958)) and thus the main term (baby equation in general wavelet theory) for the discrete wavelet transform (DWT) is: Where n and m are integers and the main wavelet function (ψ) is subjected to the Equation (2) as a compulsory united condition (Daubechies 1988;Mallat 1989Mallat , 1999;;Donoho et al. 1995;Acharya and Ray 2005).
Generally, there are several forms of wavelet functions or even a possibility to design a new wavelet function according to the desired application.Most of the existing wavelet functions have their own benefits, properties and applications (Keller 2004;Lilly and Olhede 2012;Wang and Gao 2013;Najibi and Abedini, 2013): harmonic wavelets, for example, have a special spherical harmonics expansion, while triangulation based wavelets apply DWT on a polyhedral approximation for the assumed sphere.
The most important property of the wavelets is a time-frequency localization of their basis, like a windowed Fourier transform which allows time-frequency ability for data representation (Najibi et al. 2017), but the result of the DWT, as decomposition of a superimposed time series Rt, is a matrix of wavelet coefficients as Cmn where m is the time index and n is the scale (or layer) index.Normally, the wavelet coefficients for each layer are similar to localized time series with the sampling interval of Δtn=2nΔt, while Δt may vary due to the kind of distribution status of original data.
By considering the detail of wavelet coefficients (e.g. at level 1 and level 2), the probability distribution function (PDF) of the wavelet coefficients C(t) is a distribution function symmetric to 0 and thus it is applicable to the wavelet superimposed decomposition function in DWT domain (Vidakovic 1999;Keller 2004;Wang et al. 2012;Najibi and Jin, 2013) but, in order to apply the wavelet superimposed decomposition function to the data (which in this study is the total horizontal and vertical velocity from the horizontal and vertical components of the GPS stations displacements), the significant seismic correlation wave among these stations is estimated as a function of the wavelet layer index, which in this study represents the different major frequency bands of the data caused by recent earthquakes in western Alaska.In the following sections, the main designed DWTs, including non-parametric Haar wavelet (section 2.2.1) and DMeyer wavelet (section 2.2.2) are described separately.In particular, we applied a non-parametric Haar wavelet function for total vertical velocities and a DMeyer wavelet function for total horizontal velocities and reconstructed the signals of the GPS stations displacements.

Non-parametric Haar Wavelet
The main wavelet function of Haar (mother wavelet function) (ψ(t)) and its corresponding scaling function (ϕ(t)) for DWT can be expressed as follows (Haar 1910;Chui 1992a;Cohen 1995): (3 where t is the dependent variable. There are many valuable benefits of using a non-parametric Haar wavelet; since by means of Haar DWT, any continuous function defined on the interval [0, 1] can be approximated uniformly on [0, 1] frequency domain by linear combinations of the constant functions of 1, ψ(t), ψ(2t), ψ(4t), …, ψ(2nt) and their shifted functions.
Moreover, the wavelet scaling functions with different scale n have a functional relationship (Equation ( 5)) as well as they follow the coefficients of scale n, which can be calculated by coefficients of scale n+1 (Equation ( 6)) as given by (Chui 1992b): (5) where ϕ and ψ denote the wavelet mother function and the corresponding scale factor respectively and t is the dependent variable in the integral computation.
In other words, if Equation ( 5) satisfies the Haar wavelet main function and the corresponding scale factor, then Equation ( 6) is being satisfied; n and k are the integers from Z numerical set and χ and X are defined as: (7)

Non-Parametric DMeyer Wavelet
In general, the discrete Meyer (DMeyer) wavelet function is an orthogonal wavelet which is indefinitely differentiable with infinite support and defined in frequency domain in terms of function ν which is in turn defined as (Donoho et al. 1996;Keller 2004): Thus DMeyer mother and scale wavelet function are given, based on Equation ( 8), as: where ν is defined in Equation ( 8) and ω is the frequency value of the inserted signal into the wavelet function.
By changing the auxiliary function, you get a family of different wavelets.The most important required properties of the auxiliary function ν (Equation ( 8)), is to ensure the orthogonal analysis, while the wavelet mother function (ψ) does not have finite support, but it decreases to 0 when x→∞; as: Therefore, it is straightforward to apply DMeyer wavelet function into each supposed signal, as here we use the total horizontal velocities of GPS stations in order to recognize the main shakes of the GPS stations, particularly during a destructible earthquake (Moradi et al. 2016).

Data and Observations
Half-minute observations of three continue GPS stations located in the western Alaska have been selected for this study.Figure 2, shows the locations of the selected GPS stations.These stations are installed by UNAVCO and maintained by the Plate Boundary Observatory (PBO) and the International GNSS Service (IGS) (http://www.igscb.jpl.nasa.gov/)for entire Alaska region.
The displacements of the North-South, East-West and vertical components, with their corresponding standard deviations as daily time data series for the above mentioned GPS stations (named AC31, AB09 and AB11), have been computed in a post-processing procedure (e.g. in Ehsani et al. 2016 andNajibi andDevineni, 2017).
In this study, we have considered one year of daily observations from January 1, 2012 to December 31, 2012 (366 days) for those GPS stations located in western Alaska area and aligned in East-West direction; these stations are located in a remote distant area in comparison to the stations located in central or southern Alaska, and close to the American tectonic plate boundary.The signals coming from these stations may then help to trace and investigate the response of the ground in the western part of Alaska which has not been studied sufficiently yet; Table 1 presents their coordinates.

Displacement Computations
According to Fig. 1, which describes the methodology proposed in this study, we use ground GPS stations observations (epochs) as the input data for the displacement computations (Delavar and Najibi, 2009;Piryonesi and Tavakolan, 2017).The time series displacements of the selected stations, AC31, AB09 and AB11, with respect to the North American reference frame (defined as Stable North American Reference Frame (SNARF)) are then presented in Fig. 3.  average daily and annual velocity of the mentioned GPS stations in horizontal and vertical directions (Jin and Najibi, 2014a;2014b), where the negative values denote the corresponding movements of the GPS station downwards.Moreover, the total GPS station velocities, as vertical and horizontal summation of the components are also computed (e.g.see similar examples in Ehsani and Afshar, 2011a;2011b); this is shown in Fig. 4 together with the daily trend.In general, the distribution of the seismic waves, particularly after a very large magnitude and shallow earthquake, is complex and not completely solved yet.This is due to large variations in the characteristics of the Earth's layers and materials (Najibi and Sheibani, 2013).Although the convergence between the Pacific Plate and the North American Plate in the last few decades in Alaska region has been deeply documented (Altamimi et al. 2002 and2011;Freymueller et al. 2008), and shows that Alaska is tectonically very active; however, the major earthquakes reported by National Earthquake Information Center (NEIC) occurred in an area relatively far from the GPS stations selected in this study.According to the NEIC (http://earthquake.usgs.gov/regional/neic/)database, the following earthquakes have taken place throughout entire Alaska State, during the year 2012: • August 10; Mw=6.2, 96 km southeast of Nikolski, Alaska, at 13.1 km b.s.l.
• November 12; Mw=6.4,near Gulf of Alaska, the United States, at 12 km b.s.l.Notwithstanding this, our computations (Fig. 3 and Fig. 4) show that moderately large shocks have also occurred in western Alaska region during 2012 on March 24, September 17 and December 29.These events correspond to large displacements of the selected GPS stations, especially in the total horizontal and vertical direction for the year 2012.Fig. 5 shows the exact value of horizontal and vertical displacements caused by the three aforementioned shocks.Referring to the comparison between the horizontal and the vertical components, it can be understood from Fig. 5 that the total horizontal displacements (dark green, red and blue colored bar-column) are very large on March 24, while the total vertical displacements (light green, red and blue colored bar-column) are very large on December 29.According to the dynamics of Alaska area, (e.g. in Altamimi et al. 2011), the Pacific Plate subducts under the North American Plate, enhancing the vertical and North-South displacements of the ground GPS stations, particularly for those located in western Alaska region; displacements in East-West direction seem not to be affected by the dynamics of the subduction process.

Wavelet Configuration for Total Horizontal and Vertical GPS Station Displacements
The total vectorial summation of horizontal (North-South and East-West) and vertical displacements components have been used as the input signals for the corresponding designed horizontal and vertical wavelet functions; as discussed in details in Section 2, the DMeyer wavelet function for total horizontal velocities and the non-parametric Haar wavelet function for total vertical velocities, together with their relevant reconstructed signals of GPS stations velocities, are given in Fig. 6 and Fig. 7 respectively.Considering Fig. 6 and Fig. 7 together, we can see that the sensitivity of Haar mother wavelet function is in good agreement with the vertical velocities variations, while the DMeyer wavelet function is affected completely by the horizontal velocities variation, especially in the days 84, 261 and 364 of the year 2012.The wavelets are used in two level numbers of normalized coefficients, including level 1 and level 2, where the yellow color stands for the sudden changes in the total horizontal and vertical direction, while the red color represents the static case and fewer movements originally sensed by the GPS station.Moreover, the panels (c) in Figs. 6 and 7 show the first and second number of coefficient levels in the corresponding designed wavelets of horizontal and vertical displacement of the GPS station signals and the light-green and dark-green graphs represent these two levels only for the normalized wavelets decomposition of the horizontal and vertical displacement signals (Najibi et al. 2013a;2013b;Najibi and Jin, 2015).together with a precise applied scale in two levels as level 1 and 2 (Najibi and Abedini, 2010;Jahandideh et al. 2014).However both the designed wavelets of DMeyer and Haar proposed at level 2 are more sensitive than level 1, since a small irregularity in the signal would affect it largely, considering the similar situations of the wavelet functions.However, the designed wavelets at level 1 are more stable when facing erratic movements and may neglect high abnormalities at a certain time in a certain time interval.Total Horizontal (H) . Non-parametric DMeyer wavelet function and its normalized wavelet in Level 1 and 2 for the total horizontal GPS stations velocity.The standard deviations for the vertical displacements of the selected GPS stations are very large; however, the standard deviation for the designed Haar wavelet reconstructed signal approximation decomposition in level 1 is bigger than the reconstructed signal details at the same level 1.This can also be understood from Fig. 7, where the color bars in both level 1 and level 2 are more dispersed.The bias for level 1 is smaller than level 2, but the reason is that level 1, which has been inserted in both the reconstructed approximation signal and the

Summary and Conclusions
In this paper, we analyzed and elaborated the observations (epochs) of three GPS stations located in western Alaska (AC31, AB09 and AB11).From these observations, we computed the components of the horizontal (North-South and East-West) and vertical displacements.Moreover, their corresponding total horizontal and vertical displacements variations, derived from vectorial summation of the displacements during 366 days starting from January 1, 2012 to December 31, 2012 have been also computed.Then, basing on the total horizontal and vertical components of the above mentioned GPS stations, the average daily velocity and the annual velocity of each station have been computed during 2012.Furthermore, we applied non-parametric DMeyer wavelets for the total horizontal velocity signal and a Haar wavelet for the total vertical velocity signal.Besides, non-parametric decomposition of the total vertical and horizontal normalized displacements, based on two levels, i.e. level 1 and level 2 coefficients, were achieved by using wavelet reconstructed approximation signals at level 1 and wavelet reconstructed detail at level 1 for the vertical displacement of the GPS stations.The normal histograms probability distribution and the corresponding cumulative histograms for each applied As it presented in Fig. 4 and Fig. 5, the homogenous sudden fluctuations of each GPS station, in both total horizontal and vertical velocity directions, which are caused by high magnitude earthquakes during the days 84 (March 24), 261 (September 17) and 364 (December 29) of 2012 over western Alaska region, are consistent to the relevant designed wavelets in terms of maximum vibration variations, as showed in Fig. 6 and Fig. 7 respectively, in both level 1 and level 2. The normal and cumulative histograms (e.g.Fig. 8, 9 and 10) demonstrate the power of the effective designed wavelets to retrieve sudden changes in the signal, especially for the reconstructed detail signal at level 1 for each GPS station, where the overall trend of the histogram probability distribution converges to a normal distribution rather than the vertical displacement signal or the reconstructed approximation at level 1.
Practically, the wavelet transforms (both discrete and continuous) can localize information in both time and frequency space simultaneously; and therefore it carries out a powerful technique for seismic data analysis.However, the designed wavelets have very good capabilities to retrieve the seismic waves signals for all the selected GPS stations during the year 2012, the Earth's surface vibration modes, combined with the response from the interior layers, are very complex to be examined through non-parametric solutions; then further studies and researches are necessary, in particular for western Alaska region where the seismic movements are accommodated by the existence of several faults and where the presence of volcanic activity implies a complex kinematics.
However, the good results of our approach stimulate us to look forward for the upcoming high resolution data of space geodesy and based remote sensing facilities, such as Synthetic Aperture Radar (SAR) and as aerial laser scanner space sensors, like the Light Detection And Laser Ranging (LiDAR), with whom a further better understanding of the Earth's crustal movements can be achieved.

Figure 2 .
Figure 2. Location of GPS stations in western Alaska (this study)

Figure 3 .
Figure 3. Vertical, East-West and North-South displacements of continuous GPS stations AC31, AB09 and AB11 and their corresponding standard deviations; especially the sudden fluctuation in different GPS station displacement components on March 24 (day: 84), September 17 (day: 261) and December 29 (day: 364)

Figure 4 .
Figure 4. Total (vectorial) GPS station velocities computed as horizontal summation (left side column) and vertical summation (right side column) of North-South and East-West components as well as the yearly trend (red dashed-line) for each GPS station

Figure 5 .
Figure5.Values of the total displacements computed as horizontal (H) and vertical (V) direction for three big shocks during the days 84, 261 and 364 of 2012 for AC31, AB09 and AB11 GPS station; the subsiding seismic wave for the day 84 resulted in negative values for all these three GPS stations homogenously

Fig. 6
Fig.6and Fig.7demonstrate the ability of the designed wavelet functions to represent the signals of the sudden fluctuations throughout the surface through horizontal and vertical displacements; these fluctuations are represented by the colored bar variations, among which the yellow colored bars represent the bigger fluctuations, together with a precise applied scale in two levels as level 1 and 2(Najibi and Abedini, 2010;Jahandideh et al. 2014).However both the designed wavelets of DMeyer and Haar proposed at level 2 are more sensitive than level 1, since a small irregularity in the signal would affect it largely, considering the similar situations of the wavelet functions.However, the designed wavelets at level 1 are more stable when facing erratic movements and may neglect high abnormalities at a certain time in a certain time interval.

Figure 7 .
Figure 7. Non-parametric Haar wavelet function and its normalized wavelet in Level 1 and 2 for the total vertical GPS stations velocity

Figure 8 .
Figure 8. Histogram probability for the vertical displacement and corresponding cumulative histogram for AC31, AB09 and AB11 GPS stations during the year 2012 for the vertical displacements of AC31, AB09 and AB11, in the selected area in western Alaska, has less sensitivity to the noise and thus level 2 deviates more from the mean value of zero.

Figure 9 .
Figure 9. Histogram probability for displacement reconstructed approximation at level 1 and corresponding cumulative histogram for AC31, AB09 and AB11 GPS stations during the year 2012

Figure 10 .
Figure 10.Histogram probability for displacement reconstructed detail at level 1 and corresponding cumulative histogram for AC31, AB09 and AB11 GPS stations during the year 2012 , such as vertical displacement, reconstructed signal approximation and reconstructed signal detail at level 1, have been derived and discussed in terms of statistical values and uncertainties comparisons.

Table 1 .
Positions of the GPS stations used in this study (western Alaska)