中国是受台风影响最严重的国家之一，据统计，我国大约80%的省（市、区）都会不同程度地受到台风的影响^[1]，台风带来的强风会直接影响树木的稳定性，导致树木倾斜或倒伏，由此造成的树木损失不可忽视^[2]. 2023年台风“杜苏芮”使泉州市16万多株园林树木发生倒伏、折干、断枝等. 2024年9月，台风“贝碧嘉”吹倒上海共青森林公园近千株树木. 城市内被台风损毁的行道树，不仅影响到城市道路交通的正常运行，同时也对道路使用者的安全出行造成威胁；强风在植被上施加的风荷载传入土体，风荷载致使根土区发生松动甚至树木倾倒，可能造成滑坡^[3]，由此引发的次生灾害给人民生活、生产带来重大损失. 因此，研究树木在风荷载作用下的稳定性对减少风灾损失、提升城市绿化防护能力具有重要意义.

由于根系埋在土体中并不可见，再加上其锚固机制复杂，早期的研究将树的极限承载力与树冠结构、树干的几何特性及生物量分布等直接关联^[4-5]，缺乏对根土相互作用的考量. 近年来，Cannon等^[6]、傅胤榕等^[7]通过现场倾覆试验，探究树种、树体大小、破坏模式和土体含水率对树木稳定性的影响. 由于现场试验具有不可重复性，Zhang等^[8-9]、赵俊键等^[10]分别基于根系空间结构的实测数据和土体水文力学特征建立缩尺模型，使用3D打印技术对根-土相互作用开展物理模拟. 这些试验都为根系锚固机制的研究奠定基础，但由于根系结构观察的难度较大，相关定量分析研究进展相对迟缓^[11]. Blackwell等^[12]、Achim等^[13]相继提出根系锚固的计算模型，为植被根系在水平荷载下的分析提供理论依据. 徐华等^[14]采用MechRoot程序建立根土复合体模型，阐明根系形态与层次结构对根土复合体力学特性影响及固土作用机理. Danjon等^[15]通过现场测量将根系重构以探究坡上植被根系-土相互作用机理；Yang等^[16]使用有限元软件建模，对根系进行简化重构，探究树木拔起过程中根系的连续断裂情况. 但因为根系空间结构复杂无明显对称性，重构后的模型仍较为复杂，存在整体尺寸大而单根尺寸小的特点，考虑到计算时间成本的可行性，使用传统有限元分析所划分的网格对于单根而言，往往较为粗糙，且无法捕捉整个根系在大变形条件下的响应.

Winkler地基梁法可以避免有限元建模时对土体的网格划分，为上述问题提供解决思路. Liang等^[17]、Meijer等^[18]以桩基工程中的p-y曲线为基础，结合双曲正切函数进行简化，将土体等效为弹簧对竖向单根进行约束，探究其与土体的相互作用机理，并通过与试验结果进行对比，验证Winkler地基梁法的有效性. 在此基础上，Zhang等^[19]提出双线性软化假设，对迎风区的水平单根在大位移条件下的响应进行探究. 因为根系的几何特征复杂，且梁模型难以表征根的弯曲形状，目前此方法只在对单根响应的初探阶段，尚未运用到根系的模拟.

为了解决这个问题，提出基于Winkler地基梁法根系重构的方法，即将根系中所有单根都简化为水平或竖直向的细长圆柱. 采用此方法对A和B 2种不同空间构型的根系进行重构，通过超重力倾覆试验，比较重构前后模型的弯矩-转角图和破坏形式，验证重构方法的可行性，进而为植被根系的一体化建模提供理论基础.

根据根系的空间结构可以把根按照级数进行分类，其中，中央主根为1级根，由主根衍生的根为2级根，由2级根衍生的根为3级根，依此类推. Zhang等^[8]依据复杂根系的现场测量数据，按照如下法则进行物理模型的重构：1）综合考虑边界效应、尺寸效应、3D打印机分辨率等，将根系整体进1∶20缩尺（除特别说明，下文所涉及的尺寸均为缩尺后的大小）；2）考虑到3D打印机的分辨率，忽略根系直径小于1.2 mm的部分；3）保留所有根的几何弯曲特征，但除主根外，每个单根的直径沿轴向保持不变，且只用1.6、2.0、2.4 和 3 mm 4种直径表征，其对应的原型平均直径范围是<36、36~44、44~52、>52 mm；4）主根的直径沿轴向不作调整，从根部到根尖呈逐渐减小的趋势. 选取窄而深（模型A，图1(a)）和宽而浅（模型B，如图1(b)）的根系进行重构，得到模型Am（图1(c)）和Bm（图1(d)）. 为了使模型更适配Winkler地基梁法，便于今后的数值模拟和机理分析，在前期研究基础上，将模型根系进一步简化，提出的假设如下：1）所有根都为圆柱体，即忽略根的弯曲特征；2）1级主根竖直生长，且忽略其竖向分支；3）2级根为浅根，水平生长，且深度相同，即在xoy平面内（见图1），且不考虑20 mm以下的2级根；4）3级根沿z向竖直生长；5）只考虑1、2、3级根. 所涉及的相似律可参见表1.

将现场测量得到的根系三维空间坐标缩小20倍，得到1∶20的缩尺根系节点坐标，深度方向为z向，xoy平面为横向扩展面. 根据上述假设，按根的级数顺序构建根的简化模型. 如图2所示，对于主根，根据顶部端点E和底部端点F的坐标，构建线段EF₁，其中F₁的x、y坐标与E相同，z坐标与F相同. 以EF₁为轴线构建圆柱主根，根的直径取原模型中主根直径的平均值. 可将主根EF₁大致看作是 EF的投影. 为便于树干模型的安装，在顶部处将主根进行倒角处理，以达到顶部横截面积和树干适配.

如图3所示，基于假设，把所有横向浅根（2级根）都简化到同一平面xoy上，这些根的z坐标都相同，z值取简化前所有横向根浅根起点处深度的平均值z_h. 因为2级根从主根衍生而来，简化后所有2级根的起点坐标都为（x_E，y_E，z_h）,其中x_E、y_E为E点的x、y坐标. 对于根的另一个端点，类似地，保留x、y坐标不变，将z坐标变换为z_h. 以这2个点的连线为轴，构建圆柱，即为简化后的2级根. 实际上，这包含2个方面的简化：将曲根简化为直根，进一步将该根投影到平面z = z_h上.

如图4所示，对于3级根，须分为2类进行简化. 第1类是横向分叉GH，第2类则是竖向根IJ. 衍生该根的2级根轴线方程已知，故对于3级根起点G或I，保持该点的x坐标不变，便可在2级水平根上确定其起点坐标，即2、3级根的连接位置. 原始根系中大部分3级根与平面z = z_h的夹角呈现明显的双峰分布，主要集中在小于10°和大于45°这2个区间. 因此，选择15°作为分类阈值. 若3级根和参考平面z = z_h间的夹角小于15°，则认为其为横向分叉，此时使其终点H的x、y坐标保持不变，将z坐标变为z_h；若3级根和参考平面z = z_h间的夹角大于15°，则将其当作竖向根进行简化处理，此时令其终点J的z坐标不变，变换x、y坐标使其与I相同. 当2个端点确定后，基于平均半径构建3级根的简化圆柱.

简化后的根系Am、Bm如图1(c)、(d)所示. 由于单根的分布形式对植被的抗倾覆性能具有重要影响，将根系等分成8个区域以量化根的分布，实际风向如图5所示. 对于模型A而言，其横向根在简化前后的数量Ahs保持不变（图5(a)）；竖向根的数量Avs（图5(b)）在0°~45°区域内增加了1个，在315°~0°区域内减少了1个，这是因为竖向根在简化时是基于起始点的x坐标和简化后2级横向根的投影位置，故其简化后的位置有可能发生改变. 对于模型B而言，其横向根数量Bhs（图5(c)）在45°~270°、315°~0°区域内均减少，这因为横向根在进行简化时，由于须简化到同一平面，只选取20 mm深度（对应原型中的400 mm）的浅根，进而忽略20 mm以下的横向根. 同样地，因为不考虑2级竖向根，竖向根的数量Bvs（图5(d)）在135°~225°、270°~0°区域内有所减少，且在45°~90°区域内的1个根在简化后转移到相邻区域90°~135°内. 由文献[8]可知，主根和迎风区（315°~45°）及背风区（135°~225°）的横向根对根系的抗倾覆能力贡献较大，而简化后的根系可以大致保留这些特征根的位置.

在原型尺寸下，简化前后某深度h以上的根系总体积V分布如图6所示. 由于模型的缩尺比为1∶20，深度的相似系数为20，体积的相似系数为20³. 由图6可知，随着深度的增加，总体积逐渐增大. 模型A的深度为1 800 mm，为模型B的1.4倍；模型B的体积为1.62×10⁵ mm³，为模型A的1.7倍；由于2级横向浅根和中央主根的存在，2种不同根系的体积均在0~200 mm深度内快速增加，但差距不大. 到达400 mm深度时，模型B的累计根系体积明显高于模型A；当深度大于400 mm时，2个根系的累计根系体积增长变缓. 在较浅深度范围内，简化后的模型Am、Bm累计根系体积较模型A、B更大，这是因为0~400 mm内的2级横向根被简化到同一深度. 但因为只保留1、2、3级根，在深度增加时，简化后的根系体积有所减小.

使用丙烯腈丁二烯苯乙烯（ABS）塑料将根系进行3D打印，打印机型号为 Fortus 250mc 3D. 研究表明，ABS塑料具有和现场植被根系相近的脆性特征^[9]，在单轴抗拉和三点抗弯试验中，弹性模量E和抗拉强度T均与直径呈负幂律关系，且范围与真根吻合较好（见图7），能较为准确地模拟根系的力学特性.

为了模拟现场土的颗粒级配和抗剪强度的特性，使用的试验土为HST95砂和A50粉土混合物，质量配比为7∶3，级配曲线如图8所示，试验土不均匀系数为15.1，曲率系数为2.4，级配良好. 不同正应力$ \sigma {'} $下试验土与现场土的最大摩擦角$ \varphi_{\mathrm{max}} $如图9所示，由图可知，在低应力条件下，最大内摩擦角均有显著增加，随着应力的增大，模型土的最大内摩擦角趋近于临界内摩擦角36.4°，与现场土的37.8°相近.

使用英国邓迪大学土工离心机Actidyn C67-2对根系模型A、B、Am、Bm进行倾覆试验，考虑冬季风暴条件下的工况，4次试验的土体均饱和. 模型试验盒的内尺寸为800 mm×500 mm×550 mm，如图10所示. 在容器底部铺10 mm厚的砾石层用于排水，铺设100 mm厚的HST95砂层，使用“落雨法”将试验土扬入容器内，控制其密实度为45%~50%（对应干密度为1.67±0.02 g/cm³）. 当土层深度达到模型根系底部时，将安装上树干的根系通过细线固定在图示位置，继续扬土，直到土体覆盖在根系顶部. 在扬土完成时，试验土的厚度为300 mm，再将土体自下而上进行充分饱和.

试验中，在达到20g超重力后，在距根系顶部60 mm处对树干施加水平位移以模拟风荷载，加载速率采用0.6 mm/min，以保证土体可以充分排水. 使用位移传感器记录树干的横向位移，并根据已有研究假设根系顶部为旋转点^[8]，按照下式计算弯矩：

(1)$ \theta = {\tan ^{ - 1}} \left( {\frac{{t \cdot {v_{\mathrm{H}}}}}{h}} \right), $

(2)$ M = Fh+0.5{m_{\text{t}}}NghL\sin \theta . $

式中：t为加载时间，v_H为加载速率，F为测量的水平力，h为横向位移和加载点与根顶部之间的距离（如图10所示），m_t和L分别为树干的质量和长度，N为重力加速的倍数，g为重力加速度. 所有关键变量（如土体密实度、根系结构、加载条件等）在试验中均受到严密控制.

根据相似律将离心试验结果转化为原型尺寸，其弯矩M和加载点转角$ \theta $的关系如图11所示. 当转角小于15°时，模型Am与A的弯矩-转角曲线重合较好，且弯矩都在约15°时达到峰值. 模型A的极限倾覆弯矩M_p = 8.4 kN·m，模型Am的M_p ≈ 9.6 kN·m，二者相差不大. 值得注意的是，简化模型Am的极限倾覆弯矩反而稍大一些. 简化前主根的直径随深度的增加逐渐减小，深度较浅处有效应力小而主根截面较粗，深度较大处有效应力大但主根截面较细，而采用平均直径的简化，增加了深度较大处主根的横截面积，提高了锚固效率. 将弯矩-转角图上0.5M_p处所对应点的割线斜率记为旋转刚度k_r，模型A和Am的k_r分别为0.87和0.95 kN·m/°相差5.7%，吻合较好. 当转角大于15°时，模型Am的树干和主根交界处由于应力集中而断裂，导致弯矩突然下降约65%，在转角大于15°后，简化后的根系和原根系的响应存在较大差异. 模型Am因应力集中而断裂，这从侧面揭示主根对极限倾覆弯矩的贡献占比为65%，对植被的抗倾覆起主导作用. 如图9 (b)所示，在0°~25° 转角区间内，模型Bm与B的弯矩-转角曲线大致重合，在转角约25°时弯矩达到峰值，B的极限倾覆弯矩M_p为11.2 kN·m，比Bm（M_p = 10.45 kN·m）高7.2%. 由于B的中央主根上竖向分支较多，在简化中都被忽略，且总体深度较浅，采用平均直径对锚固效果的提升不显著，简化后的模型Bm的极限倾覆弯矩有所降低. 由于 Bm的树干和主根的平均直径相差不大，未发生应力集中导致的破坏，其弯矩达到峰值后基本保持不变，与原根系的响应基本一致. 若不考虑应力集中导致根系的破坏，简化后模型的弯矩-转角曲线与原模型差异不大.

通过与Crook等^[20]的现场试验数据对比发现，在0~10°转角范围内，模型A、Am、B及Bm的弯矩-转角曲线与现场试验高度吻合，且离心试验测得的极限倾覆弯矩（9~12 kN·m）完全涵盖现场试验统计范围（7~20 kN·m），证实所提简化方法可有效捕捉根系倾覆的过程.

将4次试验中被破坏的单根在图1中用深色标出，可知，模型A、B的中央主根上的2级竖向根均易发生断裂，此特征与模型Bm相符. 由于模型Am的2级竖向根较深，在简化中被忽略，再加上主根的应力集中效应，只观测到主根的破坏. 此外，模型A、B的迎风区横向根均有1处被拔断，与现场试验中^[20]观察到的结果一致. 基于根系重构时所提出的假设，简化时忽略根的弯曲特性，进而减小了根-土相互作用时的接触面积，导致迎风区横向根均仅被拔出而未发生破坏. 同时，尽管只保留1、2、3级根，根系总体积下降了20%~30%，但对弯矩-转角图的影响不大，印证了假设的合理性.

由于离心试验的局限性，难以直接观察到根系在试验过程中的响应和断裂顺序. 结合Winkler地基梁法，采用非线性连接器替代土体单元对根进行约束，可以对根系的破坏进行量化，并解决有限元中大位移下网格严重变形导致无法计算的问题. 对横向外荷载下主根的响应进行初探，建立的模型如图12所示. 将主根离散成多个梁单元，其材料属性基于ABS塑料的抗拉和抗弯试验，参照文献[9]设置. 沿根轴线布置非线性连接器，连接器参数依据侧向受荷桩基的土体抗力p-y理论确定. 对于地表的楔形破坏模式，单位长度的极限侧阻力 p_st 采用Reese和Van Impe^[21]针对排水饱和砂土中竖向桩基开发的 p-y 理论按下式计算：

(3)$ \begin{split} {p_{{\text{st}}}} = &\sigma _{\mathrm{v}}^{'} \left[ \frac{{{K_0}X\tan \phi '{\text{ }}{\text{sin}}\; \beta }}{{\tan \left( {\beta - \phi '} \right) {\text{ cos}}\; \alpha }} + \frac{{\tan \beta \left( {d + X\tan \beta {\text{ }}\tan \alpha } \right)}}{{\tan \left( {\beta - \phi '} \right) }} \right. + \\& \Bigg. {K_0}X\tan \beta \left( {\tan \phi '{\text{ }}{\text{sin}}\; \beta - \tan \alpha } \right) - {K_{\mathrm{a}}}d \Bigg].\\[-5pt]\end{split} $

式中：$ {K_{\text{a}}} $和$ {K_0} $分别为主动和静止的侧向土压力系数，$ {K_0} = 1 - \sin \phi ' $，$ {K_{\text{a}}} = \left( {1 - \sin \phi '} \right)\left( {1+\sin \phi '} \right) $. 此外，$ \alpha = 0.5\phi ' $，$ \beta = {45^ \circ }+0.5\phi ' $，$ \phi ' $为土壤摩擦角. 在围压较大的深层区域，不会发生楔形破坏，其侧阻力为

(4)$ {p_{{\mathrm{sd}}}} = d \sigma _{\mathrm{v}}^{'} \left[ {{K_{\text{a}}}\left( ({{{\tan }}\; \beta)^8 - 1} \right)+{K_0}\tan \phi '({{\tan }}\; \beta )^4} \right]. $

最大侧阻力p_u为

(5)$ {p_{\mathrm{u}}} = \overline {{A_{\mathrm{s}}}} \min \;\left( {{p_{\mathrm{u}}},{p_{{\mathrm{sd}}}}} \right). $

式中：$ \overline {{A_{\mathrm{s}}}} $为单调加载的经验系数.

(6)$ \overline{{A}_{{\mathrm{s}}}}\approx \left\{\begin{array}{l}0.88+0.098\;8{\left(\dfrac{x}{d}-4.5\right)}^{2}\text{，}\text{ }\dfrac{x}{d}\leqslant 4.5\text{ }\text{；}\\0.88\text{，}\text{ }\dfrac{x}{d} \gt \mathrm{4.5.}\end{array}\right. $

由于需通过一定量级的根-土相对位移才能充分发挥轴向或横向阻力，根据Meijer等^[18]的研究，采用双曲正切曲线模拟该过程，连接器的横向土体阻力系数为

(7)$ {\xi _1} = \tanh \left( {\frac{{\Delta w}}{{{b_1}}}} \right). $

式中：$ \Delta w $为横向根-土相对位移，$ {b_1} $为控制连接器距离的参数.

(8)$ {b_1} = \frac{{{p_{\mathrm{u}}}}}{{{k_{{\mathrm{py}}}}{X_{{\mathrm{eq}}}}}}. $

式中：$ {k_{{\mathrm{py}}}} $为p-y曲线的初始刚度，中密砂为 24.4×10⁶ $ {X_{{\mathrm{eq}}}} $N·m⁻³；$ {X_{{\mathrm{eq}}}} $为深度. 连接器所产生的横向阻力为

(9)$ {q_1} = {\xi _1}{p_{\mathrm{u}}}. $

采用该方法对模型Am主根（图12(a)）的模拟结果如图13所示. 主根在转角33°时达到极限倾覆弯矩M_p = 3.9 kN·m，随后发生断裂，弯矩急剧下降. 对比离心试验中，模型Am在转角大于15°时因主根断裂出现4.3 kN·m的弯矩骤降，二者误差仅为9.3%，说明此方法可以较好地模拟主根的极限倾覆弯矩. 但是由于尚未考虑其他根的贡献，2条曲线刚度相差较大. 文献[22]将模型A（包含所有根）等效为直径53 mm的圆形桩基础（图12(b)），对此桩模型采用Winkler地基梁法模拟，结果见图13. 可以发现，初始刚度和离心试验吻合较好，这再次验证数值模拟方法的有效性. 但经此简化，由于各根都不再断裂，无法准确捕捉极限倾覆弯矩的大小. 仅对主根的响应进行初探，所采用的p-y曲线基于桩基理论，对于大直径的主根模拟效果较好，但对于其余竖向单根的适用性可能须考虑尺寸效应和柔性响应等因素的修正，后续工作可开展精细化建模进行更深一步的研究.

在实际工程应用中，基于植被根系的空间几何特征，建立如图14所示的梁单元简化模型，进而结合Winkler地基梁法开展根系数值模拟：将单根视为梁单元并离散为若干短根单元，沿轴线布置非线性弹簧，进行多自由度耦合约束，模拟复杂根-土相互作用. 该方法通过规避传统有限元网格划分，可以解决网格精度不足、计算效率低下及大变形模拟受限的问题，为根系锚固机理研究提供高效的数值模拟方法.

针对树木在风荷载作用下的稳定性进行深入研究，结合Winkler地基梁理论提出根系的重构方法，通过简化根系模型并结合超重力倾覆试验，对根系在风荷载下的响应进行研究，得到根系的弯矩转角图. 试验结果表明，在初始加载阶段，模型A、B与Am、Bm的弯矩-转角曲线高度吻合，表明简化重构方法能有效反映根系的初始力学响应；随着转角增大，模型Am因主根应力集中发生断裂，导致弯矩骤降约65%，凸显了主根在整体抗倾覆能力中的主导作用；而模型Bm则表现出与原始模型相近的渐进破坏特征。此外，两种重构模型均能再现迎风区横向根拔出的典型破坏模式，与现场观测结果一致. 在总旋转刚度方面，模型A、B与Am、Bm的总旋转刚度分别相差5.7%、15.3%. 这些结论说明通过这种根系重构方法得到的模型在预测树木在风荷载下的稳定性方面具有较高的准确性和可靠性.

将根系简化为水平根和竖直根规避了有限元中梁单元无法弯曲的问题，提高了计算效率，并可依托地基梁法捕捉根系在大变形条件下位移的响应，这为以后复杂根系的物理模拟和数值模拟奠定基础，也为土木工程领域中斜坡上树木稳定性的研究提供实践指导. 对于非均质土体（如碎石土边坡），由于土体颗粒级配和力学特性的空间变异性显著，尺寸效应显著，根系的简化方法、Winkler地基梁模型中的连接器布置方案、约束参数需要进一步修正. 特别是在含有大粒径碎石或显著分层特征的土体中，根系-土体相互作用机制可能发生本质变化. 未来的研究可以进一步探索不同风载条件、土体类型、土体饱和度、根系结构以及加载条件下的根系稳定性，以验证该方法在复杂工况下的适用性，以期为斜坡绿化和风灾预防提供更加全面的理论和技术支持.