Site characterization at downhole arrays by joint inversion of dispersion data and acceleration time series
SSite characterization at downhole arrays by joint inversionof dispersion data and acceleration time series
Elnaz Seylabi ∗ , Andrew Stuart † , Domniki Asimaki ‡ Abstract
We present a sequential data assimilation algorithm based on the ensemble Kalman inversion to estimate the near-surface shear wave velocity profile and damping when heterogeneous data and a priori information that can berepresented in forms of (physical) equality and inequality constraints in the inverse problem are available. Althoughnon-invasive methods, such as surface wave testing, are efficient and cost effective methods for inferring Vs profile,one should acknowledge that site characterization using inverse analyses can yield erroneous results associated withthe inverse problem non-uniqueness. One viable solution to alleviate the inverse problem ill-posedness is to enrichthe prior knowledge and/or the data space with complementary observations. In the case of non-invasive methods,the pertinent data are the dispersion curve of surface waves, typically resolved by means of active source methods athigh frequencies and passive methods at low frequencies. To improve the inverse problem well-posedness, horizontalto vertical spectral ratio (HVSR) data are commonly used jointly with the dispersion data in the inversion. In thispaper, we show that the joint inversion of dispersion and strong motion downhole array data can also reduce themargins of uncertainty in the Vs profile estimation. This is because acceleration time-series recorded at downholearrays include both body and surface waves and therefore can enrich the observational data space in the inverseproblem setting. We also show how the proposed algorithm can be modified to systematically incorporate physicalconstraints that further enhance its well-posedness. We use both synthetic and real data to examine the performanceof the proposed framework in estimation of Vs profile and damping at the Garner Valley downhole array, and comparethem against the Vs estimations in previous studies.
Downhole arrays have been extensively used as testbeds by engineers and earth scientists for validationand improvement of predictive models of site response and physics-based ground motions. Strong motionsrecorded at depth, in particular, are widely sought after boundary conditions for the validation of one-dimensional wave propagation codes, as they minimize the uncertainty associated with source and patheffects when studying problems of site response. To best serve as validation testbeds for equivalent linearand nonlinear site response analyses, however, downhole array sites should also be accompanied by reliableestimates of soil profiles–such as small strain shear modulus (or shear wave velocity), damping and nonlinearsoil properties–as well as their variability.Site characterization efforts at strong motion stations include shear wave velocity measurements at mul-tiple locations, using both invasive (e.g. PS-logging) and/or non-invasive methods (e.g. spectral analysis(Stokoe II, 1994) and multi-channel analysis of surface waves (Foti, 2000)). Nonlinear soil properties, on theother hand, are most frequently estimated from empirical correlations of published laboratory data (e.g.,Darendeli, 2001), or, in rare cases, are measured from undisturbed samples retrieved at the site. Laboratorymeasured properties, however, are not always reliable estimates of in-situ soil properties, an incompatibilityassociated with sample disturbance, measurement error and inherent field-scale spatial variability of soil ∗ Department of Civil and Environmental Engineering, University of Nevada, Reno, NV 89557 ([email protected]) † Department of Computing and Mathematical Sciences, California Institute of Technology, Pasadena, CA 91125 ‡ Department of Mechanical and Civil Engineering, California Institute of Technology, Pasadena, CA 91125 a r X i v : . [ phy s i c s . g e o - ph ] A p r Introduction properties. Downhole array ground motion records have been used instead by several researchers to inferin-situ soil parameters. More specifically, low amplitude motions have been used to constrain the Vs profileand damping (e.g., Assimaki et al., 2006) while high amplitude motions have been used to parameterizenonlinear soil behavior in moderate to large strains (e.g., Assimaki et al., 2011; Chandra et al., 2015; Seylabiand Asimaki, 2018).While it is commonly assumed that invasive methods are more reliable than non-invasive methods to retrieveVs profiles, recent studies have highlighted that the latter are more efficient and cost effective, while yieldinguncertainties comparable to the order of invasive methods (Garofalo et al., 2016; Teague and Cox, 2016).Still, one should acknowledge that site characterization using inverse analyses, that is surface wave testingor downhole array data, can yield erroneous results associated with the inverse problem non-uniqueness.Furthermore, recent studies have shown that significantly different Vs profiles that satisfy the error thresholdsof the inversion process, can still result in similar linear viscoelastic seismic site response and amplificationfactors (e.g., Foti et al., 2009).One viable solution to alleviate the inverse problem ill-posedness is to enrich the prior information and/orthe data space with complementary data. In the case of non-invasive methods, the pertinent data are thedispersion curve of surface waves, typically resolved by means of active source methods at high frequenciesand passive methods at low frequencies. Horizontal to vertical spectral ratios (HVSR) have also been usedextensively to approximate a site’s predominant frequency, and therefore constrain the velocity structure atdepth. Furthermore, the joint inversion of HVSR with dispersion data has been used successfully to improvethe site characterization accuracy (e.g., Arai and Tokimatsu, 2005; Pi˜na-Flores et al., 2016; Molnar et al.,2018; Lunedei and Malischewsky, 2015).At strong motion arrays, recorded acceleration time-series, which include both body and surface waves, canalso be used to complement dispersion data. Because of the complementary characteristics of the body andsurface waves – the former carrying information in the form of travel time and the latter in the form of near-surface dispersive characteristics – we anticipate that formulating the inversion problem using both dispersiondata and acceleration time-series should reduce the margins of uncertainty in the Vs profile estimation.Inversion methods that are used in site characterization include the stochastic direct search method, whichis based on the neighborhood algorithm (e.g., Wathelet et al., 2004); the uniform Monte Carlo method(e.g., Socco and Boiero, 2008), which finds an ensemble of models minimizing the data misfit; and the fullyBayesian Markov Chain Monte Carlo method (e.g., Molnar et al., 2010), which provides the most probableand quantitative uncertainty estimates of the Vs profile. In this paper, we present a framework based onthe ensemble Kalman inversion and use it to examine whether and how the joint inversion of dispersionand downhole array data improves site characterization estimates. We also show how one can systematicallyincorporate a priori information in terms of physical constraints in the ensemble Kalman inversion to improvethe problem well-posedness.In the rest of this paper, we first explain the ingredients of the framework in § Inverse Problem. In § SiteCharacterization–Synthetic Data and § Site Characterization–Real Data, we perform numerical experimentsusing synthetic and real data simulated or recorded at the Garner Valley Downhole Array (GVDA) site– one of the best-studied and best-instrumented sites in Southern California. More specifically, in § SiteCharacterization–Synthetic Data we use synthetic data to study how the combined data sets improve theVs profile estimation relative to the individual data. Next, in § Site Characterization–Real Data, we userecorded acceleration time-series by the array and experimental dispersion data to invert for the Vs profileand damping, and compare our results to the inverted Vs profiles from previous studies. Finally, we provideconcluding remarks in § Conclusion.
Inverse Problem We consider the problem of finding u from a series of data sets y i where y i = G i ( u ) + η i . (1) u ∈ R k consists of k unknown/uncertain parameters, y i ∈ R m i consists of m i observations spanning the i th data space, η i ∈ R m i is the noise represented as independent zero-mean Gaussian noise with covariance Γ i ;and G i is a nonlinear function (referred to as the forward model) that maps the parameter space to thei th data set. In this paper, we work with two data sets, the dispersion curve that depicts discrete phasevelocity values of surface waves as a function of frequency, and the discrete acceleration time series recordedby the array instruments at different depths. If we combine these two data sets, the inverse problem of (1)is modified as follows: (cid:20) y y (cid:21) = (cid:20) G ( u ) G ( u ) (cid:21) + (cid:20) η η (cid:21) → y = G ( u ) + η, Γ = (cid:20) Γ
00 Γ (cid:21) . (2)We define the covariance matrix of the Gaussian noise as follows:Γ = [ β diag(max | y | )] , Γ = [ β diag( y )] (3)where β and β determine the noise levels for y and y which are acceleration time series and dispersiondata, respectively.For the two data sets relevant to the problem in hand, the forward models are described below:i) For the theoretical dispersion curve we use the transfer matrix approach originally developed by Thom-son (1950) and Haskell (1953) and later modified by Dunkin (1965) and Knopoff (1964). This approachrequires the solution of an eigenvalue problem, for which we use the well known software GEOPSY(Wathelet, 2005). We should mention here that all the examples considered in this study are normallydispersive, so results presented here only require the dispersion curve for the first mode of Rayleighwaves. However, the framework is general enough to allow for higher modes to be incorporated in theinversion process.ii) For the theoretical acceleration time series, we consider wave propagation in a horizontally stratifiedlayered soil of total thickness H and shear wave velocity V s ( z ) varying with depth z . Given an acceler-ation time series at z = H (i.e., the borehole sensor depth), we compute the soil response numericallyusing a finite element model and we use the extended Rayleigh damping (Phillips and Hashash, 2009)to capture the nearly frequency independent viscous damping ξ in time domain analyses. We alsoassume that Vs below z = H is constant, which corresponds to an elastic bedrock idealization. To solve the inverse problem involving the two data sets just described, alone or in conjunction, we use a se-quential data assimilation method (Evensen, 2009) based on ensemble Kalman inversion (Iglesias et al., 2013),a methodology pioneered in the oil reservoir community (Chen and Oliver, 2012; Emerick and Reynolds,2013). In this algorithm, we first define an initial ensemble consisting of N particles. In its most basic form,the ensemble Kalman inversion can regularize ill-posed inverse problems through the subspace propertywhere the solution found is in the linear span of the initial ensemble employed (Chada et al., 2019). Then,at each iteration j , we use the forward model predictions G ( u ( n ) j ) and the observation data y j +1 to updatethese particles sequentially: u ( n ) j +1 = u ( n ) j + C uwj +1 ( C wwj +1 + Γ) − ( y ( n ) j +1 − G ( u ( n ) j )) for n = 1 , . . . , N , (4) Inverse Problem where y ( n ) j +1 can be either identical to y j +1 (the observation data) or found by adding to y j +1 identicaland independently distributed zero-mean Gaussian noise η ( n ) j +1 with distribution the same as that of η ; thematrices C uwj +1 and C wwj +1 are empirical covariance matrices that can be computed at each iteration based onpredictions and the ensemble mean ¯ u j +1 using the following equations. C uwj +1 = 1 N N (cid:88) n =1 ( u ( n ) j − ¯ u j +1 ) ⊗ ( G ( u ( n ) j ) − ¯ G j ) (5) C wwj +1 = 1 N N (cid:88) n =1 ( G ( u ( n ) j ) − ¯ G j ) ⊗ ( G ( u ( n ) j ) − ¯ G j ) (6)and ¯ u j +1 = 1 N N (cid:88) n =1 u ( n ) j , ¯ G j = 1 N N (cid:88) n =1 G ( u nj ) . (7)To enforce physical and prior knowledge systematically, Albers et al. (2019) provide an efficient procedure toimpose constraints within the ensemble Kalman filtering framework; they use the solution of a constrainedquadratic programming to update particles that violate the enforced constraints; in the absence of constraints,the optimization delivers the update formulae above. In the case of linear equality and inequality constraintsthe approach is readily implemented using standard optimization algorithms. We briefly summarize theprocedure in the case of linear inequality constraints: Au ≤ a (8)and when the problem is formulated in the range of the covariance matrix C j +1 : C j +1 = (cid:20) C uuj +1 C uwj +1 C uwj +1 T C wwj +1 (cid:21) (9)where C uuj +1 = 1 N N (cid:88) n =1 ( u ( n ) j − ¯ u j +1 ) ⊗ ( u ( n ) j − ¯ u j +1 ) . (10)It should be noted that the empirically computed covariance is the sum of rank one matrices and its rank isat most N −
1. For this problem setting, N − k + m + m and therefore the covariancematrix is not invertible. To overcome this issue, Albers et al. (2019) reformulate the problem in the rangeof the covariance in which they seek the solution as a linear combination of a given set of vectors. For moredetails see Albers et al. (2019). As such, for each violating particle we seek a vector b ( n ) that minimizes thecost function J j,n ( b ) defined as J j,n ( b ) := 12 | y ( n ) j +1 − G ( u ( n ) j ) − N N (cid:88) m =1 b m ( G ( u ( m ) j ) − ¯ G j ) | + 12 N N (cid:88) m =1 ( b m ) (11)and subject to ABb ≤ a − Au ( n ) j (12)where Bb = 1 N N (cid:88) m =1 b m ( u ( m ) j − ¯ u j +1 ) . (13)Next, we use the computed b ( n ) to update each violating particle u ( n ) j as follows: u ( n ) j +1 = u ( n ) j + 1 N N (cid:88) m =1 b ( n ) m ( u ( m ) j − ¯ u j +1 ) . (14)We have implemented this algorithm – summarized in Algorithm 1 – and have verified its accuracy in thenumerical results section of Albers et al. (2019). Site Characterization: Synthetic Data Algorithm 1
Constrained ensemble Kalman inversion algorithm formulated in range of covariance Choose { u ( n )0 } Nn =1 , j = 0 Calculate forward model application { G ( u ( n ) j ) } Nn =1 Update { u ( n ) j +1 } Nn =1 from (4) for n = 1 : N do if u ( n ) j +1 violates constraints in (12) then b ( n ) ← argmin of (11) subject to (12) Update { u ( n ) j +1 } from (14) end if end for j ← j + 1, go to 2. We first use synthetic data to evaluate the importance of joint inversion of downhole array and dispersiondata in near-surface site characterization using the proposed framework. In this numerical experiment,we use the site conditions at the GVDA in Southern California, which is also the site where we test thealgorithm using recorded ground motion data in the following sections of this work. GVDA is located in anarrow valley, and the near-surface structure consists of an ancestral lake bed with soft alluvium down to18-25 m depth, overlaying a layer of weathered granite; the competent granitic bedrock interface is locatedat 87-m depth according to Bonilla et al. (2002) and varies across the valley (Teague et al., 2018). Figure 1shows the site geology and the layout of sensors at different depths including accelerometers (red boxes) andpressure transducers (blue boxes). Several invasive and non-invasive Vs measurements have been carriedout at this site in the past; the most recent surface wave measurements were performed in October 2016 byTeague et al. (2018), who developed Vs profiles at each of the three surface accelerometer locations (i.e., atthe location of sensors 00, 12, and 21 in Figure 1).To generate ground truth downhole array data, we use acceleration time series recorded at GVDA from anevent with magnitude 4.28. We also consider the Vs profile shown in Table 1 as the target profile; and weassume mass density ρ = 1800 kg/m , Poisson’s ratio ν = 0 . ξ = 0 .
04 for all layers. Wenext use these input parameters in the forward models described above to compute the acceleration time-history at the ground surface and the first mode Rayleigh wave dispersion curve; and use these simulateddata as observations in the inversion process.
Tab. 1:
Vs profile at Garner Valley Downhole Array site used for synthetic data experiments; the values upto depth 100 m are from Gibbs (1989).Layer Thickness [m] Shear Wave Velocity (Vs) [m/s]1 0.0–18.0 2202 18.0–64.5 5803 64.5–150 13004 Half-space 2600When reliable invasive measurements are not available at the site, prior information on the thickness of thesoil layers is also unavailable. To overcome this shortcoming, we use a fine discretization of the profile shownin (15), with increasing thickness increments ∆ h with depth, ranging from ∆ h = 5 m to ∆ h = 25 m. Thisselection resulted in r = 15 layers for the H = 150 m thick profile (from the surface to the depth of boreholesensor 05 shown in Figure 1). Our forward model also considers elastic bedrock boundary conditions beyonddepth 150 m, which we characterize by a thin layer of thickness 1 m in (15).∆ h [m] = { , , , , , , , , , , , , , , } (15) Site Characterization: Synthetic Data Fig. 1:
Garner Valley downhole array cross section (cf. § Data & Resources).Assuming that the Vs is constant in each layer, and the profile is normally dispersive, we enforce monotonicbehavior by posing the following linear inequality constraints, and lower and upper bounds for the very firstand last layers. We assume that the Vs profile can change monotonically between 50 m/s at surface and5000 m/s at the bedrock, which is wide enough search space for our inversion. As we mentioned before,enforcing such constraints reduces the velocity model complexity and the inverse problem ill-posedness. V s,i ≤ V s s,i +1 for i = 1 , . . . , r − , V s, ≥ , V s,r ≤ . (16)For estimation of the small-strain damping, we consider the range 0 . ≤ ξ ≤ .
1, while the density andPoisson’s ratio for each layer are assumed constant and equal to the values we use in our forward modelsimulation: ρ = 1800 kg/m , Poisson’s ratio ν = 0 . P ( V s,i ) = (cid:114) z i H [5 + 10 U (0 , , P ( ξ ) = 0 . . U (0 ,
1) (17)where U (0 ,
1) is the uniform distribution between 0 and 1 and z i is the depth of the bottom of layer i . Thenwe use the projection in (18) to enforce constraints on each particle u ( n ) that violates the constraints in theinitial ensemble. u ( n )0 = argmin u | u − u ( n ) | subject to Au ≤ a (18)where u ( n ) is before we enforce the constraints, and u ( n )0 is its projected counterpart that satisfies the con-straints. Figure 2a shows N = 50 realizations of the projected Vs profiles drawn from the above distributionalong with the ensemble mean and the target profile. Figure 2b shows the corresponding dispersion curvesfor each particle, the mean and the target profile. In the next subsections, we use the ensemble Kalmaninversion algorithm explained in § Inverse Problem to estimate the Vs and/or damping profiles of a horizon-tally layered soil. We specifically consider three test cases to generate the data space y for the inversion.These include: only dispersion data (both a complete and an incomplete set), only downhole array data, andfusion of downhole and dispersion data. We only estimate damping for cases that involve downhole arraydata. Site Characterization: Synthetic Data (a) Shear Velocity [m/s] -150-100-500 E l e v a t i on [ m ] TargetInitial (Mean)Initial (Ensemble) (b) Frequency [Hz] P ha s e V e l o c i t y [ m / s ] TargetInitial (Mean)Initial (Ensemble)
Fig. 2: (a) Vs profiles of the initial ensemble along with the ensemble mean and the target profile – thevertical scale is elevation relative to the ground surface; (b) Dispersion curves computed using theinitial ensemble, ensemble mean an target Vs profiles.
Here, we assume that both active and passive surface wave testing results are available; the former is usuallyused to resolve dispersion data at higher frequencies (shallow layers) and the latter to constrain dispersiondata at long periods (deep layers). The size of the testing arrays and lateral variability in the depth of thebedrock across these arrays may cause difficulties in resolving the dispersion curve continuously for a widerange of frequencies. To examine how data completeness affects the uncertainty of the inverted profile atstrong motion arrays, we define two data sets: a complete and incomplete dispersion curve. In the incompletedata set, dispersion data corresponding to frequencies smaller than 0.3 Hz, and frequencies between 0.45 and2.35 Hz, are missing.We next use Algorithm 1 to iteratively update the particles. Figure 3a shows the estimated Vs profiles, whichare computed from the ensemble mean in the last iteration, using the complete and incomplete data sets, andFigure 3b shows the computed dispersion curves compared against the data sets used as observations. Weshould mention here that the stopping criterion for all performed inverse analyses is reaching 100 iterations;this number is more than enough to make the mean of the ensemble stabilized around the reported finalestimates. The solution non-uniqueness is evident in Figure 3: the dispersion curves of both profiles arein excellent agreement with the curve associated with the target profile, whereas the two inverted profilesshow significant differences. Unsurprisingly, for the case of simulated dispersion data, slight deviation is onlyobserved for the dispersion data of the incomplete set, exactly in the frequency range where information ismissing.
In this section, we repeat the numerical experiment described above (i.e. estimation of Vs and damping)using only downhole array ground motion recordings. The forward problem comprises of propagating the“recorded” (known) motion at depth z = H to the ground surface, and minimizing the misfit betweenthe ground surface motion forward predictions and the surface acceleration time series. Downhole arraysare usually instrumented sparsely (see Figure 1 for example), which may contribute to the solution non-uniqueness and uncertainty. To reflect this issue, we only use the instrument at 150 m depth as input andthe ground surface motion as output, which is the most common configuration of a downhole array (e.g. ofthe Japanese strong motion network, KiK-Net). Site Characterization: Synthetic Data (a) (b) Frequency [Hz] P ha s e V e l o c i t y [ m / s ] CompleteIncompleteEstimated (Mean) - IncompleteEstimated (Mean) - Complete
Fig. 3: (a) Final estimated ensembles along with the ensembles’ mean and target Vs profile – the verticalscale is elevation relative to the ground surface; (b) Dispersion curves computed using the ensemblemean along with the complete and incomplete datasets used in the inverse analysis.To compare with the dispersion data inversion, we first drew the initial particles from the same distributionas in (17) but the algorithm was not successful in finding the profile that can reproduce the output data.However, when we slightly shifted the initial ensemble to Vs values closer to the target solution (see (19)and Figure 4a), the algorithm successfully converged to a Vs profile and damping ratio that captures theground surface acceleration. We should mention here that in this case, where the downhole recorded motionis used as prescribed boundary condition, our forward model is appropriately adjusted to a layered soil onrigid bedrock, in which the thickness of the last layer is 25 m. P ( V s,i ) = (cid:114) z i H [4 + 10 U (0 , , P ( ξ ) = 0 . . U (0 , . (19)Figure 4b shows the estimated Vs profile after 100 iterations and Figure 4c shows the computed surfaceaccelerations using this profile. The final estimate for damping is ξ = 0 . So far we have tried using the two heterogeneous data sets separately: first, we used the dispersion data thatprovided constraints on the dispersive characteristics of surface waves; and then we used downhole arraydata, which provided constraints on the travel time of body (shear) waves through the soil layers. As shownin the last two examples, the inversion algorithm could not recover the target profile in either case whereasall forward models captured the target observations exceptionally well.In this section, we study the effects of fusing these two complementary data sets. We create the initialensemble from (17), and consider y as the combination of the incomplete dispersion data and surface accel-eration time series. We use β = 0 .
01 and β = 0 .
01 to define the Gaussian noise covariance matrix. Usingthe same inversion algorithm, Figure 5 shows the estimated Vs profile and the computed dispersion results.As can be readily seen, by combining the two data sets, we are able to recover the Vs profile with depth;as expected, the profile matches both the incomplete dispersion curve and the acceleration time series onground surface. The latter is almost identical to the time series shown in Figure 4c and therefore is notrepeated here. Furthermore, the algorithm recovers a constant damping ratio of ξ = 0 . Site Characterization: Synthetic Data (a) Shear Velocity [m/s] -150-100-500 E l e v a t i on [ m ] TargetInitial (Mean)Initial (Ensemble) (b)
Shear Velocity [m/s] -150-100-500 E l e v a t i on [ m ] TargetEstimated (Mean)Estimated (Ensemble) (c)
Time [s] -505 S u r f a c e A cc . [ g ] -3 TargetEstimated (Mean)
Fig. 4: (a) Initial ensemble of particles along with the ensemble mean and target profile – the verticalscale is elevation relative to the ground surface; (b) Estimated ensemble of particles along with theensemble mean and target profile; (c) Computed surface acceleration time series using the target Vsprofile and damping as well the mean of the estimated ensemble.
Site Characterization: Synthetic Data synthetic example uses a constant damping ξ = 0 .
04 for all layers. This example shows the effectiveness ofjoint inversion to successfully recover the target Vs profile and damping ratio among possible solutions thatcan fit the observations well individually in the inverse problem setting. In the next subsection, we test thealgorithm effectiveness for more complex profiles and noise-contaminated synthetic data before proceedingto an example using measured dispersion data and recorded downhole array ground motions. (a)
Shear Velocity [m/s] -150-100-500 E l e v a t i on [ m ] TargetEstimated (Mean)Estimated (Ensemble) (b) Frequency [Hz] P ha s e v e l o c i t y [ m / s ] CompleteIncompleteEstimated (Mean)
Fig. 5: (a) Final estimated ensemble along with the ensemble mean and target profile – the vertical scaleis elevation relative to the ground surface; (b) Dispersion curve using the estimated ensemble meanalong with both complete and incomplete datasets used in the inverse analysis.
As shown in the previous subsection, joint inversion of downhole array and dispersion data helps improvethe estimation of the Vs profile. So far, however, we have only used noise-free synthetic data. Prior tointroducing field recorded data, we test the performance of the framework using more complex profiles andnoise contaminated synthetic data. It should be noted that the complex profile we will use in the nextexample is not intended to be representative of the Vs profile at GVDA; rather, we generate it to assess theframework’s capability in dealing with more complex cases.
Here, we use the same framework to estimate the Vs profile for a more complex case. To capture the complexVs, we use a more refined discretization of the profile, i.e., ∆ h = 5 m, namely 30 unknown Vs parameters.We again consider four types of data sets: • complete dispersion data (Case 1); • incomplete dispersion data where phase velocity values for frequencies f ≤ . . ≤ f ≤ . • downhole array data only at the surface z = 0 (Case 3); • fusing data sets in Case 2 and Case 3 (Case 4).For Case 1 and Case 2 we estimate 31 parameters including the elastic bedrock Vs and the Vs from zero to H . For Case 3, we estimate 31 parameters including the Vs from zero to H and damping, and for Case 4,we estimate 32 parameters including the Vs from zero to H , elastic bedrock Vs, and damping. Figure 6aand 6b show the estimated Vs profiles after 100 iterations along with the resulting dispersion curves. On the Site Characterization: Synthetic Data other hand, Figure 6c shows the smoothed relative error of the surface acceleration computed using differentestimated profiles. It should be noted that the absolute error for Cases 3 and 4 is very small. Again, wenotice that joint inversion of dispersion and downhole array data can significantly improve the estimatedprofile. Furthermore, increasing the number of soil layers does not affect the performance of the algorithmpresented here and it can correctly resolve the impedance contrasts in the target Vs profile. (a) Shear Velocity [m/s] -150-100-500 E l e v a t i on [ m ] TargetCase 1Case 2Case 3Case 4 (b) Frequency [Hz] P ha s e V e l o c i t y [ m / s ] CompleteIncompleteCase 1Case 2Case 3Case 4 (c)
Time [s] -2 R e l a t i v e E rr o r [ % ] Case 1Case 2Case 3Case 4
Fig. 6: (a) Finale estimate of the ensembles’ mean along with the target Vs profile – the vertical scale iselevation relative to the ground surface; (b) Dispersion curves computed using the ensembles’ meanalong with the complete and incomplete datasets used in the inverse analysis; (c) Smoothed relativeerror between the surface acceleration dataset and those computed using the estimated ensembles’mean.
To assess the framework robustness when the recorded data is noisy, we next consider Case 4 where weuse the combined data sets to estimate Vs profile and damping. For generating the noise contaminateddata sets y and y , we add a zero mean Gaussian noise where we use (3) to define the covariance matrixconsidering β = 0 .
05 and β = 0 .
05. Then, within the ensemble Kalman inversion iterations, we consider β = 0 .
025 and β = 0 .
025 to define Γ. Figure 7a shows the estimated Vs profile, and Figures 7b and 7cshow the theoretically computed dispersion curve and surface acceleration time series compared to the data(both clean and noisy); this exercise again shows the successful recovering of the Vs profile and damping inpresence of the noisy data.
Site Characterization: Synthetic Data (a) Shear Velocity [m/s] -150-100-500 E l e v a t i on [ m ] TargetEstimated (Mean)Estimated (Ensemble) (b) Frequency [Hz] P ha s e v e l o c i t y [ m / s ] Noise FreeNoisyEstimated (Mean) (c)
Time [s] -6-4-2024 S u r f a c e A cc . [ g ] -3 Noise FreeNoisyEstimated (Mean)
Fig. 7: (a) Final estimated ensemble along with the ensemble man and target Vs profile – the verticalscale is elevation relative to the ground surface; (b) Dispersion curve computed using the estimatedensemble mean along with the noise free and incomplete noisy dispersion data used in the inverseanalysis; (c) Surface acceleration time series computed using the estimated ensemble mean alongwith the noise free and noise contaminated acceleration data used in the inverse analysis.
Site Characterization: Real Data In this section, we use the same inversion algorithm to estimate the Vs profile and damping at GVDA usingfield dispersion data and recorded acceleration time series. For the downhole array data, we use a curateddataset (cf. § Data & Resources), comprising 30 events recorded from 2006 to 2016 with ground surface peakground acceleration (PGA) greater than 10 gal. The events are recorded by vertical and aligned horizontalaccelerometers from 501 m depth to the surface. From these records, we here focus on those with magnitudeMw ≤ (a) -40 -20 0 20 40 W-E (km) -40-2002040 S - N ( k m ) GVDA
Mw = 2.5Mw = 3.5Mw = 4.5 (b)
Distance [km] P G A [ g ] (c) Distance [km] P G V [ c m / s ] Fig. 8: (a) Epicentral distance, (b) PGA and (c) PGV of the events considered for Vs estimation at GVDA.As dispersion data set, we use the mean field measurements associated with the north accelerometer location00 by Teague et al. (2018), who performed both active-source multi-channel analysis of surface waves andpassive source microtremor array measurement testing at GVDA. Using the downhole array data a e =( a e , a e , a e , a e , a e ) for each event e and dispersion data V r , the observation can be formed as y e = ( a e , V r )for e = 1 , . . . ,
23. For each event, we use the inversion algorithm to estimate both V s,e and damping ξ e .Please note that a , a , a , a , a , a are accelerations recorded at z = 0 , , , ,
50 and 150 metersrespectively (see Figure 1). All acceleration time series were filtered using second order bandpass Butterworthfilter with frequency range [0.1, 10] Hz.To start the inversion, we consider three initial ensembles with 50 particles each (see Figures 9a, 9b, and9c); and we seek to estimate the Vs and damping for 23 × e and each initial ensemble i along with the average Vs profile V avg ,is . V avg ,is = 123 (cid:88) e =1 V is,e , i = 1 , , . (20)With the exception of small discrepancies at depths larger than 100 m, the profiles recovered by consideringeach prior and the corresponding ensemble are almost identical. Figure 10 shows the box plot of the estimatedVs at each layer, one for each of the prior distributions that we used to select our particle ensembles. Again,the trend in recovered Vs and associated uncertainty is nearly identical for all three prior Vs distributions.The resulting damping ratios (i.e., the ensemble mean for each e and i ) are shown in Figure 11 as a functionof the event magnitude Mw, PGA and PGV; the average damping is ξ avg = 0 .
049 over 69 events and itappears that higher damping estimation is associated with higher ground motion intensities, which is verylikely the manifestation of nonlinear response. To evaluate the effectiveness of joint inversion, we next usethe second initial ensemble ( i = 2) to estimate the Vs profile using only the dispersion data and only the Site Characterization: Real Data (a)
500 1000 1500 2000 2500 3000 3500 4000
Vs [m/s] -150-100-500 E l e v a t i on [ m ] Initial (Mean)Initial (Ensemble) (b)
500 1000 1500 2000 2500 3000 3500 4000
Vs [m/s] -150-100-500 E l e v a t i on [ m ] Initial (Mean)Initial (Ensemble) (c)
500 1000 1500 2000 2500 3000 3500 4000
Vs [m/s] -150-100-500 E l e v a t i on [ m ] Initial (Mean)Initial (Ensemble) (d)
500 1000 1500 2000 2500 3000 3500 4000
Vs [m/s] -150-100-500 E l e v a t i on [ m ] V savg,1 V s,e1 (Ensemble Mean) (e)
500 1000 1500 2000 2500 3000 3500 4000
Vs [m/s] -150-100-500 E l e v a t i on [ m ] V savg,2 V s,e2 (Ensemble Mean) (f)
500 1000 1500 2000 2500 3000 3500 4000
Vs [m/s] -150-100-500 E l e v a t i on [ m ] V savg,3 V s,e3 (Ensemble Mean) Fig. 9:
Initial ensemble of particles along with the ensemble mean: (a) i = 1, (b) i = 2, (c) i = 3; Finalestimated ensemble mean for events e = 1 , . . . ,
23 along with the average Vs profile: (d) i = 1, (e) i = 2, (f) i = 3. The vertical scale is elevation relative to the ground surface.downhole array data. Figure 12a shows the estimated Vs profiles V dsp s and V dh s for dispersion and downholearray data, respectively; V avg s is the average over all 69 events that we used to search the fused data space.As shown, V dh s closely follows the V avg s while V dsp s starts deviating after z around 30 m. Figures 12b and 12cshow the effects of these inverted Vs profiles on the site signatures, i.e., dispersion curves and site transferfunctions (TFs). The empirical TF for each event is obtained by dividing the Fourier transform of therecorded acceleration signal on the surface a e, † by that of sensor a e, † at depth z = H = 150 m. Figure 12cshows the median TF and its standard deviation among the considered 23 events. It is interesting to notethat:1. although V dh s is following the same trend as V avg s , the resulting dispersion curve cannot capture theexperimental data while they capture the joint inversion TF quite well.2. while V dsp s captures the dispersion data, the associated TF is different from those cases that downholearray data is incorporated in the inversion.3. Joint inversion of the dispersion and downhole acceleration data changes the inverted dispersion dataonly in the frequency range where phase velocity data are not available.4. While peaks in the theoretical TFs are well-aligned with the empirical TF, it is worth noting thatdeviations at higher modes are possibly due to the three dimensional effects.To assess the effects of the inverse problem parameterization on inverted Vs profiles, we consider two morecases of discretizing the soil height with ∆ h = 10 m and ∆ h = 20 m and repeating the joint inversion. Fig-ure 13 shows the estimated Vs profiles and associated dispersion curves and transfer functions. These resultssuggest that decreasing the layer thickness, i.e., increasing the number of Vs values to be estimated, doesnot necessarily have adverse effects on the well-posedness of the inverse problem in hand. Furthermore, ourprevious numerical experiments suggest that fine discretization is successful to capture impedance contrastsat least when dealing with synthetic data. That being said, it is worth noting that, in many cases, it may Site Characterization: Real Data (a)
1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31
Layer Number [-] S hea r V e l o c i t y [ m / s ] (b)
1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31
Layer Number [-] S hea r V e l o c i t y [ m / s ] (c)
1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31
Layer Number [-] S hea r V e l o c i t y [ m / s ] Fig. 10:
Box plot of the estimated Vs at each layer for the considered a priori distributions: (a) i = 1, (b) i = 2, (c) i = 3. The central mark indicates the median, and the bottom and top edges of thebox indicate the 25th and 75th percentiles, respectively. The whiskers extend to the most extremedata points not considered outliers, and the outliers are plotted individually using the ‘+’ symbol. (a) Magnitude Mw D a m p i ng [ % ] (b) PGA [g] D a m p i ng [ % ] (c) PGV [cm/s] D a m p i ng [ % ] Fig. 11:
Final estimated ensemble mean of damping ratio for events e = 1 , . . . ,
23 and as function of (a)event magnitude, (b) event PGA, (c) event PGV. (a)
500 1000 1500 2000 2500 3000 3500 4000
Vs [m/s] -160-140-120-100-80-60-40-200 E l e v a t i on [ m ] V sdsp V savg V sdh (b) Frequency [Hz] P ha s e v e l o c i t y [ m / s ] V s,e ExperimentalV sdsp V savg V sdh (c) Fig. 12: (a) Estimated average Vs profile – the vertical scale is elevation relative to the ground surface; (b)Dispersion curves computed using different average Vs profiles compared against the experimentallyavailable dispersion data used in the inverse analysis; (c) Theoretical transfer functions computedusing different average Vs profiles compared against the median of empirical transfer functions andtheir standard deviation.
Conclusion not be clear which parameterization produces a “better” answer than others. Therefore, to fully understandthe uncertainties, it may be necessary to consider multiple parameterizations. Also, it is worth noting thatour future research is aimed at addressing this issue systematically by estimating both layer thicknesses andVs values. (a) Vs [m/s] -150-100-500 E l e v a t i on [ m ] V savg - h = 5mV savg - h = 10mV savg - h = 20m (b) Frequency [Hz] P ha s e V e l o c i t y [ m / s ] V savg - h = 5mV savg - h = 10mV savg - h = 20mExperimental (c) Fig. 13: (a) Estimated average Vs profile – the vertical scale is elevation relative to the ground surface; (b)Dispersion curves computed using different average Vs profiles compared against the experimentallyavailable dispersion data used in the inverse analysis; (c) Theoretical transfer functions computedusing different average Vs profiles compared against the median of empirical transfer functions andtheir standard deviation.Site characterization of Garner Valley has also been the subject of several geotechnical and geophysicalstudies, and therefore we compare our results to those available in open literature. More specifically, shallowand deep PS suspension logging results were provided by Stellar (1996). Gibbs (1989) used a three componentgeophone in a 100 m deep borehole, and determined the P and S wave velocities using the conventionalmethods of travel-time plots and straight line segments. On the other hand, Bonilla et al. (2002) used atrial-and-error approach to find the best Vs profile that minimizes the misfit between the synthetic and realdownhole accelerometer time series for weak motions; they used the velocity models by Gibbs (1989); Peckerand Mohammadioum (1993) as the initial guess. Finally, Teague et al. (2018) used the active and passivesurface wave measurements to invert for the Vs profile using the dispersion data. They also used the HVSRcurves to further constrain the inversion results by comparing the first mode frequency of theoretical TFsobtained from the inverted Vs profiles against the mean first mode frequency obtained from the experimentalHVSR curves. These results along with those computed in this study are shown in Figure 14 in terms ofVs profile, theoretical dispersion curves, and TFs. In all cases except for Teague et al. (2018) we assumethat ν = 0 .
45 and ρ = 1800 kg/m . For Teague et al. (2018) results, we use the data provided by thefirst author to compute the dispersion curves and TFs considering layering ratio Ξ ranging from 1.5 to 7.Furthermore, for all computed TFs, we assume that the damping ratio is equal to the estimated averagedamping ratio, i.e, ξ avg = 0 . , . , , In this paper, we introduced a sequential data assimilation approach based on ensemble Kalman inversion fornear-surface site characterization. Our method was based on heterogeneous data set fusion and assimilationwith prior knowledge, which we introduced in terms of inequality constraints to the Vs and/or small straindamping. To characterize the general trend of the Vs profile, we used the piece-wise constant function witha known number of layers. Through a series of synthetic experiments, we demonstrated the inverse problemsolution non-uniqueness when dispersion data or acceleration time series were used in isolation and showed
Conclusion (a) Vs [m/s] -150-100-500 E l e v a t i on [ m ] Bonilla et al. (2002)Gibbs (1989)V savg
Teague et al. (2018)Shallow PS (Steller 1996)Deep PS (Steller 1996) (b) Frequency [Hz] P ha s e V e l o c i t y [ m / s ] Bonilla et al. (2002)Gibbs (1989)V savg
Teague et al. (2018)Exp. (Mean STD) (c)
Fig. 14:
Comparison of (a) Vs profiles (the vertical scale is elevation relative to the ground surface), (b)dispersion curves and (c) transfer functions at GVDA.
Data and Resources how the joint inversion of these complementary data could improve the Vs estimation. We also showed thatincreasing the number of layers can help capture more complex profiles without affecting the performance ofthe algorithm. Lastly, we tested the algorithm on real data using the Garner Valley Downhole Array site asour testbed, and compared the inverted Vs profile against previous site characterization studies. Our studyshowed that inversion uncertainties, such as the ones described by Teague et al. (2018), may be attributedto incomplete dispersion data in the medium frequency range. Future non-invasive testing that will helpcomplete the available dispersion data across the entire frequency range of interest will help refine inversealgorithms such as the one presented here. The GVDA downhole array data are available at http://nees.ucsb.edu/curated-datasets and Figure 1 isdownloaded from http://nees.ucsb.edu/facilities/GVDA (last accessed March 2020). The theoretical dis-persion curves are computed using GEOPSY software package gpdc
References
Albers, D. J., Blancquart, P.-A., Levine, M. E., Seylabi, E. E., and Stuart, A. M. (2019). Ensemble kalmanmethods with constraints.
Inverse Problems, http:// iopscience.iop.org/ 10.1088/ 1361-6420/ ab1c09 .Arai, H. and Tokimatsu, K. (2005). S-wave velocity profiling by joint inversion of microtremor disper-sion curve and horizontal-to-vertical (h/v) spectrum.
Bulletin of the Seismological Society of America ,95(5):1766–1778.Assimaki, D., Li, W., and Kalos, A. (2011). A wavelet-based seismogram inversion algorithm for the in situcharacterization of nonlinear soil behavior.
Pure and applied geophysics , 168(10):1669–1691.Assimaki, D., Steidl, J., and Liu, P. C. (2006). Attenuation and velocity structure for site response analysesvia downhole seismogram inversion. pure and applied geophysics , 163(1):81–118.Bonilla, L. F., Steidl, J. H., Gariel, J.-C., and Archuleta, R. J. (2002). Borehole response studies at the garnervalley downhole array, southern california.
Bulletin of the Seismological Society of America , 92(8):3165–3179.Chada, N. K., Stuart, A. M., and Tong, X. T. (2019). Tikhonov regularization within ensemble kalmaninversion. arXiv preprint arXiv:1901.10382 .Chandra, J., Gu´eguen, P., Steidl, J. H., and Bonilla, L. F. (2015). In situ assessment of the g– γ curve forcharacterizing the nonlinear response of soil: Application to the garner valley downhole array and thewildlife liquefaction array. Bulletin of the Seismological Society of America , 105(2A):993–1010.Chen, Y. and Oliver, D. S. (2012). Ensemble randomized maximum likelihood method as an iterativeensemble smoother.
Mathematical Geosciences , 44(1):1–26.Darendeli, M. B. (2001).
Development of a new family of normalized modulus reduction and material dampingcurves . PhD thesis, University of Texas at Austin.Dunkin, J. W. (1965). Computation of modal solutions in layered, elastic media at high frequencies.
Bulletinof the Seismological Society of America , 55(2):335–358.Emerick, A. A. and Reynolds, A. C. (2013). Investigation of the sampling performance of ensemble-basedmethods with a simple reservoir model.
Computational Geosciences , 17(2):325–350.
Data and Resources Evensen, G. (2009).
Data assimilation: the ensemble Kalman filter . Springer Science & Business Media.Foti, S. (2000).
Multistation methods for geotechnical characterization using surface waves . na.Foti, S., Comina, C., Boiero, D., and Socco, L. (2009). Non-uniqueness in surface-wave inversion andconsequences on seismic site response analyses.
Soil Dynamics and Earthquake Engineering , 29(6):982–993.Garofalo, F., Foti, S., Hollender, F., Bard, P., Cornou, C., Cox, B., Dechamp, A., Ohrnberger, M., Perron,V., Sicilia, D., et al. (2016). Interpacific project: Comparison of invasive and non-invasive methods forseismic site characterization. part ii: Inter-comparison between surface-wave and borehole methods.
SoilDynamics and Earthquake Engineering , 82:241–254.Gibbs, J. F. (1989). Near-surface p and s-wave velocities from borehole measurements near lake hemet,california.Haskell, N. A. (1953). The dispersion of surface waves on multilayered media.
Bulletin of the seismologicalSociety of America , 43(1):17–34.Iglesias, M. A., Law, K. J., and Stuart, A. M. (2013). Ensemble kalman methods for inverse problems.
Inverse Problems , 29(4):045001.Knopoff, L. (1964). A matrix method for elastic wave problems.
Bulletin of the Seismological Society ofAmerica , 54(1):431–438.Lunedei, E. and Malischewsky, P. (2015). A review and some new issues on the theory of the h/v technique forambient vibrations. In
Perspectives on European earthquake engineering and seismology , pages 371–394.Springer, Cham.Molnar, S., Cassidy, J., Castellaro, S., Cornou, C., Crow, H., Hunter, J., Matsushima, S., Sanchez-Sesma,F., and Yong, A. (2018). Application of microtremor horizontal-to-vertical spectral ratio (mhvsr) analysisfor site characterization: State of the art.
Surveys in Geophysics , 39(4):613–631.Molnar, S., Dosso, S. E., and Cassidy, J. F. (2010). Bayesian inversion of microtremor array dispersion datain southwestern british columbia.
Geophysical Journal International , 183(2):923–940.Pecker, A. and Mohammadioum, B. (1993). Garner valley: Analyse statistique de 218 enregistrementssismiques. In
Proc. of the 3eme Colloque National AFPS .Phillips, C. and Hashash, Y. M. (2009). Damping formulation for nonlinear 1d site response analyses.
SoilDynamics and Earthquake Engineering , 29(7):1143–1158.Pi˜na-Flores, J., Perton, M., Garc´ıa-Jerez, A., Carmona, E., Luz´on, F., Molina-Villegas, J. C., and S´anchez-Sesma, F. J. (2016). The inversion of spectral ratio h/v in a layered system using the diffuse field assump-tion (dfa).
Geophysical Journal International , page ggw416.Seylabi, E. E. and Asimaki, D. (2018). Bayesian estimation of nonlinear soil model parameters: Theory andfield-scale validation at kik-net downhole array sites. In
Proc. 16th European Conference on EarthquakeEngineering , Thessaloniki, Greece.Socco, L. V. and Boiero, D. (2008). Improved monte carlo inversion of surface wave data.
GeophysicalProspecting , 56(3):357–371.Stellar, R. (1996). New borehole geophysical results at gvda.
NEES@ UCSB Internal Report .Stokoe II, K. (1994). Characterization of geotechnical sites by sasw method, in geophysical characterizationof sites.
ISSMFE Technical Committee .Teague, D. P. and Cox, B. R. (2016). Site response implications associated with using non-unique vsprofiles from surface wave inversion in comparison with other commonly used methods of accounting forvs uncertainty.
Soil Dynamics and Earthquake Engineering , 91:87–103.
Data and Resources Teague, D. P., Cox, B. R., and Rathje, E. M. (2018). Measured vs. predicted site response at the garner valleydownhole array considering shear wave velocity uncertainty from borehole and surface wave methods.
SoilDynamics and Earthquake Engineering , 113:339–355.Thomson, W. T. (1950). Transmission of elastic waves through a stratified solid medium.
Journal of appliedPhysics , 21(2):89–93.Wathelet, M. (2005). Array recordings of ambient vibrations: surface-wave inversion.
PhD Diss., Li´egeUniversity , 161.Wathelet, M., Jongmans, D., and Ohrnberger, M. (2004). Surface-wave inversion using a direct searchalgorithm and its application to ambient vibration measurements.