The Role of Polarization Analysis in Reducing Natural Hazard

Abstract

The complexity of the regions where the landslide occurred requires a detailed survey of the site response properties of these regions. Ambient noise analysis is a common and nondestructive approach that provides more detailed information on site resonance properties characterized by directional variations. Obtaining geological information from seismic data motivates researchers to innovate and improve efficient tools for seismic wave processing. Polarization-based methods have received much attention regarding their capability to discriminate between different phases of the seismic wave based on their polarities. Using polarization filtering to extract Rayleigh wave ellipticity provides more detailed information on site resonance properties characterized by directional variations. The outputs showed better performance of the method in terms of stability and reliability of the results. Indeed, it enabled the detection of site resonance characteristics that were previously undetectable by classic Nakamura method.

1. Introduction

The annual destruction caused by landslides inspires near-surface engineers to discover strategies to prevent them. Among the mechanisms responsible for slope instability, earthquakes require special attention due to their ability to trigger a large number of simultaneous mass movements, complicating the management of seismic crises. Slopes susceptible to landslides exposed to earthquakes can be more easily mobilized in the presence of ground motion amplification due to subsoil physical properties [1] or topography [2], characterized by directional variations of resonance phenomenon. The impacts of directional resonance on soil slopes have been investigated by analyzing seismic ground motion at landslide-prone slopes [3]. For instance, Burjánek et al., in [4] found that the largest ground motion amplification was directed along the steepest slope in a survey conducted at the Randa rock slope in Switzerland.

The standard spectral ratio (SSR) approach, which records seismic events at research sites and a nearby reference site to compare site response characteristics, is the most straightforward method for studying site response. The lack of long-term accelerometers on the rough topography of landslide-prone areas to record ground motion has increased the need for alternative approaches to studying site effects. The high expense and limitations of studying subsurface structure by boreholes and other geophysical methods motivated the researchers to choose ambient noise analysis as a feasible and cost-effective alternative tool for monitoring landslide-prone slopes [5, 6, 7]. The ambient noise analysis provides valuable information regarding the effect of subsurface structure on site response at observation sites [8]. Nakamura method [9], also known by the acronym HVSR (horizontal-to-vertical spectral ratio), is the most popular method for ambient noise analysis. It analyzes the spectral ratios, H/V, between the horizontal and vertical components of ambient noise recordings, disclosing site-specific resonance frequencies via H/V curve peaks. Directional resonance can also be defined by the analysis of azimuthal changes in spectral ratios [10].

It has been shown that Nakamura method cannot precisely define resonance frequencies when amplification occurs at several frequencies, particularly in sites with complex geology [7, 8, 11]. To effectively analyze local site amplification, a complete understanding of local soil characteristics such as shear-wave velocities and sediment thickness is required. Surface wave dispersion curves generated from active or passive seismic array data can provide valuable information on soil properties. However, it often renders reliable information for frequencies higher than the resonance frequency, and the inverted velocity profiles do not give information on deeper structures. Several studies have shown that the S-wave velocity profile may be calculated by inverting the H/V spectra [12, 13] or by inverting the H/V spectra and dispersion curves [13, 14]. Since both Rayleigh and Love waves contribute to the H/V spectrum, the extraction of these waves can be effective in the inversion process [15, 16].

Polarization methods are the efficient approaches that have the ability to discriminate between different phases of the waves and extract Rayleigh waves depending on their polarities [17]. Polarization is a seismic wavefield property that describes the particle motion, such as linear polarization in P-waves or the elliptical particle motion in Rayleigh waves. The ellipticity curve calculated by the polarization methods is more similar to the theoretical curve than the classic Nakamura H/V curve [18]. In the following, some of these polarization analysis methods are explained with the aim of improving the reliability of the results in seismic hazard assessments.

2. Polarization analysis methods

Polarization techniques are mostly used in the time or frequency domains, depending on the wave content and event type. Time-domain approaches are more useful in the event of the near arrival of waves in time, such as wave conversions at crustal discontinuities. The direction of particle motion and rectilinearity of the three-component (3C) earthquake data were defined using the correlation approach in time-domain polarization filtering proposed by Montalbetti and Kanasewich in [19]. Although their nonlinear polarization filter enhanced the shear and compressional phases, it ignored the fact that the waveforms have distinct polarizations for various frequencies. Addressing this problem, Samson in [20] studied particle polarization as a function of frequency, and after that, Samson and Olson in [21] developed a frequency domain polarization filter. In real-world seismograms, events with different frequencies and polarization characteristics may overlap, while events with identical frequencies but distinct polarization content may occur at different times. Therefore, the researchers became interested in time-frequency (TF) methods that consider polarization variations of 3C signals [22, 23, 24, 25, 26].

2.1 Time-domain polarization analysis

According to certain research (e.g., [27, 28]), ambient noise contains a mix of several types of waves (body, Rayleigh, and Love waves). Their proportional percentage is determined by site circumstances and source characteristics, which may differ from one example to the next. As a result, assumptions about this fraction make the inference of the site resonance features from H/V curves more challenging. Generally, extracting Rayleigh waves from a recording can yield the following results: (I) an estimation of the S-wave resonance frequencies; (ii) identification of site response directivity and orientation from preferential polarization directions; and (iii) some constraints for velocity modeling using the ellipticity curves. Since horizontal components of ambient noise can contain a large or even dominating fraction of Love waves (cf. [29]), the H/V ratios can be different from the anticipated one for Rayleigh waves. Therefore, some studies tried to separate Rayleigh waves from Love waves in noise recordings by selecting signal parts with considerable energy in the vertical component [30, 31]. On the same lines of inquiry, Del Gaudio in [11] provided a technique for identifying wave packets with Rayleigh-type particle motion inside ambient noise recordings by identifying instantaneous polarization features.

2.1.1 Methodology

Considering the multicomponent signal, its analytic transformation has the form

$${\mathbf{u}}_c\left(t\right)=\mathbf{u}\left(t\right)+\mathrm{j}\widehat{\mathbf{u}}\left(t\right)=\mathbf{A}\left(t\right)e^{\mathrm{j}\mathbf{\phi }\left(t\right)}$$E1

where u(t) represents the real component of a vector on a complex plane, $j$ is the imaginary unit, and $\widehat{\mathbf{u}}\left(t\right)$ is the Hilbert transform of u(t). The vector modulus A(t) is the instantaneous amplitude, and phase φ(t) is the real phase. Consider a(t) and b(t) as the major and minor semi-axes of the instantaneous elliptical trajectories, respectively: these two vectors can be derived by calculating the phase shifts, ${\mathrm{\Phi }}_0$, as defined in [23].

$${\mathbf{a}}_i\left(t\right)=Re\left[{\mathbf{A}}_i\left(t\right).e^{\mathrm{j}\left({\mathbf{\phi }}_i\left(t\right)-{\varphi }_0\left(t\right)\right)}\right]\mathit{\text{for}}i=E,N,Z$$

$${\mathbf{b}}_i\left(t\right)=Re\left[{\mathbf{A}}_i\left(t\right).e^{\mathrm{j}\left({\mathbf{\phi }}_i\left(t\right)-{\varphi }_0\left(t\right)-\frac{\pi }{2}\right)}\right]\mathit{\text{for}}i=E,N,Z$$E2

In 3D, the orientation of the instantaneous elliptical trajectory can be characterized by the planarity vector, p(t), which is defined as the vector product of a(t) and b(t) [32]. When the a(t) amplitude exceeds the b(t) amplitude, the elliptical trajectory degenerates into a rectilinear ground motion. According to Schimmel and Gallart in [33], the rectilinearity (Eq. 3) can be used to assess the proximity to such a state.

$$\mathbf{rl}\left(t\right)=1-\frac{\left|\overrightarrow{\mathbf{b}}\left(t\right)\right|}{\left|\overrightarrow{\mathbf{a}}\left(t\right)\right|}$$E3

The rectilinearity takes values ranging from 0 to 1, which relate to precisely circular and linear polarizations, respectively.

To define Rayleigh wave packets and their features (ellipticity and polarization direction), an instantaneous polarization analysis can be employed on time series of ambient noise recordings. The analytic transformation can be applied to time series by filtering them through several narrowband filters with varying central frequencies (f_c) to see the variation of Rayleigh wave characteristics with frequency.

The analytic transformation of the two horizontal components permits calculating the maximum amplitude Hmax of ground motion on the horizontal plane and its orientation for each frequency. The analytic transformation of the vertical component yields a measure of its amplitude V, and dividing Hmax by V yields an estimate of the ratio between the amplitudes of the horizontal and vertical components of ground motion, as well as the relative azimuth, at each instant and for the analyzed frequency.

Generally, a considerable dispersion of H/V values may be predicted as a result of the overlapping of various types of waves arriving from a geographically spread distribution of distinct noise sources. However, if Rayleigh waves can be extracted from the noise recording, they should exhibit more coherent characteristics, which can be correlated to site response properties. To detect such signals, the analytic transform is applied to all 3C, yielding the vectors a(t) and b(t), as well as the planarity vector p(t) and the rectilinearity vector rl. Parts of a noise recording containing Rayleigh waves are then recognized based on the presence of a Rayleigh-type polarization of ground motion, that is, if (i) the instantaneous trajectories lie on a vertical plane (which indicates a horizontal p(t)) and (ii) one of a(t) and b(t) is horizontal and the other is vertical. If a(t) is horizontal, H/V will be more than 1; otherwise, H/V will be less than 1.

Since Rayleigh waves can be contaminated by noise effects from other wave types, it is important to consider deviations of p(t), a(t), and b(t), which should not surpass angular thresholds. Theoretically, instantaneous polarization analysis can recognize single Rayleigh wave samples, and the existence of isolated samples matching the identification requirements among a substantial amount of data can be completely coincidental. Therefore, if a coherent form of polarization is observed in more consecutive samples, more reliable detection of Rayleigh waves may be achieved.

Finally, HVIP values derived from Rayleigh-type wave packets (denoted as HVIP_R) are averaged within each filtered time series, and a curve of such HVIP values as a function of filtering central frequency fc represents an estimate of Rayleigh wave ellipticity as a function of frequency. The scatter of instantaneous HVIP_R values around their average represents the uncertainty impacting ellipticity estimations. Analyzing the distribution of HVIP_R azimuths reveals the presence of preferred signal polarization, which might indicate the presence of site response directivity.

2.1.2 Synthetic example

A 900-s synthetic ambient noise recording with a sampling interval of 0.01 s on a model introduced in [34] was generated by the synthetic seismogram codes [35]. The recording was analyzed using classic Nakamura and HVIP methods to compare the H/V curves derived from these methods. The code gpell [36] was applied to calculate the theoretical Rayleigh wave ellipticity through 1D modeling and hence to verify the reliability of the results. Based on the model used in this example, it is expected to observe two main peaks on the H/V curve for the fundamental and first higher mode of Rayleigh waves at 2 and 6 Hz.

Comparing the results represents that both methods detect the fundamental mode of the Rayleigh waves, while only HVIP detects the first higher mode peak. Multiple peaks can be seen on the curves derived from HVIP, which can be related to different vibration modes. It can be inferred that H/V curves obtained by the use of the geometric mean, similar to HVSR method, can underestimate the H/V ratio, and ignore the secondary peaks (Figure 1 and Table 1).

Figure 1.

a) Synthetic H/V curves derived from a) HVIP, and b) Nakamura methods.

	H (km)	Vp (km/s)	Vs (km/s)	RHO (GM/CC)	Qp	Qs
Layer 1	0.025	1.35	0.2	1.9	50	25
Layer 2	0	2.0	1.0	2.5	100	50

Table 1.

Synthetic crustal model.

It is known that Nakamura method can work properly under the condition that both site response and noise source are distributed isotropically. In some cases where noise sources insufficiently energize the frequency band of maximum site amplification, HVIP outperforms HVSR since the former at least can show some reliable peaks in some particular directions, while Nakamura method shows small peaks not only at frequencies that are not significantly amplified but also with the underestimate of amplification. The main discrepancy between HVIP and Nakamura method is that the HVSR averages the H/V values over the length of the recording, in which different types of waves can come along with Rayleigh waves.

The main objective of applying HVIP is to study the site response more accurately due to the ambiguous nature of ambient noise. There are some similarities between the HVIP method and other techniques of polarization analysis. For instance, Jurkevics in [22] applied polarization analysis to site response studies based on the analysis of earthquake recordings. The instantaneous properties of the particle motion are calculated in [37], which show the most similarities to the HVIP approach, to analyze earthquake recordings, while Burjánek et al., performed the method on ambient noise data [4]. Applying polarization analysis simultaneously to the 3C of a recording to identify Rayleigh wave packets makes the HVIP approach different from [4].

2.1.3 Site response directivity

Site response directivity is a directional resonance characteristic that has been shown to accurately reflect resonance frequencies associated with topographic amplification, anisotropy in slope materials, and mechanical features [38]. In studying the site response directivity in landslide areas, Del Gaudio et al., in [8] were the pioneers who presented a comprehensive directional analysis of ambient noise data to study the differences in H/V ratios estimated along orthogonal directions. Despite the numerous case studies in landslide areas, the origin of directional resonance has remained uncertain due to the difficulty in recognizing the involved factors. The anisotropy of shear wave velocity generated by landslide mobilization might represent directed amplification; however, the directional amplification is not always in the direction of the landslide, and other factors, such as topography, geology, lithology, and tectonic features, can be effective in generating this phenomenon [7, 39, 40, 41]. It is important to note that some sites are more likely to experience stronger shaking in a specific direction and frequency range, which is of interest to earthquake engineers. For instance, the similar direction of slope and maximum amplification can amplify the influence of directivity on mass movement [7, 41]. For these reasons, studying site response directivity has become increasingly important in seismic hazard assessments.

2.1.4 Real example

On June 1, 2014, an earthquake with a magnitude of 3.6 followed by five aftershocks occurred on two sides of an active fault in the vicinity of the Nargeschal village, located in the north of Iran, and then a landslide occurred 3 days after this seismic sequence [7]. Based on the landslide hazard zonation map of the province given by the Geological Survey and Mineral Exploration of Iran, Nargeschal village lies in a high-risk zone, which makes landslides the most serious geological hazard in the study region. Slope steepness (nearly 30°) and the presence of soft and sticky soils covering relatively hard and brittle rocks create conditions for gravity-driven motions. Furthermore, the region’s tectonic activity, which retains the dynamic stresses discharged by seismic activity, can operate as both a triggering and predisposing factor for slope failures.

To evaluate the landslide and prevent future damages, ambient noise data were recorded and analyzed by the HVIP method. The results confirmed the role of landslides and structural features in the directional resonance detected in the study region. As an example, the polar diagrams that show the results of HVIP analysis at stations located on the head and foot of the landslide are depicted in Figure 2. The results represent a preferential directivity that is perpendicular to the tension cracks. The reason for the observed directivity is a mass movement, which induces the opening of fissure or fracture systems transverse to the sliding direction and makes the material more deformable in such a direction.

Figure 2.

Polar diagrams depicting the results of the ambient noise analysis processed via the HVIP method for the stations located on a) the head and b) the foot of the landslide region (the black arrow denotes the slope direction, the red and white arrows are the direction of tension cracks and old landslides, respectively).

This example confirms the necessity of an advanced ambient noise analysis technique on sites with complex lithology, where the classic Nakamura method could have difficulty revealing site resonance properties. The detailed HVIP analysis allows us to have a more reliable hazard zonation map to identify the areas that are more susceptible to future seismically induced mobilization. More details about this example can be found in [41].

2.2 TF polarization analysis

TF decompositions have a significant role in designing TF polarization filters. Although several TF decompositions are available in the literature [42, 43, 44, 45, 46, 47], the resolution of these decompositions is remarkably effective in designing polarization filters.

Wavefields are typically recorded in geophysical research using 3C sensors that transfer the same information with various energies. Since real geophysical data are typically contaminated by noise, signals with lower energy may be weakened while being decomposed by high-resolution methods. This is critical in analyzing all 3C of a signal at the same time since the support for any specific part of a wavefield should be consistent in all 3C TF maps and just vary by a scaling factor. Kakhki et al., in [48] proved that employing a single-component sparse adaptive S transform for 3C data results in suboptimal decompositions in the TF domain. Therefore, they proposed the 3C sparse adaptive ST (3C-SAST), which employs the group sparse constraint to lock the TF supports of all 3C of a signal together and retrieves the spectral paths of weak parts of waveforms without deformation. The invertibility and resolution flexibility make the proposed method a viable tool for integrating the TF relationship in polarization analysis. Then, the polarization attributes can be calculated directly from the TF maps derived from the 3C-SAST without using Eigen analysis or other complicated methods [25, 49]. These attributes, which have a higher resolution compared with other available decompositions, can be used to construct an efficient TF filter to remove the undesired waves.

2.2.1 Methodology

Polarization has an intuitive and straightforward description of a monochromatic wave propagating in 3D space. Assume the 3C seismic data, $\left[N\left(t\right),E\left(t\right),Z\left(t\right)\right]$, and the TF transform of these components are $\left[N\left(t,f\right),E(t,f),Z(t,f)\right]$, whose real and imaginary parts are N_R(t,f), N_I(t,f), E_R(t,f), E_I(t,f), Z_R,(t,f), and Z_I(t,f), respectively. The ellipse parameters, schematically depicted in Figure 3, can be calculated directly from the real and imaginary parts of the 3C TF spectra. For the sake of brevity, only the results are presented, and more explanation can be found in [49].

$$a=\frac{1}{\sqrt{2}}\sqrt{A+\sqrt{B^2+C^2}}$$

$$b=\frac{1}{\sqrt{2}}\sqrt{A-\sqrt{B^2+C^2}}$$

$$I=arctan\left\{\frac{{\left[{\left(Z_R{Y}_I-Z_I{Y}_R\right)}^2+{\left(Z_R{X}_I-Z_I{X}_R\right)}^2\right]}^{1/2}}{\left(Y_R{X}_I-Y_I{X}_R\right)}\right\}$$

$$\mathrm{\Omega }=arctan\left\{\frac{\left(Z_R{Y}_I-Z_I{Y}_R\right)}{\left(Z_R{X}_I-Z_I{X}_R\right)}\right\}$$E4

where

Figure 3.

The ellipse of prograde motion is plotted together with the polarization attributes. D and D′ are the intersection of the ellipse with the horizontal plane, a and b are the ellipse major and minor semi-axes, respectively, the angle I is the inclination of the ellipse plane concerning the horizontal plane. The angle Ω is that between DD’ and the x axis, measured on the horizontal plane.

$$A=X_R^2+X_I^2+Y_R^2+Y_I^2+Z_R^2+Z_I^2$$

$$B=X_R^2-X_I^2+Y_R^2-Y_I^2+Z_R^2-Z_I^2$$

$$C=-2\left(X_R{X}_I+Y_R{Y}_I+Z_R{Z}_I\right)$$E5

2.2.2 Earthquake example

To evaluate the role of a high-resolution TF decomposition in calculating the polarization attributes and to represent the efficiency of the TF polarization method in detecting different wave types, we chose an earthquake example that occurred on June 3, 2003, and was recorded at the PEMO station, in Canada (45.68° N, 77.25° W), with a magnitude of 6.9 and an epicentral distance of 7124 km (Figure 4).

Figure 4.

The earthquake recorded at PEMO station. a) Radial, b) transverse, and c) vertical components.

Figure 5 depicts the polarization attributes estimated by employing the 3C-SAST and SAST, which represent the higher resolution of the 3C-SAST in providing polarization attributes in the TF domain. The resolution of the attributes is an important factor in designing the polarization filters since their lower resolution leads to inaccuracies in the distinction between various wave types, particularly in the presence of noise. In this example, the polarization attributes were used to detect Rayleigh and Love waves. For instance, Rayleigh waves can be observed more obviously in the semi-minor axis since this attribute is representative of the circular motion. On the other hand, the dip attribute shows a stable value of π/2 for the Rayleigh waves, while it has an unstable value between 0 and π for the Love waves, which is 0 in this example. Therefore, the initial high-energy arrival has Love wave features, but the second one has Rayleigh wave characteristics. Finally, after detecting the aimed wave type, which is Rayleigh wave in this example, a filter can be designed based on the polarization attributes to extract it.

Figure 5.

The spectra of semi-major, semi minor, and the inclination of the seismogram in Figure 4 by 3C-SAST transform (upper plots) and SAST (lower plots).

3. Conclusions

Evaluating seismically induced vulnerability to slope failure has been a critical problem for geotechnical engineers, especially in metropolitan areas near landslide-prone zones. Slopes with a deep-seated landslide can be an amplification factor for seismic ground motion, with a maximum along the sliding direction. In such cases, slopes with high impedance contrast of layers throughout the 3D geometry of a landslide body may trap seismic waves and prolong shaking, making slopes more susceptible to seismically induced landslide reactivation. This phenomenon is more possible if shaking energy is polarized along the hill slope facing direction. The presence of more factors generating directional spectral peaks is a common characteristic of site response directivity phenomena. Ambient noise analysis is a practical and cost-effective method for studying site response on landslide-prone slopes. Despite the popularity of Nakamura method for ambient noise analysis, it has limitations at sites with complex geology. This limitation can be addressed by extracting Rayleigh waves from ambient noise recordings using polarization methods in the time, frequency, or TF domain. The synthetic and real-world examples proved the superiority of polarization methods over classical methods in studying site response in landslide areas.

Acknowledgments

The authors appreciate the funding for this study from the Fundação de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ) with grant numbers E-204.210/2021, and E-26/201.124/2021, Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) with grant numbers 305949/2019-5, and 422589/2021-7, and Centro de Pesquisas Leopoldo Américo Miguez de Mello (CENPES) under grant number SAP 4600571711.

Conflict of interest

The authors declare no conflict of interest.