Virtual Shake Robot: Simulating Dynamics of Precariously Balanced Rocks for Overturning and Large-displacement Processes

Understanding the dynamics of precariously balanced rocks (PBRs) is important for seismic hazard analysis and rockfall prediction. Utilizing a physics engine and robotic tools, we develop a virtual shakerobot(VSR)tosimulatethedynamicsofPBRsduringoverturningandlarge-displacementprocesses. We present the background of physics engines and technical details of the VSR, including software architecture, mechanicalstructure, controlsystem, andimplementationprocedures. Validationexperimentsshowtheme-dianfragilitycontourfromVSRsimulationiswithinthe95%predictionintervalsfrompreviousphysicalexper-iments,whenPGV/PGAisgreaterthan0.08 s . Using a physical mini shake robot, we validate the qualitative consistency of fragility anisotropy between the VSR and physical experiments. By overturning cuboids on flat terrain, the VSR reveals the relationship between fragility and geometric dimensions (e.g., aspect and scaling ratios). The ground motion orientation and lateral pedestal support affect PBR fragility. Large-displacement experiments estimate rock trajectories for different ground motions, which is useful for understanding the fate of toppled PBRs. Ground motions positively correlate with large displacement statistics such as mean trajectorylength,meanlargestvelocity,andmeanterminaldistance. Theoverturningandlargedisplacement processes of PBRs provide complementary methods of ground motion estimation.


Introduction
Precariously balanced rocks (PBRs) are boulders balanced on and not fixed to a sub-horizontal pedestal.The balance configuration and contact physics define PBR fragility-probability for overturning by a stimulus, usually earthquake ground motions.Seismologists have studied the overturning responses of PBRs from ground motions for seismic hazard analysis (Housner, 1963;Brune, 1996;Shi et al., 1996;Anooshehpoor et al., 2004;Rood et al., 2020Rood et al., , 2022)).PBR fragility provides an upper bound on the strength of the ground motions in the time interval since the PBRs became frag-ile (Brune et al., 2006;Anderson et al., 2014).In southern California, most PBRs have been fragile for thousands of years or longer (Brune et al., 2006;Rood et al., 2022).Studying PBRs allows ground motion estimation with long return times, which are much longer than the modern instrumental earthquake catalogs.Such longhistory ground motion estimation is important for assessing hazards for critical facilities such as large dams, nuclear power plants, and nuclear waste repositories (Rood et al., 2020).In principle, PBRs allow a partial test of hazard curves obtained from other information, including geological appraisals of earthquakes and nearby faults (Rood et al., 2020).Hazard curves, outputs of probabilistic seismic hazard analysis (PSHA), ex-press the rate at which ground motions are equaled or exceeded as a function of the amplitude of the motion.PBRs constrain the hazard curves at very long return times.
Seismic hazard analysis typically considers the PBR overturning responses, which are immediate binary results (balanced or overturned, Anderson et al., 2014).However, their motions after overturning are complex and informative.Overturned PBRs can slide, rotate, rock, and bounce.The large displacements of these rocks contribute to understanding the fate of PBRs and the development of rocky slopes.Additionally, many discovered PBRs present rockfall hazards (Anderson et al., 2014).As a serious natural hazard, rockfall poses a major threat to infrastructure, transportation lines, and people (Dorren, 2003;Leine et al., 2013).Predicting rockfall trajectories in complex terrains is essential for implementing protective measures.
Using PBRs for seismic hazard analysis presents challenges in PBR mapping, PBR dynamics, PBR dating, and hazard modeling (Fig. 1).PBRs are not everywhere.PBR mapping locates and obtains PBR geometries and contact properties.Hundreds of PBRs have been manually located in southern California, and their metadata is archived at Southern California Earthquake Center (SCEC researchers, 2022).However, many of them were discovered near developed roads because of accessibility.The heterogeneous distribution indicates a sampling bias for ground motion estimations.A recent study developed unpiloted aerial vehicles (UAV) and onboard machine learning to search for PBRs (Chen et al., 2024, in prep).The PBR geometry is a critical factor affecting PBR dynamics.PBR geometry is often modeled by minimal contact angles from a 2D side-view photo, where a contact angle is an angle between the gravity vector and the connection from the mass centroid to a rocking point (Haddad et al., 2012;Shi et al., 1996).3D PBR models are reconstructed using terrestrial laser scanning, structure from motion (SfM), or robotic realtime mapping technologies (Veeraraghavan et al., 2016;Rood et al., 2020;Chen et al., 2024, in prep).PBR dynamics focus on the response of PBRs from known ground motions (a forward dynamics problem).PBR forward dynamics are nonlinear and have been studied for over a century (Milne and Omori, 1893;Housner, 1963).The methods of PBR dating include cosmogenic, rock varnish, and quantitative geomorphic models (Bell et al., 1998;Rood et al., 2020).With the forward dynamics and ages of PBRs, hazard modeling integrates such information into a PSHA to test or rectify hazard curves.Despite the importance of all four challenges in PBR usage for hazard analysis, this paper specifically concentrates on modeling PBR dynamics through simulation tool development and experiments.
The dynamics of PBRs for seismic hazard analysis are aimed to model an overturning process, which is a mapping from ground motions to PBR response of overturned or not.Ground motions are characterized by intensity measures (Anderson et al., 2014;Purvance et al., 2008), such as peak ground acceleration (PGA), peak ground velocity (PGV), peak ground displacement, SA(1Hz), and SA(2 Hz), where SA(n Hz) Figure 1 Workflow of using precariously balanced rocks (PBRs) for seismic hazard analysis, PBR fate, rocky slope development, and rockfall prediction.PBR dynamics involve overturning process and large-displacement process.
represents the peak acceleration response of a singledegree-of-freedom oscillator with undamped natural frequency of n Hz and 5% damping to the ground motions (Baker et al., 2021).Purvance et al. (2008) found that PGV/PGA and PGA were the strongest indicators of the overturning potentials among the above intensity measures.Since then, PGV/PGA and PGV have been commonly used as ground motion inputs to study the overturning problem (e.g., Veeraraghavan et al., 2016;Rood et al., 2020).The PBR overturning response can be described as a function of PGV/PGA and PGA, OR = f (P GV /P GA, P GA) (1 where OR is a binary variable with value 1 indicating overturned response and value 0 indicating balanced response.Note that in the real world the PBR dynamics are deterministic, which means, given a ground motion, a PBR response can only be either overturned or balanced.In this case, if we uniformly discretize (P GV /P GA, P GA) ground motion space and assume each ground motion has the same probability, a PBR is more fragile than another when it has more overturned responses.To consider uncertainties in PBR dynamics and to integrate the PBR dynamics model into PSHA, a probabilistic model that indicates the probability of being overturned can be obtained from Monte Carlo simulation and logistic regression (Purvance, 2005).The overturning dynamics of PBRs have been studied for over a century.Assuming sufficiently large friction, previous studies modeled 2D rotation motion of rectangles (vertical plane) (Milne and Omori, 1893;Kirkpatrick, 1927;Housner, 1963;Yim et al., 1980;Purvance et al., 2008).Other studies explored the models and criteria for more complicated motions such as rotating, sliding, bouncing, and rocking in 2D (Ishiyama, 1982;Scalia and Sumbatyan, 1996;Shenton and Harry, 1996;Pompei et al., 1998).Purvance et al. (2012) used the discrete element method (DEM) to analyze the overturning response of 3D PBR models.This method required many repeated experiments to calibrate parameters of stiffness and damping (Purvance et al., 2012;Saifullah andWittich, 2022, 2021), and the simulation for each experiment is computationally expensive.Veer-araghavan et al. (2016) presented an alternative method to analyze the 3D PBR overturning dynamics using a constraint-based model.The constraint-based model formed a linear complementarity problem from contact constraints and solved contact impulses by an iterative numerical algorithm (Chapter 2 in Veeraraghavan, 2015).From the contact impulses, contact forces were computed to update object velocity.Their study (Veeraraghavan et al., 2016) validated the constraint-based model in agreement with the physical shake table experiments from Purvance et al. (2008).The constraintbased model from Veeraraghavan et al. (2016) is similar to the collision response stage in physics engines.However, physics engines directly apply the contact impulse to change object velocity instead of computing intermediate contact forces (see Section 2.1).Modern physics engines are also enhanced with efficient algorithms and Graphic Processing Unit (GPU) hardware accelerations.
Besides the overturning dynamics, we also investigate large-displacement dynamics of PBRs, which are important for PBR fate study, rocky slope development understanding, and rockfall prediction.The goal of the large-displacement dynamics is to predict trajectories of PBRs after being overturned.PBR trajectories are affected by factors including PBR initial state, PBR physics properties, terrain morphology, and terrain physics properties.Given the same configurations for all the other factors, a PBR trajectory is distinguished from its initial state, where T is the trajectory (position and orientation with time), h is a function that maps an initial state to a trajectory, and s is the initial state such as position, orientation, and initial velocity.Compared with the overturning dynamics that have been widely explored in previous PBR dynamics studies, general large-displacement dynamics of rocks were more studied in rockfall hazard applications.Early studies restricted rockfall motions to 2D vertical planes and built mathematical models to describe discrete motion modes (Bozzolo and Pamini, 1986;Kobayashi et al., 1990;Azzoni et al., 1995).3D rockfall models were developed to simulate particle interactions with digital elevation models and digital terrain models (Gascuel et al., 1999;Agliardi and Crosta, 2003;Lan et al., 2007;Guzzetti et al., 2002;Caviezel et al., 2019).Recently, Hao et al. (2021) used a physics engine to simulate rockfall trajectories on a terrain model that was reconstructed by aerial photographs from unpiloted aircraft systems (UAS).By integrating advanced technologies from physics engines and robotics, we have developed a virtual shake robot (VSR) to facilitate the study of the overturning and large-displacement dynamics of PBRs.Our VSR relies on three core technologies: Robot Operating System (ROS), Gazebo simulation toolbox, and Bullet physics engine.ROS is a software platform that provides a set of libraries and tools for robot control and perception (Stanford Artificial Intelligence Laboratory, 2018).Gazebo is a simulation toolbox that provides a simple way to build a virtual world including virtual robots and environments (Koenig and Howard).The dynam-ics of the virtual world are managed by a physics engine.Gazebo supports four high-performance physics engines: Open Dynamics Engine, Bullet, Simbody, and DART.Because the Bullet physics engine has shown reliable simulation results in many scientific studies (Zhu and Zhao, 2019;Ma et al., 2018;Sun et al., 2019), we select Bullet as the physics engine for the VSR.
This study is aimed to advance the simulation of PBR dynamics by seamlessly addressing both overturning and large-displacement processes.Simulation tools are important for science studies, but technical nuances may affect experiment results.For example, rock behaviors such as rotating, sliding, jumping, and rocking affect the dynamics of a shake table.However, such details are unclear in previous studies (Purvance et al., 2012;Veeraraghavan et al., 2016).To reduce the effects of the coupling dynamics, the VSR implements a hierarchical control system.
Previously, simulation models for the overturning and large-displacement dynamics were built independently.Our VSR is the first tool for both dynamics studies.Additionally, complex terrain-anything other than a flat pedestal-can either increase or decrease the fragility of a PBR, depending on the specific characteristics of the surrounding terrain and the contact.Our VSR supports arbitrarily complex terrains, e.g., mesh models from UAS and SfM, to advance PBR dynamics studies.From shake simulations, we demonstrate that surrounding pedestals are critical to reduce PBR fragility.Lastly, we expect to see more robotic applications to PBR mapping, and thus the VSR can be integrated into autonomous mapping systems where ROS has widely been used.
To enhance the clarity of the article structure, we provide an outline as follows.Following the introduction in Section 1, we present technical advancements and review related work in Section 2. Because a physics engine is a new approach to simulating the PBR dynamics, we begin by reviewing the technical details in the Bullet physics engine and compare it with the DEM to establish the necessary context for our research.We also introduce relevant applications of physics engines.Moving on to Section 3, we provide a comprehensive description of the VSR system, covering its software architecture, mechanical structure, control system, and implementation procedures.For validation purposes, Section 4 compares the simulation results from the VSR with physical experiment data from a previous study, as well as with data collected from a physical, mini shake robot.In Section 5, we present our simulation experiments, with a focus on scientific applications related to the overturning and large-displacement processes.In Section 6, we discuss the validation experiments, the overturning and large-displacement experiments, as well as the limitations and prospects for future work.Finally, we summarize this study in Section 7. Throughout this paper, we use the terminology 'robot' to refer to both the VSR and mini shake robot because of their reliance on robotic concepts, including control and perception modules, as well as the use of robotic tools such as ROS.It is worth emphasizing that both robots facilitate the automation of data collection, making them valuable tools for obtaining data in earthquake studies.

Bullet Physics Engine
The Bullet physics engine simulates rigid body dynamics by computing rigid body states (poses and velocities) in a discrete simulation loop with a fixed small time step in the range of 30 Hz to 10 kHz (Coumans and Bai, 2016).
The Bullet rigid body simulation loop includes four stages, as shown in Fig. 2a.The collision detection stage predicts contact points (where to contact) and time of impact (when to contact) using efficient algorithms (van den Bergen, 2003;Williams et al., 2014).When collision is predicted within a simulation time step, the collision response stage computes the impulse from the collision.The forward dynamics stage computes external force, torque, and inertia.Finally, the numerical time integration stage updates position and linear velocity of each object using semi-explicit Euler integration method, where J [kg • m • s −1 ] is the total impulse from the contact, which is computed at the collision response stage.
Similarly to the linear equations, the angular equations update angular velocity and orientation by considering torque, inertia, and angular acceleration.Such a method of computing contact impulse was known as impulse-based dynamics (IBD), introduced by Mirtich and Canny (1995) and improved by Bender et al. (2013).
Bullet computes the contact impulse by modeling collision dynamics as equality and inequality constraints to form a mixed linear complementarity problem (MLCP).
The constraints from a collision include contact constraints, friction constraints, and joint constraints.The projected Gauss-Seidel algorithm (Erleben et al., 2005) solves the MLCP by iteratively approximating an impulse until all the constraints are satisfied (converged) or a maximum number of iterations is met.

Physics Engine versus Discrete Element Method (DEM)
Classical DEM discretizes each object as a collection of small spheres (Cundall and Strack, 1979).For each fixed small time step in discrete simulation loops, as illustrated in Fig. 2b, DEM computes the states of all spherical particles.The rationale of DEM is that the time step chosen is so small that during a single time step disturbances cannot propagate from any spherical particles further than its immediate neighbors.The total force on each spherical particle is only determined by external force defined by the user and contact force imparted by its neighbors with which it is in contact.DEM allows penetration (geometric overlaps) between spherical particles, and thus does not need collision detection to predict when and where to collide.When spherical particles contact, DEM uses a force-displacement model to compute contact force from penetration depth.The processes of the forward dynamics stage and numerical time integration stage in DEM are similar to the processes in Bullet.While originally intended to represent granular media with spherical particles, DEM has since been extended to simulate the behavior of arbitrary rigid bodies (meshes) using polyhedral particles (Cundall, 1988;Itasca Consulting Group, Inc., 2020).However, DEM computes the contact forces at all particles' penetrating faces and nodes, which is different from Bullet where the collision computation is based on each individual mesh object.In DEM, a more highly discretized surface within an object provides a more accurate representation of the contact force distribution.However, increasing the number of discretized particles significantly increases the simulation time.
The difference in the contact models of the Bullet physics engine and DEM is a major factor affecting their computational efficiency.Bullet applies a hard contact model where penetrations are not allowed between any two rigid bodies.Thus, in its simulation loop, Bullet needs the collision detection stage to predict when and where to contact, and then collision impulses are computed in the collision response stage.The hard contact model is formalized as a MLCP such that iterative numerical methods (e.g., projected Gauss-Seidel algorithm) are leveraged to efficiently compute collision impulses.The Bullet physics engine requires users to provide macro physics parameters, including restitution, Coulomb friction, and rolling friction coefficients (Chen, 2022).On the other hand, DEM adopts a soft contact model (force-displacement model) where penetrations are allowed between two directly contacting particles.The soft contact model requires user-defined parameters such as stiffness, damping constant, and friction coefficient.To simulate rigid bodies, the stiffness is usually set at a very large value to reduce macroscopic deformation.Because of the high stiffness, the simulation time step must be small to yield small elastic rebound at each iteration to ensure numerical stability, which significantly increases the simulation time.Additionally, the hard contact model treats each object as an individual entity (e.g., using polygon mesh model), whereas DEM discretizes each object into small particles.The computational time and memory of DEM increase markedly with the number of particles.With the same computational hardware and settings, previous studies found physics engines were 10-250 times faster than DEM to achieve similar results (He et al., 2020(He et al., , 2021)).Given the fact that the soft contact model in DEM is a hypothetical model, the user-defined parameters within this model lack direct connections to macro physics properties.These parameters can be calibrated through an iterative process of adjusting parameter values and matching experimental observations.In contrast, the parameters required in Bullet represent macrophysics properties and may be directly measured through experiments.
Physics engines and DEM are numerical methods, neither of which is a true representation of reality, and both of which need calibration.Our goal here is to promote physics engine applications in PBR dynamics studies, which provides one more option with some advantages relative to DEM.More research is needed to com-pare the performance of physics engines and DEM on this topic.

Physics Engine Applications
In this subsection, we review physics engine applications related to this study.Physics engines have succeeded in simulating behavior of granular assemblies.Izadi and Bezuijen (2014) used Bullet to simulate the behavior of granular materials subjected to pluviation and vibration.Their simulation results were within the range of repeated laboratory experiments.Toson and Khinast (2017) applied Bullet to study quasi-static granular flows of non-spherical particles.While their Bullet simulation results for spherical particles were in agreement with DEM simulation results, implementing non-spherical particle simulation in Bullet was easier because simulations of non-spherical particles in DEM required an advanced discretization method (e.g., Lu et al., 2015).Their Bullet simulation results of nonspherical particles agreed with the prediction from the empirical Beverloo equation (Beverloo et al., 1961).He et al. (2020) compared physics engine and DEM simulations in granular soil behavior and showed that the physics engine achieved similar results to those of the DEM in a significantly shorter time.Komaragiri et al. (2021) demonstrated that the compaction behavior of asphalt mixture from Bullet simulation was very similar to the compaction behavior recorded from laboratory measurements.
Zhu and Zhao (2019) demonstrated the benefits of utilizing physics engines in material analysis.Their study employed the Bullet physics engine and integrated peridynamics to simulate crushable granular materials under mechanical loading.It was difficult to analyze particle breakage using DEM because traditional spherical particles were incapable of approximating particles with sharp corners and edges (Zhu and Zhao, 2019).The Bullet simulation results were consistent with experimental observations on normal compression line, particle size distribution, fractal dimension, and particle morphology.
The Bullet physics engine demonstrated reliable capability of simulating overturning and largedisplacement dynamics (Ma et al., 2018;Sun et al., 2019;Hao et al., 2021).Ma et al. (2018) applied Bullet to simulate rocking dynamics of cuboids with various sizes and aspect ratios (height-to-width ratios).Their study demonstrated that the response of rocking blocks in Bullet simulation was consistent with analytical solutions of Housner equations (Housner, 1963).Sun et al. (2019) used Bullet to analyze overturning dynamics of agricultural tractors on a bank slope and on a uniform slope.The Bullet simulation results were similar to previous reports and also in reasonably good agreement with experimental results.Hao et al. (2021) simulated rockfall trajectories on terrains using PhysX physics engine.They showed the advantage of using the physics engine to simulate rock interactions with a high-resolution, realistic terrain model that was reconstructed by UAV aerial photographs.

Virtual Shake Robot
As shown in Fig. 3, we developed the VSR using ROS, Gazebo simulation toolbox, and Bullet physics engine.Utilizing various libraries and tools available in the ROS ecosystem, we built robot software composed of control and perception modules.The control module computed actuation forces for the VSR; the perception module monitored PBR states and ground motion.Using Gazebo, we designed the mechanical structure of the VSR and defined the general physics properties, including gravity and lighting.Gazebo simulation toolbox supports four physics engine options to manage the dynamics of the virtual world.We selected the Bullet physics engine for its reliable performance in many previous scientific studies (see Section 2.3; Zhu and Zhao, 2019; Ma et al., 2018;Sun et al., 2019).Implementing such a software architecture provides two main advantages.First, rather than directly interacting with the details in a physics engine, Gazebo only requires XML format configuration files where the user can select the desired physics engine and configure physics parameters.Once the configuration files are properly set, Gazebo passes the parameters to the physics engine, simplifying the user's experience and making it easy to switch between physics engines if necessary.Second, as illustrated in Fig. 3, leveraging ROS in the development of the robot software allowed us to reuse the software for both virtual and physical shake robots (Chen et al., 2022).The same robot software ensures that ground motion generation processes are consistent, which is critical to compare simulation and physical experiments.
The VSR has a straightforward mechanical structure, as illustrated in Fig. 4, composed of a base, linear rail, and pedestal.All three components are rigidly fixed, with the base anchored in the virtual world and the linear rail attached to the base.A prismatic joint, functioning as a prismatic motor, links the linear rail and pedestal and generates translational force to actuate the pedestal.This translational force is calculated from the control module.The VSR enables one-dimensional prismatic, horizontal ground motions under the constraint of the prismatic joint (Fig. 4).The VSR supports various pedestal models, including a flat terrain (Fig. 4a) and realistic terrains that were mapped from UAS and SfM (Fig. 4b, c).Switching a pedestal model is as simple as configuring the model path in a Gazebo configuration file.
We developed a hierarchical control module for the VSR (Fig. 5a).The control module provides two ground motion options.The first one is a single-pulse cosine ground displacement, where d(t) is the ground displacement function, A is the amplitude, f is the frequency, and t ∈ [0, 1/f ] is time.A and f are derived from PGV/PGA and PGA, where r is PGV/PGA, a is PGA, and g is the gravitational acceleration.As shown in Fig. 5a, the motion interpreter converts PGV and PGA to f and A using Eq. 8 and Eq. 9. Based on A and f , the trajectory planner takes the derivative of the cosine displacement function (Eq.7) to obtain ground velocity function, where v(t) is the ground velocity function, and t ∈ [0, 1/f ] is time.We uniformly discretized v(t) to sample a set of velocities, {v}, as the input of the velocity controller.The sampling frequency is a user-defined parameter, usually between 100 and 200 Hz.
In addition to the cosine ground displacement function, the VSR supports ground motion simulation from real seismometer records.As shown in Fig. 5a, a lowpass filter first removes the high-frequency noise in the raw acceleration data.The numerical integration produces velocities from the accelerations.Then a highpass filter removes the low-frequency noise in the velocities, because the low-frequency velocity noise may accumulate displacement errors.Note that the output of the high-pass filter is a set of velocities {v}, which has the same format as the output from the trajectory planner.We utilize a shared velocity controller to process the desired velocity commands, simplifying the control software.
We implemented a PID velocity controller to generate force commands for the prismatic joint.The velocity commands for the PID controller are derived from either a cosine ground placement function or a seismometer record.The output force from the PID velocity controller actuates the pedestal.At the same time, Gazebo measures the actual velocity of the pedestal and feeds it back to the PID velocity controller.The actual velocity (velocity measurements) of the pedestal may be different from the desired velocity (velocity commands), resulting in velocity error, where v m is the actual velocity from measurement, and v d is the desired velocity from the higher level of the control module.One reason for the velocity error is that PBR overturning behaviors may produce collision impacts that affect the pedestal dynamics.The objective of the PID velocity controller is to generate force for the prismatic joint to minimize the velocity error, where F is the output force command, and K p , K i , K d are non-negative coefficients for the proportional, integral, and derivative terms, respectively.For example, if the desired velocity is greater than the actual velocity, the velocity error becomes positive.Thereby, the proportional component K p e(t) increases to generate larger force, which increases the actual velocity and reduces the velocity error.The integral and derivative components make the response of the controller more stable and responsive.To further reduce the effects of the coupling dynamics between PBR and pedestal, we set a large mass for the pedestal.A larger mass pedestal with a larger inertia can absorb more energy from the collision, resulting in a smaller collision-caused velocity change.
We developed an automation program to repeat overturning experiments with different single-pulse cosine ground motions.Fig. 5b illustrates the automation workflow.The objective of the automation process is to obtain PBR overturning responses (Eq. 1) to a linear mesh of PGV/PGA and PGA.From the set {(PGV, PGA)}, the automation program popped every pair of (PGV, PGA) and sent it to the controller to actuate a single-pulse cosine ground displacement until the set was empty.To ensure the consistent initial PBR position and orientation across all experiments, the automation program loaded and deleted the PBR model before and after every experiment, respectively.The automation program recorded the PBR states during the With the Bullet physics engine, the VSR supports meshed models of terrains mapped in the real world.For example, Fig. 6 shows the process of creating 3D geometric models of terrain and PBR.The terrain and PBR were mapped by UAS and SfM at a study site of Double Rock, which is located close to both the population center of San Luis Obispo and the critical infrastructure at Diablo Canyon in south-central coastal California (The PBR is named DRE2 in Rood et al., 2020).We first separated the terrain and PBR in the mesh model reconstructed by SfM.The separated terrain, however, lacked supporting and lateral surfaces, which were not reconstructed because of their invisibility without physically removing the PBR.To address this, we completed the missing surfaces by manually adding planes such that the slopes of the added planes were close to the slopes of the local jointing surfaces on the terrain.Similarly, we added planes to the PBR geometric model and built a closed-surface mesh model using Poisson reconstruction, as shown in Fig. 6.From a closed-surface mesh model, we used CAD software (e.g., Autodesk Fusion 360) to measure the mass and moment of inertia of the PBR, which is made of chert with a density of 2,110 kg/m 3 .The height of the PBR is approximately 12 m Geometric simplifications like this, particularly at the base or interface, introduce potentially significant uncertainties in overturning fragility.Complex basal conditions, as is the case for many PBRs, effectively introduce multiple points of rocking or potential uplift.The resulting increase in fragility was evident in the shake table tests of Purvance et al. (2008) and the analytical model of Wittich and Hutchinson (2017).More recent shake table testing by Saifullah and Wittich (2021) quantified that the overturning demand can vary up to ±50% because of small modifications in the basal geometry from surveying techniques.The basal contact can be made arbitrarily complex if necessary in the VSR, but that was not the goal of this paper.Although the geometric modeling approach taken herein does not fully capture the basal interface, this paper aims to provide a demonstration of a first-generation technology for modeling the dynamics of PBRs.

Velocity Controller
We evaluated the performance of the PID velocity controller in the VSR using a single-pulse cosine displacement ground motion and realistic ground motion.The objective of the PID velocity controller was to generate the translation force to actuate the pedestal (flat surface or realistic terrain), such that the actual velocity measured from the pedestal (ground motion velocity) in the simulation closely matched the desired velocity.Fig. 7a presents an example where the desired velocity and actual velocity from a single-pulse cosine displacement ground motion with PGA of 0.2 g and PGV/PGA of 0.8 s.The overlap of the two velocities in Fig. 7a demonstrates the good performance of the PID velocity controller.We observed similar matches between the desired velocity and actual velocity in other single-pulse cosine displacement ground motions.To evaluate the realistic ground motion, we calculated the desired velocity based on raw acceleration data collected in an accelerometer on a physical shake table.Fig. 7b shows the desired and measured velocities from the VSR.During the shake test, the pedestal experienced a rapid disturbance in the actual velocity, as highlighted in the box in Fig. 7b, caused by the overturned PBR.The PID controller was able to quickly correct the velocity disturbance, demonstrating the robustness of the controller.Additionally, we used this realistic ground motion to shake the Double Rock PBR on a flat surface and on the realistic terrain (in the yaw 0°direction).This realistic ground motion overturned the Double Rock PBR on the flat surface but did not on the realistic terrain with the surrounding pedestals, demonstrating the effects of surrounding pedestals on PBR fragility.

Previous Overturning Experiments
We validated the overturning dynamics by comparing the shake experiments from the VSR and Purvance et al. (2008).Referencing the known height of a nearby steel I-beam section, we estimated the dimension of a wooden block (labeled W2) from Figure 4 in Purvance et al. (2008) as 5.5 × 1.1 × 1.1 cm.We modeled a cuboid with the same dimension and a density of 1,500 kg/m 3 (wooden density) in Gazebo.The cuboid was placed on flat terrain with coefficients of Coulomb friction, rolling friction, and restitution of 0.6, 0.6, and 0.2, respectively.The two friction coefficients were selected based on dry rock friction (Byerlee, 1978), and the restitution coefficient of a wood block was measured from a normal drop test (Haron and Ismail, 2012).Using the automation program (Fig. 5b), the VSR conducted 2,500 overturning experiments on a linear mesh space of single-pulse cosine displacement ground motions where PGV/PGA ∈ [0.05, 0.35] and PGA ∈ [0.05, 0.5].For each PGV/PGA, the VSR increased the PGA from 0.05 g to search for the first PGA that overturned the cuboid.Fig. 8 shows the results from the VSR and Purvance et al. (2008).The blue curve delineates the first overturning PGAs of the cuboid from the VSR experiments.The squares represent the results of the wooden block W2 from the physical overturning experiments by Purvance et al. (2008).The gray-filled region, showing prediction intervals about the fragility contours within which 95% of the overturning responses occur, was calculated from empirical equations based on the square data (Purvance et al., 2008).The results from the VSR were within the 95% prediction intervals when PGV/PGA is greater than 0.08 s.When PGV/PGA is smaller than 0.08 s, the 95% prediction intervals were below the VSR curve (Fig. 8), indicating that the cuboid in the VSR was less likely overturned by ground motions with small periods.Note that the VSR curve was resulted from 2,500 overturning experiments with different parameters densely sampled from the ground motion space.However, the 95% prediction intervals were calculated based on a small number of data points (the squares), and only one data point was collected for PGV/PGA smaller than 0.08 s.Future work should collect more data points to examine the 95% prediction intervals, especially on the small PGV/PGA space.

Mini shake robot
We developed a physical, mini shake robot to validate the simulation (Chen et al., 2022), as shown in Fig. 9 The mechanical hardware design of the mini shake robot adopts a closed-loop stepper motor for actuation and a toothed belt for transmission.We developed the software of the mini shake robot based on ROS, which had a similar robotic software architecture as the VSR (Fig. 3).The mini shake robot uses the same high-level and middle-level control programs as the VSR.The primary difference is the low-level velocity controller: the mini shake robot employs a closed-loop stepper motor with toothed belt transmission, whereas the VSR uses a prismatic motor (prismatic joint).The closed-loop stepper motor consists of a regular stepper motor (specifically, NEMA 34HS31) and an encoder that measures the shaft orientation, forming the feedback loop for the lowlevel PID controller.With a toothed pulley of 53.98 mm outside diameter, the stepper motor enables the bed to reach a maximum horizontal acceleration of 1.2 g and a Figure 8 Overturning experiment results of a wooden block (W2 from Purvance et al., 2008).Squares and curve indicate the experimentally observed PGA at which the block overturns for the first time, as PGA is gradually increased from a small value for each PGV/PGA.Gray-filled region delineates prediction intervals about the fragility contours within which 95% of the overturning responses occur.The mini shake robot provides a reverse method for simulation validation.Using the mini shake robot, we conducted experiments with small-scale, free-standing blocks such as 3D-printed PBRs.We then repeated these experiments in simulation using the VSR to compare the physical and simulation results.Specifically, we downscaled the Double Rock PBR from the height of 151.0 cm to 12.8 cm and used polyethylene terephthalate glycol (PETG) material with a density of 1,240 kg/m 3 to 3D print the PBR.The PETG PBR weighted 403.5 g.Following the instructions in Anooshehpoor et al. (2004), we applied grip tapes on the bed to increase friction (no sliding friction during the experiments).Previous studies used a crane machine to lift and rest a heavy PBR after each overturning experiment (e.g., Saifullah and Wittich, 2021).Because of the small size and light weight of the PETG PBR, we were able to reset its pose precisely without the use of a crane.Based on the geometry and physical properties of the PETG PBR, we created a free-standing block with the same dimension in simulation and validated its overturning dynamics using the VSR.Fig. 10 presents the overturning validation results.We used the logistic regression model (Purvance, 2005) to approximate the boundary curves between the overturned and balanced responses.Despite the observation that the real PETG PBR was more fragile, the boundary curves from the simulation and physical experiments were close.In the second experiment, we downscaled the Double Rock PBR to a height of 12.0 cm and 3D-printed it with polylactic acid (PLA) material, which has a density of 1,250 kg/m 3 , resulting of a weight of 105 g.To validate fragility anisotropy, we placed the PLA PBR with two initial orientations (yaw angles of 0°and 270°) on the bed.Using the mini shake robot, we obtained the response diagrams from a set of ground motions, as shown in Fig. 11.The response diagrams show that PLA PBR oriented at a yaw angle of 0°is more fragile than at yaw angle of 270°.This result is consistent with analysis from the previous studies, which found PBRs with smaller minimal contact angle along the motion direction is more fragile (Purvance et al., 2008;Haddad et al., 2012).The resulting boundary curves have a similar pattern to the simulation results of the Double Rock PBR with original dimensions and chert density (see Section 5.2).

Experiments
We conducted a set of experiments to demonstrate the applications of the VSR, test the limitations of certain PBR assumptions, and investigate the PBR overturning and large-displacement processes.
Fig. 12b-e shows the overturning response to singlepulse cosine ground motions.In the cuboid's coordinates, the horizontal ground motions were along the x axis (as the red arrow indicated in Fig. 12a).By comparing Fig. 12b and c, the cuboid with a larger heightwidth ratio was more fragile.With the same heightwidth ratio, the smaller cuboid was more fragile by comparing Fig. 12b and d.These findings are consistent with the previous reports on rectangle experiments on a 2D plane (Purvance et al., 2008;Anderson et al., 2014;Yim et al., 1980).The cuboids with 1×1×3 meters and 1×2×3 meters had similar overturning response diagrams except for a small number of noises (Fig. 12c,  e), which suggested the overturning was not affected by the cuboid materials on the y axis (as the green arrow indicated in Fig. 12a).For an object with more complex geometry, however, extending materials on the y-axis is anticipated to lead to a more intricate overturning response, such as twisting behaviors, which would need further investigation to be confirmed.

PBR Fragility Anisotropy on Flat Terrain
Because geometries of natural PBRs are often asymmetrical, the overturning responses should vary from different ground motion directions-fragility anisotropy.However, most previous studies simplified PBR geometries and only considered the minimal contact angle with the pedestal (e.g., Purvance et al., 2008;Haddad et al., 2012).Veeraraghavan et al. (2016) studied fragility anisotropy based on a rigid-body dynamics algorithm (Chapter 2 in Veeraraghavan, 2015).Using the VSR, we examined the fragility anisotropy of the Double Rock PBR (Fig. 6) simply by placing the PBR on a flat pedestal in 12 different initial orientations (spaced every 30°on yaw).The orientation of PBR was defined in Fig. 13a.
By placing the PBR with different orientations, we simulated the varying ground motion directions.We set the Double Rock PBR with the same physical properties as the cuboids described above and applied the automation program (Fig. 5b).
Fig. 13b-f depict the overturning response diagrams for the fragility anisotropy study on a flat pedestal.From the opposite directions such as Fig. 13b, d or Fig. 13c, e, the same single-pulse cosine ground displacements produced different overturning responses.By comparing Fig. 13b-d, the PBR was less fragile to the 90°ground motions than 0°or 180°ground motions.In Fig. 13f, some balanced responses separate a small cluster of toppled responses on the left (green polygon) and the major boundary curve.For this ground motion direction (yaw 300°), the logistic regression method proposed by Purvance (2005) would not be the ideal model to compute the overturning probability, because the logistic regression model, calculated based on the data near the first PGA for a PGV/PGA, would fail to take into account the balanced responses between the two clusters of toppled responses.

PBR Fragility Anisotropy on Realistic Terrain
Using the VSR, we investigated the effects of surrounding pedestals on the overturning responses of the Double Rock PBR.Most previous studies on the dynamics of free-standing blocks assumed that the blocks had no interaction with adjacent objects.Konstantinidis (2008) explored the overturning dynamics of a free-standing block anchored to a wall via chains.Bao and Konstantinidis (2020) investigated 2D analytical models to study the dynamics of a free-standing block considering impact with an adjacent wall.However, no previous PBR studies have directly accounted for lateral supports in 3D.In this experiment, we set up the terrain model in different orientations to study fragility anisotropy with surrounding pedestals (Fig. 14a, c).The Double Rock PBR had the same physical properties as described earlier.We set the same contact parameters (friction and restitution) for the terrain model.Fig. 14b and d show the overturning response diagrams of the Double Rock PBR on the realistic terrain model.Note that the ground motion ranges in Fig. 14 are larger than those in Fig. 13.By comparing Fig. 13b and Fig. 14b, the surrounding pedestals significantly reduced the PBR fragility, when the ground motion was along the yaw 0°direction.From Fig. 13c and Fig. 14d, the surrounding pedestals only slightly reduced the PBR fragility, when the ground motion was along the 90°direction.The effects of the surrounding pedestals varied for different ground motion directions.

PBR Large Displacement on Realistic Terrain
In this experiment, we investigated the largedisplacement dynamics induced by ground motions.
Using the VSR, we recorded the PBR states during the experiments of PBR fragility anisotropy on the Double Using the recorded trajectories, we analyzed the relationship between ground motions and large displacements.As shown in Fig. 15c-f, the number of trajectories increases with PGA (increase in overturned PBRs).Comparing the trajectory plots (Fig. 15c-f) reveals how trajectory and velocity change with PGAs.Fig. 16 highlights the relationship between trajectory lengths with ground motions.Concurrently large PGA and PGV/PGA result in a long trajectory.Only one large value in either PGA or PGV/PGA is insufficient to produce a long trajectory.Fig. 16b, c show an increasing trend of trajectory length as PGA or PGV/PGA increases.The trajectory data points around 6 m (Fig. 16b, c) correspond to situations where the overturned PBR lands on the niche position shown as the white box in Fig. 15a.
We conducted a correlation analysis to quantify the relationship between ground motions and large displacements.The correlation analysis used the 692 trajectories resulting from the overturned states (red data points in Fig. 14b).Each trajectory is characterized by its trajectory length, largest velocity, and terminal distance.The terminal distance is the Euclidean distance between the start and terminal positions.For each PGA, we ensembled the trajectories from all corresponding PGV/PGAs (all red data points on a PGA column in Fig. 14b).Similarly, for each PGV/PGA, we ensembled the trajectories from all corresponding PGAs.In each ensemble, we calculated the trajectory statistics such as mean trajectory length, mean largest velocity, and mean terminal distance.Fig. 17 plots the trajectory statistics and ground motions.The PGA and PGV/PGV positively correlates with mean trajectory length, mean largest velocity, and mean terminal distance, as the R 2 and p value summarized in Table 1.The correlation analysis results reject a null hypothesis that the trajectory statistics and ground motions are uncorrelated.

Validation
The validation experiments examined the performance of the VSR in terms of the velocity controller and overturning dynamics.The disturbance in the actual velocity in Fig. 7b shows the coupling dynamics between the PBR and pedestal.However, the effects of such coupling dynamics have often been neglected in previous studies.The quick correction in the actual velocity (Fig. 7b) indicates that the PID controller was robust to such disturbance.The overlay between the desired velocity and actual velocity in Fig. 7 shows the good performance of the PID velocity controller.Additionally, the realistic ground motion (0.4 g PGA and 1.2 s PGV/PGA) overturned the Double Rock PBR on the flat surface but did not on the realistic terrain with the surrounding pedestals, which agrees with the overturning response diagrams in Fig. 13b and Fig. 14b.
We compared the overturning results from the VSR and Purvance et al. (2008).When PGV/PGA was greater than 0.08 s, the overturning results from the VSR were consistent with Purvance et al. (2008).The reasons for the difference from the low PGV/PGAs remain to be explored.However, we observed that the first overturning PGAs in the VSR were less scattered from the median fragility contour than Purvance et al. (2008).More physical experiment data points in the low PGV/PGA region are needed to refine the 95% bounds.Additionally, low PGV/PGA motions typically exhibit higher-frequency  movements.On physical free-standing structures, highfrequency movements may trigger a slight rocking or wobbling because of imperfections in the basal contact.Because rocking behavior is nonlinear to orientation, even such a slight rocking from higher-frequency move-ments can increase susceptibility to overturning, especially when compared to VSR predictions that do not account for physical imperfections at the base.Thus, validating the rocking behaviors resulting from highfrequency movements would require more physical ex- We built the mini shake robot to examine the performance of the VSR on PBR overturning dynamics.Because of using ROS, the mini shake robot reused the control software from the VSR, ensuring consistent ground motion generation processes in the simulation and physical experiments.In the first validation experiment, the boundary curves from the simulation and physical experiments were close.We noticed that the real PETG PBR was more fragile and are working to identify the causes of the difference.During the experiment, we observed high-frequency mechanical vibrations along the vertical direction (perpendicular to the ground motion direction) on the mini shake robot.Additionally, because the PETG PBR was downscaled and we know that smaller PBRs are more fragile (see Section 5.1), the fragility of the downscaled PBR may be more easily affected by ground motion noise.Therefore, this scale concern warrants further experimentation to validate the VSR using full-scale testing and rock material.To verify the consistency of the overturning results regardless of the PBR scale, future work should repeat the same experiment with 3D printed models at various scales.Because iterative algorithms for solving MLCP may not always guarantee convergence or unique solutions, resulting in variations in rock responses (Veeraraghavan et al., 2020), future work should investigate solutions to mitigate such uncertainties.For example, following Purvance et al. (2008), we can use the Monte Carlo method to construct a probabilistic model and to quantify the uncertainties.
In the second mini shake robot experiment, the fragility anisotropy results from the mini shake table were consistent with the results from the VSR, providing strong evidence of qualitative consistency in fragility analysis.The same pattern of boundary curves was observed in both simulation and physical experiments-the initial PBR orientation of 0°was more fragile than 270°.Given that the edge of the Double Rock PBR appeared well-defined and perpendicular to the Table 1 Large displacement statistics and ground motion correlation analysis mass center, inertia matrix, and geometry.In our experiments, it took only a few seconds up to a minute to finish one overturning experiment in the VSR.In DEM, the computational time increases with the complexity of geometry, discretization, and contact stiffness.Finishing one overturning experiment in DEM usually takes a few minutes to hours.The rapid deployment of the VSR facilitates qualitative analysis of PBR fragility, which aids in field data collection and assessment such as searching for the most fragile PBRs in the field.

Overturning and Large Displacement Studies
The VSR offers an integrated solution for studying the dynamics of both overturning and large-displacement processes.
From the overturning experiments of cuboids, the relationships between the dimension parameters and fragility were consistent with the previous analytical solutions (Purvance et al., 2008;Anderson et al., 2014).The VSR advances the study of PBR overturning dynamics in many aspects.For example, the overturning experiments indicated that the PBR overturning dynamics were more complex than previously understood.Specifically, the PBR overturning responses were found to vary from ground motion directions and the presence of surrounding pedestals.
Most of our experiments employed single-pulse cosine displacement ground motions.First, such ground motions were easy to synthesize, enabling us to comprehensively cover the configuration spaces of PGA and PGV/PGA.Second, while PGA and PGV/PGA are commonly used as representative ground motion descriptors in PBR studies, recent concerns have emerged regarding their use in PBR studies.When two seismometer data records share the same PGA and PGV/PGA values but differ in other properties, they may produce different PBR overturning responses.To mitigate such ambiguity, we adopted single-pulse cosine displacement ground motions.
The VSR's ability to support realistic terrain models allowed for the seamless study of large-displacement dynamics.This functionality of large-displacement analysis enabled the trajectory prediction of overturned PBRs, with potential applications in rockfall hazard zoning and the study of rocky slope development.Additionally, the large-displacement analysis revealed the relationship between large displacement statistics and ground motions for PBRs.The PGA and PGV/PGV positively correlated with mean trajectory length, mean largest velocity, and mean terminal distance.When the overturning dynamics provide upper-bound ground motion constraints by studying fragile PBRs, the large displacement dynamics provide lower-bound ground motion constraints by studying overturned PBRs.For example, given a known trajectory of an overturned PBR, lower bounds of PGA and PGV/PGV can be inferred from Fig. 17b, c.Together, the overturning and large displacement dynamics form complementary methods to refine ground motion estimation.

Limitations and Future Work
In future research, we will investigate several critical aspects to enhance both the VSR functionality and scientific insights.First, the previous study by Veeraraghavan (2015) demonstrated that 3D PBR fragility results are more sensitive, indicating higher fragility, compared to their 2D counterparts.Looking ahead, the development of the VSR that incorporates 3D ground motions will be a significant advancement beyond the existing 1D ground motion assessment methods.Second, our future work should align with the workflow delineated by Rood et al. (2020), focusing on the exploration of additional PBRs within the Double Rock site.Expanding our scope to encompass a broader array of PBRs has the potential to enhance the accuracy of ground motion estimation.Because of the significant role of contact geometry in rock response, future work should thoroughly investigate the VSR's ability to model dynamic processes involving complex contact surfaces.Such modeling would include complex interface geometry, variable properties and rheology, and evolution of geometry and properties with continued loading.Fourth, the validation of the large displacement process is not addressed in this study.Future research should focus on exploring this aspect.
Additionally, the values of physics parameters, such as friction and restitution coefficients, play a pivotal role in dynamics simulations within physics engines.Our experiments observed that these parameters produce nonlinearity in the PBR responses to ground motions.For example, we found several thresholds for friction coefficients.Within these thresholds, the friction coefficient displays various nonlinear properties.Generally, a significantly high friction coefficient made PBRs more fragile compared to a very low coefficient.However, within certain threshold ranges, the influence of varying friction is less pronounced.To quantitatively measure this nonlinearity, we recognize the necessity for additional experiments in our future research endeavors.
Future work could bypass Gazebo and directly build a VSR in the Bullet physics engine.Gazebo simplified the simulation configuration and provided perception and control packages compatible with ROS.However, when passing configuration parameters to the Bullet physics engine, Gazebo reduced the number precision of some parameters, such as the bits for floating point numbers.This reduction in the number precision presents challenges in calibrating contact properties.Building a VSR directly in Bullet allows complete control of the configuration parameters and aids in contact physics calibration.Additionally, as a previous study has used Bul-let to simulate the crushing process of granular materials (Zhu and Zhao, 2019), modeling PBR simulations in Bullet allows for the study of how overturned rocks are crushed along a trajectory, providing insights into more complex modeling of rocky slope development.Simulating the crushing process also enables the examination of the impact of PBR deformation on fragility (Saifullah and Wittich, 2021).Generally, physics engines hold promise in various physics-based scenario simulations, including testing of dynamic rupture model outputs (e.g., Lozos et al., 2015).
Future work should use the VSR to build ground motion models for PSHA.For the Double Rock site, for example, we can examine the effects of the surrounding pedestals and ground motion directions on hazard curves.This study has demonstrated a forward model of large displacement dynamics, which are useful for rockfall prediction and ground motion study.At the same time, we are conceptualizing the idea of using the tool to trace the origin of overturned rocks.For example, we can simulate a large number of trajectories from various initial conditions to build a probability heatmap of the origin of the overturned rock denoted in the blue box in Fig. 15a.Using the VSR, we can also simulate a large number of PBRs with various shapes and dimensions to study the effects of ground motions on PBR distributions.
Combined with our research of using robots and machine learning for rock detection and mapping, the VSR presents a paradigm of rock detection-mappinganalysis for automated geoscience.In the future, the VSR could be installed on a companion computer of a UAV that is also developed using ROS.Once the UAV detects and maps a PBR, the VSR can rapidly analyze the PBR fragility, facilitating field data collection.

Conclusion
The development of the VSR has demonstrated the potential of using robotics and physics engines for studying PBR dynamics.The advances in simulating PBR overturning and large-displacement processes by the VSR provide valuable information for seismic hazard analysis, PBR fate study, rocky slope development, and rockfall prediction.Validation experiments confirmed the good performance of the velocity controller in the VSR.To validate the overturning dynamics, we compared the overturning results from the VSR with those from previous experiments and the mini shake robot.The VSR produced consistent results with previous 2D studies in the cuboid overturning experiments.The overturning response diagrams suggested that ground motion directions had complex effects on PBR fragility.We investigated the effects of surrounding pedestals on overturning responses, which had been under-explored in previous studies.The effects of the surrounding pedestals generally reduced the PBR fragility compared with flat terrains and varied with ground motion directions.For the study of the large-displacement process, we conducted 2500 experiments with different ground motions and plotted 692 trajectories where the PBR was toppled.Correlation analysis showed that the ground motions positively correlated with large displacement statistics such as mean trajectory length, mean largest velocity, and mean largest terminal distance.The rapid deployment of the VSR facilitates qualitative analysis of PBR fragility, aiding in field data collection and assessment.As a result, the VSR provides a screening method to identify the most fragile PBRs in the field for more detailed dynamics analysis.Overall, the VSR represents a significant step forward in studying PBR dynamics, providing valuable insights for researchers and practitioners.

Figure 2
Figure 2 Simulations loops of (a) Bullet physics engine and (b) discrete element method (DEM).In both Bullet physics engine and DEM, states of objects are updated in each sequential discrete simulation loop with a fixed small time step.Their mechanisms to process collision are fundamentally different.

Figure 3
Figure3Software architecture of the virtual shake robot (VSR), composed of Robot Operating System (ROS), Gazebo simulation toolbox, and Bullet physics engine.Developing robot software based on ROS allows the reuse of robot software for both virtual and real shake robots, ensuring the ground motion generation processes remain consistent between the two environments.
overturning experiment, which included the final status of overturned or balanced after the single-pulse ground motions.The automation program was at the highest level in the hierarchical control module.The automation program passed high-level control signals {(PGV, PGA)} to the middle-level motion interpreter and trajectory planner.Then the middle-level trajectory planner passed control signals (velocity commands) to the lowlevel PID velocity controller, which generated force for the prismatic joint.

Figure 4
Figure 4 VSR with (a) flat and (b, c) realistic terrains.Unpiloted Aircraft System (UAS) and Structure from Motion (SfM) produced the full-scale realistic terrains in (b) Double Rock, California and (c) Granite Dells near Prescott, Arizona.Arrows indicate ground motion directions.

Figure 5
Figure 5 (a) Control module and (b) automation workflow of the VSR.

Figure 6
Figure 6 3D geometric modeling of terrain and PBR mapped by UAS and SfM.The site of Double Rock is located in southcentral coastal California.

Figure 7
Figure 7 Velocity plots from (a) single-pulse cosine displacement ground motion and (b) realistic ground motion.The PID velocity controller generates translational force to actuate the pedestal of the VSR based on the desired velocity from a trajectory planner or a record of a real-world accelerometer.The actual velocity is measured from the pedestal of the VSR.The desired velocity and actual velocity are overlain in panel (a).(b) A sudden disturbance in the actual velocity, as highlighted in the box, was caused by the overturned PBR but quickly rectified by the PID controller.
maximum velocity of 0.5 m/s with up to 2 kg PBR payload.Overall, the mini shake robot provides a low-cost, open-hardware, and open-software platform for earthquake research and education.

Figure 9
Figure 9 Mini shake robot from (a) side view and (b) top-down view.(b) 3D-printed PLA PBR is placed on the pedestal.Reprinted from Chen (2022) with permission.

Figure 10
Figure 10 Overturning response diagram of 3D-printed Double Rock PBR.Blue curve represents logistic regression results from the mini shake robot.Red curve represents logistic regression results from the VSR.Reprinted from Chen (2022) with permission.

Figure 11
Figure 11 Overturning response diagrams of 3D printed Double Rock PBR from initial orientations of yaw 0°and 270°.The red and blue dots represent overturning responses of being toppled and balanced after single-pulse cosine ground motions, respectively.The curves indicate the boundaries of the overturning responses.Reprinted from Chen (2022) with permission.

Figure 12
Figure 12Overturning response diagrams from cuboids with different dimensions on a flat pedestal.Ground motions are along x direction in the x×y×z dimension.Red and blue dots represent overturned response and balanced response, respectively.Horizontal axis represents the peak ground acceleration (PGA) in gravity constant.Vertical axis represents the ratio of peak ground velocity to peak ground acceleration (PGV/PGA).Unit of PGV/PGA is second.The blue dots within the toppled zones represent the PBRs being pulled back to balanced states by the returning ground motions.

Figure 13
Figure 13 Overturning response diagrams from Double Rock PBR on a flat pedestal.(a) Arrows indicate the equivalent ground motion directions of the different initial PBR orientations.

Figure 14
Figure 14 Overturning response diagrams from Double Rock PBR on the realistic pedestal mapped from UAS and SfM.(a, c) Arrows indicate directions of the single-pulse cosine ground displacement motions.Inset at (a) shows a zoom-in of the PBR on the realistic terrain.

Figure 15
Figure 15 Large-displacement experiment on the Double Rock site from ground motions in the yaw 0°direction.Arrow in panel (a) shows the direction of the single-pulse cosine ground displacement motions.The PBR trajectories in panel (b) are relative positions to the terrain.Colors indicate absolute velocity of the PBR.(c-f) PBR trajectories from different PGA ranges, where N represents the number of trajectories.

Figure 17
Figure 17 Correlation analysis between large displacement statistics and ground motions.(a, c, e) Correlation between PGA and mean trajectory length, mean largest velocity, and mean terminal distance.(b, d, f) Correlation between PGV/PGA and mean trajectory length, mean largest velocity, and mean terminal distance.Red dash-dotted lines result from least-square linear regression.