第27卷 第12期
岩石力学与工程学报 V ol.27 No.12
2008年12月 Chine Journal of Rock Mechanics and Engineering Dec.，2008
Received date ：2008–08–07；Revid date ：2008–09–16
Corresponding author ：JIANG Jingcai(1962–)，male ，Ph. D.，graduated from Hohai University in 1983. He is now an associate professor ，and his main rearch interests include static and dynamic slope stability analys ，shear strengths of landslide soils ，risk asssment of earthquake-induced landslides ，and soil bioengineering. E-mail ：kushima-u.ac.jp
IDENTIFICATION OF DEM PARAMETERS FOR ROCKFALL SIMULATION ANALYSIS
JIANG Jingcai 1，YOKINO Kazuyoshi 2，YAMAGAMI Takuo 1
(1. Department of Ecosystem Design ，Institute of Technology and Science ，University of Tokushima ，Tokushima 770–8506，Japan ；数据寄存器
2. Kiso-kentsu Consultants Co .，Ltd .，Tokushima 779–3120，Japan )
Abstract ：Distinct element method(DEM) provides a powerful analytical tool for modeling rockfall behaviors ，i.e. trajectories ，velocities and energies of falling rocks. Better DEM prediction of the motion of the falling rock requires realistic input parameter values ud for the analysis. A methodology is developed in which the DEM parameters are determined from back analysis of a rockfall trajectory obrved in field rockfall experiments or traced by site investigation of natural rockfall events. A number of techniques and strategies are propod to ensure the efficiency and robustness of the solution procedure. The method is also extended into non-homogeneous slope cas ，by applying repeatedly the prented procedure to each of the homogeneous gments along a given rockfall trajectory. Laboratory rockfall experiments on small model slopes are carried out to verify the propod method. Ca studies on well-documented rockfall records are also prented to show good performance of the method for both homogeneous and non-homogeneous slopes.
Key words ：slope engineering ；rockfall ；distinct element method(DEM)；back analysis ；DEM parameters ；rockfall trajectory
CLC number ：P 642.22 Document code ：A Article ID ：1000–6915(2008)12–2418–13
滚石离散元数值模拟的参数反演
蒋景彩1，能野一美2，山上拓男1
(1. 徳岛大学 研究生院技术科学研究部建设工学科，日本 徳岛 770–8506；
2. 基础建设设计咨询有限公司，日本 徳岛 779–3120)
摘要：离散元法(DEM)是数值模拟崩塌落石运动的一个有力方法，但确定离散元输入参数是个难题。提出一个基于已知落石轨迹的离散元参数的反演计算方法。落石轨迹可以通过现场落石试验或者野外实地调查实际落石过程获得。该方法的特征为反演得到的离散元参数可以间接反映边坡与落石冲击时的变形以及植被缓冲效应。建议的几个保证解的有效性得到验证，稳定性措施效果良好。室内落石试验结果和现场落石工程实录验证证明了该算法的有效性和实用性。
关键词：边坡工程；落石；离散元法；反演计算；DEM 参数；落石轨迹
1 INTRODUCTION
The purpo of this paper is to develop a practical
and uful method for determining input parameter values for distinct element method(DEM) simulation of two-dimensional rockfall motion on slopes. DEM provides a powerful analytical tool for modeling
第27卷第12期JIANG Jingcai，et al. Identification of DEM Parameters for Rockfall Simulation Analysis • 2419 •
rockfall behaviors，such as trajectories，velocities and
energies of falling rocks[1]. A distinct element analysis
of falling rocks requires known values of five
parameters as shown in Fig.1. The parameters are
normal spring constant K
nich liebe dich
，tangential spring constant
K
s ，normal damping ratio η
n
，tangential damping ratio
ηs，and friction coefficient µ. Although DEM numerical simulation of the rockfall behavior has been widely carried out in practice，an effective and efficient method has not been established as yet for determining the input parameters which govern rockfall motion and may differ greatly from one slope to another due to the difference in geometry and soil properties involved in the slope[2].
Fig.1 DEM contact model between slope and falling rock
Asssment of the dynamic DEM parameters is rather difficult becau rockfall motion depends on many factors including the slope geometry and falling blocks，deformation and strength characteristi
cs of slope materials，and vegetation. Another difficulty associated with the determination of above-mentioned parameters is that they cannot be obtained directly from rockfall experiments. A rockfall experiment can produce a trajectory of the ud falling rock and its velocity distribution，but it cannot directly provide DEM parameters. In particular，it is difficult to asss the magnitude of the parameters considering the effects of deformation of the slope ground when the falling blocks collide against the slope and tho of trees and other plants being on the slope.
In this paper，a methodology is prented in
which the DEM parameters(K
n ，η
n
，K
s
and η
s
) between
a falling block and a given slope are identified from
back analysis of either natural or artificial rockfall
events. The propod method，though in an indirect
manner，can evaluate the parameters accounting for
the effects of deformation of the slope ground and/or
collision between trees and blocks. In a distinct
element analysis，friction coefficient µ is also needed；
and basically this parameter should be included in the
back analysis procedure as well. However，the effect
of µ on the back-analyzed parameter values is found
so small that it is excluded from the parameters to beassault
back analyzed. Usually，the u of an empirical value
for µis enough for determining the other four
parameters(K
n
，η
n
，K
s
and ηpropo是什么意思
s
) in the propod method.
The central idea of the prent study is to conduct
a rockfall experiment at the site in consideration and to
monitor its falling trace using video cameras. After
digitizing the video records，the location of the falling
block ud in the experiment is measured at different
instants. The back analysis method is then contrived to
identify the DEM parameters bad on the trajectory
by using a nonlinear optimization technique. It is
noticed that the rockfall trajectory obtained indirectly
reflects the deformation of slope ground and thevision
effects of rockfall impacts to ground and trees on
behavior of the falling rock. It must also be
emphasized that given a rockfall trajectory as a result
of the site investigation of a natural rockfall event，the
propod method can also be ud without performing
in-situ rockfall tests.
2 METHODOLOGY
As mentioned in the introduction，the prented
back analysis method requires a known rockfall
trajectory. This trajectory of the falling rock can be
obtained from two ways. One is carrying out in-situ
full-scale rockfall experiments and catching the
position of the falling rock by digitizing the video
records. The other is tracing a natural rockfall event by
detailed site investigations.
A substantial difference in the two cas is that in
the first one，the elapd time to reach each obrvation
• 2420 • 岩石力学与工程学报 2008年
point on the trajectory from the start of a rock-block falling test is known(Fig.2)，whereas the elapd time is unknown at all in the cond ca(Fig.3).
Ca 1 Time-kown trajectory
(X i+1,Y i+1) at t=t i +
∆t (X i ,Y i ) at t=t ipax
Fig.2 Description of a rockfall trajectory in ca that
elapd time is known
Ca 2 Time-unkown trajectory
) at point j+1
Fig.3 Description of a rockfall trajectory in ca that elapd
time is unknown
Either of the following objective functions is minimized in order to arch for a t of parameters which reproduce the same behavior of the rock fall trajectory. For cas in which the elapd time as well as rock fall trajectory is known ，the objective function can be expresd ：
22
ob cal ob cal 1()()n
i i i i i U X X Y Y =⎡⎤=−+−⎣⎦∑，，，， (1)
where U is the objective function ；ob i X ， and ob i Y ， are obrved coordinates of the center of gravity of the falling rock at time i (i = 1，
2，…，n )；and cal i X ， and cal i Y ， are computed coordinates at the same time i
from DEM numerical simulation.
For time-unknown cas ，the objective function can be expresd as follows ：
()
2
ob cal 1m
j j
j U Y Y ==−∑，， (2)
where ob j Y ， is measured y coordinate of the center of gravity of the falling rock at point j (j = 1，2，…，m )；and cal j Y ， is computed y coordinate at the same point
j from DEM numerical simulation.
At first ，the four parameters are tried to be back analyzed at once using a single optimization technique bad on nonlinear mathematical a
modified Marquardt method [3]). However ，such a
simple solution procedure can not yield convergent ，satisfactory back analysis results. Presumably ，the number of parameters to be back analyzed ，four ，is considered to be too many in such complex phenomena as rockfall. To solve this problem ，it is of interest to note the study of T. Ohmachi and Y . Arai [4]，who introduced the repulsion coefficient e ，as shown in the following equation ：
exp e ⎛
⎞
⎜=⎜⎝
(3) where M is the mass of the falling rock block.
The authors examine the effects of this coefficient
on computed trajectories of the falling rocks. It is found that the computed trajectories are almost the same when the values of e ，K n /K s and ηn /ηs are the almost identical respectively ，even if the value of each individual parameter is different.
As an example ，Fig.4 shows four trajectories of a global falling rock with a mass of 1 680 kg and a radius of 1.0 m ，which are computed by DEM using almost the same t of e ，K n /K s and n s /ηη but different combinations of K n ，K s ，n η and s η values
(Table 1). It can be en from Fig.4 that the trajectories and the velocity distributions are almost the same ，though the combinations of K n ，K s ，ηn and
ηs ud for the four cas are different.
Using the study results shown in Fig.4 and having examined various situations by trial and error ，an
(X i +1，Y i +1) at t i + ∆t (X i ，Y
i ) at t i ) at point j + 1
第27卷 第12期 JIANG Jingcai ，et al . Identification of DEM Parameters for Rockfall Simulation Analysis • 2421 •
X /m
(a) (b)
Fig.4 Rockfall trajectories and velocity distributions simulated
using the parameters values in Table 1
Table 1 Values of DEM parameters ud for cas 1 to 4
Ca No.
K n /(N ·m －1) K s /(N ·m －1) ηn /(N ·s ·m －1) ηs
/(N ·s ·m －
1)
e
1 6.00×107 9.50×107 1.60×104 1.00×104 0.955 0
2 3.00×108 4.80×108 4.00×104 3.00×104 0.955 3
3 5.00×108 8.00×108 5.00×10
4 4.00×104 0.956 74 4.00×109 6.50×109 1.50×10
5 1.00×105 0.954 1Ca No.
K n /K s
ηn /ηs µ ∆t /s
1 0.63
2 1.600 0.2 1.0×10－5
2 0.625 1.33
3 0.2 1.0×10－5
3 0.625 1.250 0.2 1.0×10－5
4 0.61
5 1.500 0.2 1.0×10
－5
effective and practical solution procedure that can provide stable and appropriate solutions ，who flow- chart is shown in Fig.5，is achieved. This procedure includes an artificial neural network(ANN)
[5，6]
for
evaluating e ，K n /K s and ηn /ηs values and a double optimization technique for determining individual values of K n ，K s ，ηn and ηs . It should be noted that friction coefficient µ is assumed to be known ，as stated previously.
In order to evaluate e ，K n /K s and ηn /ηs values corresponding to a given rockfall trajectory ，four-layer networks shown in Fig.6 are ud in this study. Input data for the ANN are the coordinates of the center of gravity of the ud falling rock block at a number of points on the trajectory ；and output data are the values of e ，K n /K s and ηn /ηs . The implementation of a neural network model mainly involved training data generation and network construction.
Training data generation is an esntial step of the
Fig.5 Flowchart for the back analysis procedure bad on
ANN and double optimization technique
Fig.6 Neural network ud to estimate values of e ，K n /K s and
ηn /ηs for a given trajectory of rockfall
ANN model. Firstly ，a t of K n ，K s ，ηn and ηs values is appropriately specified ；and using the values a
Y /m
Ca 1
Ca 2 Ca 3 Ca 4
F a l l /m
Velocity/(m ·s －1
)
• 2422 • 岩石力学与工程学报 2008年
distinct element analysis is carried out to obtain a rockfall trajectory. Then ，the coordinates of a number of points on the obtained trajectory and the corresponding values of the parameters(e ，K n /K s and
ηn /ηs ) are taken as input data and the desired output
data. The above procedure is repeated to obtain a sufficient number of training data ts.
樱知叶留学
As a four-layer neural network has been specified above ，network construction is then carried out to determine the number of units of two hidden layers and to train the network using the training data genera- ted above. The number of units of the input layer is decided by the number of the lected points on the rockfall trajectory ；and only three units are necessary for the output layer. The number of hidden layer units is adjusted in the training process in order to make the total error minimal ，and to make the network capable of predicting a reliable output for an input included in the training data. However ，at prent ，there is no rule for determining the optimal number of units in the hidden layers. A common practice is to u a trial and error method in which training is repeated veral times ，starting with a few units for each hidden layer ，and then the unit number is incread. In this process ，the total error between the current output and the desired output is monitored ；and the training continues until the error is no longer reduced significantly by increasing the number of hidden layer units. The above approach is also employed in this study.
提出用英语怎么说
In order to verify the performance of the constructed neural network ，some of the training data ts 用纯牛奶洗脸
are usually lected as the testing data which do not enter into the data ud for the training. The satisfactory completion of training for the network signifies that the network solutions from the testing data should agree well with the desired outputs. After validation ，the trained network can be ud to predict values of e ，K n /K s and ηn /ηs ，i.e. prenting the given trajectory into
the trained ANN will yield the corresponding values of e ，K n /K s and ηn /ηs at the output layer.
Using the values of e ，K n /K s and ηn /ηs obtained from the ANN model and assuming an value of K n ，the
老友记第十季
associated values of K s ，ηn and ηs can be calculated through Eq.(3). The t of four parameters so obtained is employed as initial values for the optimization problem in terms of either Eq.(1) or (2). The modified Marquardt method is then ud to solve the problem. This process is enclod with a dotted line ，indicated as Optimization loop I in Fig.5.
In general ，the accuracy of solution of optimization problems depends significantly on the initial values. It is therefore common to solve the problem with veral different initial values ，and to adopt the best result among them. However ，this approach is rather cumbersome. Thus ，to avoid this tedious task ，a double optimization method has been constructed by introducing an additional l
oop in which preferable initial values are arched automatically. In this method ，the solution of the optimization problem can be regarded as a function of initial values ，and a simplex method [7] is adopted to solve the prent problem. The process of the double optimization technique is included in Fig.5 as Optimization loop II.
3 APPLICATION TO LABORATORY
ROCKFALL TEST RESULTS
Results of rockfall test on small scale model slopes are ud to verify the propod method. Two types of model slopes ，deformable and rigid(non-deformable)，are prepared in order to e whether or not the prented method can consider and reflect the deformability of the slope. The two model slopes ，a concrete slope and a decompod granite soil slope ，are both straight line patterns with an angle of 40°.
Tests are carried out on cubic blocks with sizes of 4 cm×4 cm×6 cm. Falling blocks are videotaped and rockfall trajectories are established by digitizing the records. An attempt is made so that almost two- dimensional rockfall motion within a vertical ction of the slope takes place ，as shown in Fig.7.
After sliding on the slope for a certain distance ，the block collides against an obstacle fixed on the slope and then bounces off the slope. Trajectories
