1. 引言
我国西部地区地质灾害频发，以四川为例，每年都会发生滑坡、泥石流、崩塌落石等灾害，频发的落石灾害严重危害公路、铁路行车安全。近年来，我国多采用隧道明洞结构以及各种柔性防护网对落石灾害进行防护，在防落石明洞结构设计过程中开始考虑落石冲击力的影响，但是落石冲击力目前没有一个非常明确的计算公式来定量计算。目前全世界共有5种落石冲击力计算方法，分别来自隧道设计技术手册 [1] 、公路路基设计规范 [2] 、西南交通大学(杨其新等) [3] 、日本道路协会 [4] 和瑞士学者 [5] ，都是基于试验的基础上提出了落石冲击力计算方法，每种方法的计算原理和考虑的因素不一样，导致落石冲击力计算差异较大；越来越多的研究学者开始用离散元方法来研究落石冲击力，并用颗粒流软件PFC3D来模拟落石冲击过程。裴向军，刘洋，王东坡 [6] 等通过有限元数值模拟落石冲击砂土垫层，得出砂土垫层可以有效减小落石对棚洞顶板冲击效应的结论，并经过优化得到合理垫层厚度；杨建荣、白羽、杨晓东 [7] 等提出了一种由金属柔性网和弹簧撑杆组成的柔性棚洞并采用有限元方法进行数值模拟，计算动力响应；王静峰、赵鹏、袁松 [8] 等基于大尺度模型试验和有限元数值模拟相结合的方法，对EPS和砂土垫层组成的复合垫层钢棚洞抗落石冲击性能展开了研究，得出复合垫层可以有效提高钢棚洞对落石的抵抗能力的结论；吴佳楠 [9] 同样以钢棚洞模型冲击试验和有限元数值模拟相结合的方法，对钢棚洞在EPS-细砂复合垫层对落石的冲击效应展开研究，得出峰值加速度、最大应变与下落高度成正比的结论；王永东，周天跃，柴伦磊 [10] 等基于能量分析方法对落石冲击棚洞结构进行研究，考虑垫层弹塑性、垫层与落石之间摩擦力、棚洞自身耗能三方面影响因素，得到了冲击深度以及冲击压力的计算公式，还得出垫层成坑变形耗能最大达到50%、垫层与落石摩擦力耗能次之达到30%以及棚洞自身耗能占比20%的结论；闫帅星，何思明，李新坡 [11] 等基于冲击动力学、修正的混凝土塑性特性及永久变形分阶段弹性接触理论，构建了钢筋混凝土棚洞在落石冲击作用下动力响应的一般解析方法，并通过有限元进行验证。
作者以离散元理论为基础，用ITASCA公司开发的颗粒流PFC3D软件建立落石冲击结构顶板的颗粒流模型，通过内置fish语言提取结构顶板在落石冲击下的最大冲击力。在文献 [12] [13] 中绘制落石冲击力时程曲线。本文在此基础上采用SPSS软件进行数据处理，拟合出落石冲击力在多因素影响下的多元线性计算公式。
2. 落石冲击模型介绍
2.1. 落石冲击结构顶板的离散元模型
落石采用单个颗粒小球模拟，回填土采用随机分布的颗粒小球模拟，如图1所示。
Figure 1. Calculation model of rockfall impact
图1. 落石冲击计算模型
2.2. 计算工况如表1所示
本次数值模拟工况组合共计5*6*6*7 = 1260种，如表1所示。
Table 1. Combination of working conditions [12] [13]
表1. 落石冲击顶板结构离散元数值模拟工况组合 [12] [13]
2.3. 离散元数值模拟微观参数设定
本次数值模拟微观参数设定采用文献 [12] [13] 中的数据，具体信息见表2。
Table 2. Micro-properties of the model element [12] [13]
表2. 模型单元细观参数 [12] [13]
2.4. 离散元数值模拟微观参数标定
本次离散元数值模拟微观参数标定采用文献 [12] [13] 中的方法和结果，落石、回填土、结构的微观参数标定结果如表3所示。
Table 3. Macro-parameter calibration results of backfill [12] [13]
表3. 回填土宏观参数标定结果 [12] [13]
3. 多元线性回归理论
3.1. 多元线性回归的定义和模型
在回归分析中，如果一个因变量受两个或两个以上的自变量的影响，那么因变量与这些自变量就存在多元线性关系。研究因变量和多个自变量的线性关系，就称为多元线性回归 [14] 。
多元线性回归的数学模型是：
(3-1)
式中，P是拟合变量，Xi是各影响变量即自变量，B0是多元线性回归常数项，Bi是各个自变量的偏回归系数。
3.2. 线性回归数据处理
多元线性回归方法介绍如下：
1) 建立多元线性回归模型；
2) 求解各个自变量的偏回归系数Bi；
3) 回归模型检验，包括相关性检验、相关系数R、R2检验、回归系数t检验，显著性F值检验等等 [13] [14] 。
4. 落石冲击力线性回归
落石冲击力深入研究的最主要目的有两个：第一，确定落石冲击力的计算公式，为防落石结构设计提供落石冲击荷载的设计依据；第二，确定最佳合理回填土厚度，既减轻结构自重，经济合理，又能有效缓解落石冲击效应，达到防灾减灾的目的，确保山区落石灾害地区公路、铁路运营安全。
文献 [12] [13] 得出0.6 m厚回填土就能降低50%以上的冲击力和最佳回填土厚度是2~3 m的结论，因此，本文以回填土厚度为主要自变量来划分分段回归区段，将回填土厚度t = 0、(0,2)、[2,3]、(3,5] m作为回归分段。本文只介绍分析回填土厚度为2~3 m时的线性线性回归方程。
4.1. 回归模型分析
回归模型R2数据分析：
从表4中可以看出，模型1只考虑常数项和自变量一次项可以解释落石冲击力变化的70.3%；模型2引入自变量二次交互影响项，解释落石冲击力增加18.4%增大至88.6%，说明自变量的二次交互项对落石冲击力的影响很明显；同理模型3和模型4分别引入三次交互项，R2分别增加1.9%和1%，说明三次交互项和四次交互项对落石冲击力影响不太显著。后续结合各交互项显著性进行进一步分析。
Table 4. Regression model analysis [13] [14]
表4. 回归模型分析 [13] [14]
4.2. 方差分析
方差分析结果解析如下。
Table 5. Anova analysis of regression model [14]
表5. 回归模型方差分析 [14]
从表5可见，四个模型Sig值都小于0.05，说明四个回归模型都是显著的其中模型二的F最大，说明模型二拟合度最好。
4.3. 回归系数显著性分析
表6中列出了各个模型中各个变量的回归系数以及t检验P值，从表中数据可以看出，几个模型中各个变量的回归系数P值均小于0.05，说明回归系数是有效的。
4.4. 回归系数相关性检验
表6中列出了各个模型中回归系数相关性检验结果，表中数据显示，各个自变量与因变量之间是相关的，而且各个自变量之间是相互独立的，不会相互影响，回归模型有效。
4.5. 落石冲击力回归方程
结合4.1~4.4分析结果可知，回填土厚度为2~3 m时，四个模型都是显著的，四个自变量对落石冲击
Table 6. Significance analysis and correlation test of regression coefficient
表6. 回归系数显著性分析及相关性检验
力变化的影响都很明显，二次交互项全部是落石冲击力的显著性影响因素，三次交互项只有HθW显著，四次交互项不显著，对落石冲击力影响较小可以忽略。模型2拟合度最好，模型3的精度R2最大，且保留HθW会使精度提高19%，综合考虑选取模型3并且保留三次交互项HθW，最后的拟合精度R2 = 0.905。落石冲击力计算公式如下：
即
5. 结论
本文在离散元颗粒流PFC3D数值模拟落石冲击结构顶板并提取顶板结构最大冲击力数据的基础上，采用SPSS软件进行数据处理，基于多元线性回归理论与方法，对1260组计算结果进行多元线性回归处理，经过仔细对比分析各回归模型及其回归系数的显著性和相关性，最后选取最佳模型，得出落石冲击力的计算公式。同时，考虑了各影响因素对落石冲击力的交互影响作用，回归模型显著，拟合度较好，精度较高。最后可以得出如下结论：
1) 自变量H、W、θ、t对落石冲击力的影响都是显著的；这个结论与文献 [12] [13] 的结论是相一致的；
2) 四个自变量一次项、二次交互项对落石冲击力影响显著，个别三次交互项显著，充分说明各个自变量在对落石冲击力影响显著的同时，还对落石冲击力存在交互影响作用，交互项在落石冲击力计算公式研究过程中要予以考虑，不可以忽略不计。
参考文献