无模型数据驱动的计算力学方法，用于断裂岩体的变形分析

《Journal of Rock Mechanics and Geotechnical Engineering》：Model-free data-driven computational mechanics for deformation analysis of fractured rock masses

　　冯帆丁|杨殿森|姜庆辉|刘耀来|刘东
武汉大学土木工程学院，中国武汉430072

摘要：在外部载荷作用下，断裂岩体的变形和应力分布仍然是关键的科学挑战，这主要是由于它们固有的异质性、不连续性和各向异性。为了解决这一问题，本研究提出了一种无模型数据驱动的计算力学方法，用于断裂岩体的变形分析，该方法能够更准确地捕捉它们在外部载荷作用下的力学响应。该方法通过直接结合实验数据（包括岩石应力-应变数据、节理法向应力-位移数据和切向应力-切向位移数据）来绕过传统的本构建模。通过最小化这些数据点与满足守恒方程（即平衡方程和几何方程）的数据点之间的距离来获得解。具体来说，断裂岩体被视为由岩石和节理组成的二元系统，其中岩石用固体元素表示，节理通过Goodman元素模拟以准确捕捉不连续变形行为。通过将数据驱动的计算力学方法与基于本构模型的参考解进行典型数值计算比较，验证了该方法的有效性和计算精度。此外，使用真实实验数据的验证表明，该方法在分析断裂岩体的非线性力学行为方面具有显著优势，为研究断裂岩体的变形响应提供了一种统一且实用的数值方法。

1. 引言
长期的复杂地质过程通常会在岩体内产生不连续性，如层理面、节理、断层和裂缝（Jing, 2003; Shang et al., 2016; Li et al., 2024）。这些结构表面的分布不均匀，具有不同的尺寸（例如长度）、形状和尺度（例如从小到大），并且相互交叉形成复杂的网络拓扑结构（如图1所示），从而赋予岩体异质性、不连续性和各向异性（Barton and Quadros, 2015; Tang et al., 2020; Wang et al., 2024）。断裂岩体的变形和破坏通常与这些复杂的结构特征密切相关，尤其是在存在大量不连续性时。结构表面对岩体力学性能的影响远大于完整岩体。断裂岩体的力学行为主要表现为在法向载荷下的正常变形和在切向载荷下的剪切滑移。准确计算断裂岩体在外部载荷作用下的变形和应力分布仍然是一个需要解决的关键科学挑战。

随着计算机科学技术的发展，数值模拟方法已成为描述和预测断裂岩体变形行为的有效工具。这些方法可以分为三类：（1）连续介质力学方法，包括有限元方法（FEM）（Yang et al., 2015; Jia and Zhou, 2023）、有限差分方法（FDM）（Yang et al., 2018; Eremin, 2020）和边界元方法（BEM）（Islam and Islam, 2016; Shi and Shen, 2025）；（2）不连续介质力学方法，如不连续变形分析（DDA）（Do and Wu, 2020; Xia and Yang, 2024）和离散元方法（DEM）（Regassa et al., 2018; Chen et al., 2024）；以及（3）连续-不连续耦合方法，如有限-离散元方法（FDEM）（Xu et al., 2023; Sharafisafa et al., 2024）和数值流形方法（NMM）（Wei et al., 2018; Wei et al., 2020）。每种方法都有其独特的优点和局限性，在选择和应用时必须仔细评估其适用性和计算能力。例如，连续介质力学方法在模拟大规模岩体的变形方面具有较高的计算效率。然而，它们通常需要通过数值实验确定代表性基本体积（REV）的大小，并依赖于特定的假设来推导等效模量参数。此外，这些方法无法准确模拟岩体的实际结构，忽略了节理之间以及节理与岩体之间的相互作用。相比之下，不连续介质力学方法可以更准确地模拟岩体的结构特征。然而，它们在模拟复杂的断裂网络方面面临挑战，包括高建模复杂性、低计算效率和收敛性差。耦合的连续-不连续介质力学方法结合了连续介质和不连续介质方法的数值模拟技术，实现了互补的优势。然而，连续介质和不连续介质域之间的力传递算法复杂，常常导致数值不稳定。

在用于断裂岩体的数值模拟方法中，有限元方法（FEM）因其能够高效分析各种载荷条件下的复杂变形响应而被广泛使用。在断裂岩体的有限元分析中，节理裂缝对岩体的变形行为和整体稳定性有显著影响。因此，节理裂缝是变形和应力演化分析中的关键因素。考虑这些效应的常见方法可以分为隐式方法和显式方法（Shin and Santamarina, 2019; Chen et al., 2022）。在隐式方法中，通过调整大体积岩体的强度和变形特性来将断裂网络的弱化和软化效应纳入计算。尽管这种方法计算效率高，但在高度无序的条件下难以准确捕捉断裂网络的不连续性以及实际的临界滑移表面。相比之下，显式方法将断裂表面视为岩体的离散元素，准确表示每个断裂表面的物理特性。这使得能够更准确地描述单个断裂的力学行为及其对岩体整体变形的影响（Li et al., 2022）。常见的显式断裂模拟方法包括薄层固体元素（Desai et al., 1984）、摩擦接触界面元素（Katona, 1983）和Goodman节理元素（Goodman et al., 1968）。其中，Goodman节理元素因其清晰的概念框架、参数确定的便利性和明确的物理意义而在断裂岩体模拟中得到广泛应用，有效模拟了不连续变形（Fischer and Diederichs, 2023）。

在使用显式有限元方法对断裂岩体进行变形分析时，采用适当的岩石和节理本构模型对于准确表示力学行为（如结构表面的滑动、闭合和开启）至关重要。如图2a所示，本构方程通常是通过拟合原材料的实验数据建立的。然后，它们与守恒方程（包括平衡微分方程和几何方程）同时求解，以寻找最优解（即它们的数学交集）。必须严格满足从基本物理原理推导出的守恒定律；然而，本构模型的选择具有一定的主观性。例如，使用线性或非线性本构关系可能导致不同的数值解，从而影响计算结果的准确性以及实验数据的原始性和可靠性。此外，开发和验证新的本构模型需要大量的智力努力和时间投入。在研究岩石本构模型时，学者们提出了多种基于岩体力学特性的模型，包括线性各向同性、非线性各向同性（Ren et al., 2023; Ru et al., 2023）、线性各向异性（Xu et al., 2017; Jin et al., 2024）和非线性各向异性（Kou et al., 2023; Li et al., 2024）本构模型，以更准确地描述不同岩体在复杂载荷条件下的力学响应。关于节理法向变形，研究人员开发了各种经验模型。经典模型包括Goodman双曲模型（Goodman, 1976）、Kulhawy双曲模型（Kulhawy, 1975）、Barton-Bandis（BB）模型（Bandis et al., 1983; Barton et al., 1985）、幂函数模型（Sun et al., 1985）和指数函数模型（Malama and Kulatilake, 2003）。最近，Rong et al.（2012）、Huang et al.（2018）和Li et al.（2022）提出了新的模型，以更好地捕捉节理的非线性变形特性。关于节理切向变形，经典模型包括Barton-Bandis（BB）模型（Bandis et al., 1983; Barton et al., 1985）、Saeb-Amadei节理模型（Saeb and Amadei, 1992）和Seidel-Haberfield模型（Seidel and Haberfield, 2002）。最近，Bahaaddini et al.（2015）和Feng et al.（2024）提出了新的模型，以更好地捕捉节理的非线性剪切变形特性。

最近，Kirchdoerfer和Ortiz（2016）提出的数据驱动计算力学（DDCM）方法为解决这些问题提供了新的视角。如图2b所示，DDCM试图通过直接利用大规模实验数据、守恒定律和边界条件来绕过传统的本构模型开发过程。它将耦合的平衡微分方程、几何方程和本构方程的数值解过程转化为一个优化问题。通过迭代求解，该方法从数据库中识别每个数值积分点的最佳材料数据对，从而显著减少了材料建模的经验性、易出错性和不确定性，同时保留了实验信息。由于其无模型和高保真度的特点，数据驱动方法受到了越来越多的关注，并已逐步应用于动力学（Kirchdoerfer and Ortiz, 2018; Zheng et al., 2020; Garcia-Suarez et al., 2023）、大变形（Nguyen and Keip, 2018; Platzer et al., 2021）和弹塑性问题（Eggersmann et al., 2019; Ciftci and Hackl, 2022; Dandin et al., 2024）。在数据驱动计算力学的基础框架内，岩石应力-应变数据库可以替代传统的材料本构模型。在此基础上，作者提出了一种用于层状岩石变形分析的数据驱动方法（Feng et al., 2025），其核心过程是通过实验构建岩石应力-应变数据库，并直接利用离散应力-应变数据点来计算特定载荷条件下各向异性和非线性层状岩体的变形响应。然而，传统的DDCM框架受到固体元素离散化策略的限制，只能处理连续介质力学问题或将断裂网络的弱化和软化效应隐式地耦合到计算模型中。这使得难以实现断裂岩体不连续变形行为的显式和高精度建模。为了解决这一限制，本研究进一步改进了DDCM框架，将断裂岩体建模为由岩石和节理组成的二元系统，其中完整岩石通过固体元素离散化表示，节理使用Goodman元素准确模拟。该方法消除了对岩石和节理预定义本构关系的需求，直接利用岩石应力-应变数据库、节理法向数据库（法向应力-法向位移关系）和切向数据库（切向应力-切向位移关系）来计算解，从而能够准确模拟断裂岩体的不连续变形特性。在实际应用中，岩石应力-应变数据库可以使用数字图像相关（DIC）（Musial et al., 2019）技术和数据驱动识别（DDI）（Leygue et al., 2018）方法获得，而节理法向和切向数据库可以从节理法向加载试验（Malama and Kulatilake, 2003）和节理切向直剪试验（Asadollahi et al., 2010; Li et al., 2022）中得出。这种方法不仅能够更准确地模拟断裂岩体的复杂变形行为，还绕过了选择和参数化传统本构模型所带来的挑战，为断裂岩体的力学分析提供了一种更高效和可靠的方法。

本研究首先回顾了经典的断裂岩体显式有限元方法分析框架。在此基础上，提出了一种无模型的数据驱动计算力学方法用于断裂岩体的变形分析，并详细解释了其理论基础和求解过程。随后，通过对几个代表性数值示例的模拟分析，将数据驱动解与基于传统本构模型的参考解进行比较，从而验证了所提方法的有效性和鲁棒性。此外，通过结合真实实验数据，进一步突出了该方法在实际工程应用中的优势及其强大的非线性计算能力。基于这些结果，对方法的优点、局限性和未来研究的方向进行了批判性评估。总之，本研究的主要研究结果已总结如下。2. 无模型数据驱动计算力学的理论基础2.1. 裂缝岩体的连接有限元方法框架在经典的裂缝岩体连接有限元方法（JFEM）分析框架中，裂缝岩体被视为岩石和节理的复合系统，其中岩石使用固体元素离散化，而节理、裂缝和其他不连续性则使用Goodman节理元素进行建模（如图3所示）。在求解过程中，首先推导出固体元素和节理元素的刚度矩阵，然后将它们组合成全局刚度矩阵。接下来应用边界条件来求解未知的自由度，随后计算元素内的应变和应力分布，从而便于对具有不连续结构（如节理和裂缝）的岩体进行数值模拟和分析。下载：下载高分辨率图像（423KB）下载：下载全尺寸图像图3. 裂缝岩体和JFEM原理示意图。Goodman节理元素由Goodman等人（1968年）提出，是一种四节点、零厚度的元素，其特征是两条平行直线。它为节理元素的研究奠定了重要基础。Goodman节理元素引入了沿剪切滑动面的切向刚度ks和沿法向滑动面的法向刚度kn。假设接触面上的法向应力和切向应力与相应的相对法向位移和切向位移无关。应力与相对位移之间的关系表示为：(1)其中σs是切向应力，MPa；σn是法向应力，MPa；ks是切向刚度，MPa/mm；kn是法向刚度，MPa/mm；ws是切向位移，mm；wn是法向位移，mm。Goodman节理元素准确模拟了节理与岩石之间的不连续变形，包括接触面的滑动、张开和闭合。该方法概念清晰，法向刚度kn和切向刚度ks是主要参数，具有明确的物理意义，并且可以通过对节理进行正常加载试验和切向直剪试验获得。预测裂缝岩体中应力和变形的准确性在很大程度上取决于法向刚度kn和切向刚度ks的精度。目前，刚度系数的测量或估计面临若干挑战，包括实验条件的复杂性、节理表面粗糙度的影响以及与尺寸效应相关的不确定性。此外，节理特性测试表明，法向应力与法向变形以及切向应力与切向变形之间的关系通常表现出显著的非线性行为（Li等人，2022年）。2.2. 无模型数据驱动计算力学的裂缝岩体框架如2.1节所述，裂缝岩体的JFEM需要预先定义的本构模型和固体元素及节理元素的力学参数来准确捕捉其力学行为。然而，这些本构模型通常是通过实验数据的经验拟合得到的，这可能导致原始信息的丢失，从而影响数值模拟的准确性和可靠性。为了克服传统本构模型的局限性，Kirchdoerfer和Ortiz（2016年）引入了一种新方法，该方法利用应力-应变数据库直接驱动力学模拟。该方法将耦合材料的本构定律（CL）、平衡微分定律（EL）和运动定律（KL）的求解过程转化为一个优化问题。优化目标是确定最接近的数据分配方案（实验材料数据点）、允许变形的应变场（一组受几何方程约束的可行解）以及允许静态的应力场（一组受平衡方程约束的可行解）。数据驱动计算力学的计算框架如图2b所示。优化目标是在材料数据和守恒定律之间获得最佳近似，而无需明确定义应力与应变之间的本构关系。传统的DDCM框架受到固体元素离散化方法的限制，只能处理连续介质力学问题，或者需要隐式耦合来近似裂缝网络的弱化和软化效应。在本研究中，DDCM框架被扩展为显式构建块-节理二进制离散化系统，从而能够精确模拟裂缝岩体的不连续变形行为。这种方法消除了对岩石和节理预先定义本构关系的需求，而是利用实验数据（如岩石应力-应变数据库和节理法向及切向应力-位移数据库）直接驱动力学模拟。这种方法不仅消除了基于经验假设的不确定性，还显著提高了数值模拟的准确性和可靠性，为裂缝岩体力学行为的数值分析提供了创新的研究方向。在这种方法中，计算不再依赖于显式的材料本构关系。相反，它通过引入类似于应变能量密度和应变残余能量密度的度量来定义空间距离。对于固体元素，空间距离定义为满足平衡微分方程和几何协调方程的力学允许解与材料数据库中的实验数据对之间的度量距离。对于节理元素，空间距离类似地定义为满足相同方程的力学允许解与实验数据对之间的度量距离。(2)(3)通过引入适当的权重系数，可以将每个固体元素和节理元素的空间距离组合起来，构建一个表示全局优化目标的惩罚函数G：(4)相应的约束包括：(5)(6)(7)其中G是全局优化惩罚函数； 是固体元素积分点索引； 是节理元素积分点索引； 是固体元素积分点的积分权重（指元素面积）； 是节理元素积分点的积分权重（指节理元素长度）；Ce和Cje分别是固体元素和节理元素的平衡因子矩阵，用于归一化独立分量之间的尺寸差异，并防止某些分量在计算过程中被过度放大或忽略； 是固体元素积分点的应变变换矩阵； 是节理元素积分点在全局坐标系中的变换矩阵； 是节理元素积分点在局部坐标系中的变换矩阵；T是坐标变换矩阵； 是满足位移边界条件的可行位移向量列矩阵； 是外力载荷向量列矩阵； 是固体元素的应变-应力数据库； 是节理元素的位移-应力数据库。在二维数值分析中，固体元素的平衡因子矩阵Ce表示为一个3×3的对称正定矩阵；节理元素的平衡因子矩阵Cje表示为一个2×2的对称正定矩阵。通过引入拉格朗日乘数，我们构建了增广目标函数来求解上述优化问题：(8)节点位移向量、固体元素应力分量向量、节理元素应力分量向量和拉格朗日乘数向量可以如下推导：(9)可以从方程(9)中得出以下两个线性控制方程：(10)(11)在求解过程中，必须考虑这两个线性方程的边界条件：节点位移向量在受限位移自由度处固定，而拉格朗日乘数向量在这些约束处为零。然后可以处理从上述方程中得到的系数 和 ，以获得相应的允许变形的应变场 和 ，以及相应的允许静态的应力场 和。(12)(13)在每个数值积分点获得的 和 值通常是独立的，可能不符合宏观材料本构关系。因此，必须在相应的材料数据库中搜索与当前结果最接近的数据点，然后将它们重新分配为最优解 和 。这个优化问题可以在每个数值积分点独立分解和求解。当所有积分点都从数据库中检索到最近的数据点时，全局优化目标最小化。具体来说，所有积分点必须满足以下约束条件。(14)数据驱动计算力学中的优化问题可以简化为一个整数二次规划问题。在本研究中，我们采用可分离凸优化理论中的交替最小化算法来解决问题。该算法将原始优化问题分解为两个耦合的最小化子问题，分别称为子问题PM和子问题PC。子问题PM可以描述如下（对应于方程(10)、(11)、(12)、(13)：从实验获得的材料数据集 和 开始，目标是求解在满足静态平衡和变形兼容性的同时最小化全局目标函数的力学可行解 和 。子问题PC可以描述如下（对应于方程(14)）：从满足静态平衡和变形兼容性的力学可行解 和 开始，任务是寻找最小化目标函数G的最优实验数据对 和 。通过交替求解子问题PM和子问题PC，当满足方程(15)的条件时，迭代过程终止，表明算法已收敛到最优解。(15)图4a展示了数据驱动交替最小化迭代过程。左圆圈标记为EL/KL，表示力学状态（Ms），由满足平衡和运动方程的材料数据点 和 组成。右圆圈标记为CL，表示本构状态（Cs），由从数据库中选择的材料数据点 和 组成。数据驱动算法的优化目标是确定每个积分点的解，以最小化力学状态（Ms）和本构状态（Cs）之间的差异，即两个圆圈之间的最接近解。此外，考虑到数据库采样过程中可能存在噪声，并且认识到平衡和运动方程是确定性的，我们选择满足静态平衡和变形兼容性的力学状态作为可行的应力-应变解，表示为 和 。下载：下载高分辨率图像（489KB）下载：下载全尺寸图像图4. 数据驱动算法。(a) 交替最小化迭代示意图，(b) 算法流程图。图4b显示了数据驱动算法的流程图，迭代步骤如下：步骤1：从数据库中随机为每个单元分配 ；步骤2：使用 方程(10)、(11)和边值问题来找到ukl(k)和λ(k)；步骤3：从方程(12)、(13)中找到 ；步骤4：通过方程(14)找到 ；步骤5：重复步骤2-4直到 。3. 数值示例本节通过三个典型的数值示例系统验证了所提出的无模型数据驱动计算力学在裂缝岩体变形分析中的应用：一个贯穿的裂缝岩体、一个随机裂缝岩体和一个工程规模的裂缝岩体。需要注意的是，该方法所需的材料数据库包括三个关键组成部分：岩石应力-应变数据库、节理法向应力-法向位移数据库和节理切向应力-切向位移数据库。在实际应用中，这些数据库可以通过相应的实验获得。然而，为了验证该方法的计算准确性，本节基于岩石和节理的本构模型构建了一个合成数据库，详细构建方法在3.1节中有描述。如果数据驱动解收敛到参考解，则验证了数据驱动计算力学的有效性。使用相同的元素网格在每个元素上将数据驱动解与参考解进行视觉比较。为了最小化数据库中每个积分点的初始分配对实验结果的影响，在本研究中呈现的数值实验中应用了相同的初始值。本节还研究了不同规模数据库对计算准确性、效率和收敛性的影响。数据驱动计算程序是在Visual Studio 19.0平台上用C++开发的，而基于本构模型的参考解则是使用JFEM计算的。3.1.数据库的构建

(1) 关于岩石应力-应变数据库的构建方法：基于本构矩阵和适当的应变范围（预先估计以确保数据库覆盖研究问题的解空间），生成了一个用于数据驱动计算的合成岩石数据库。首先，针对具体的研究问题，确保应变分量在指定范围内均匀分布，例如 εx∈[0, 1×10-3]；εy∈[-2×10-3, 0]；γxy∈[-1×10-4, 1×10-4]。接下来，将应变张量代入刚度矩阵中，计算出相应的应力张量，从而构建应力-应变数据库。对于二维问题，岩石数据库包括三个独立的应变分量及其对应的三个应力分量。通过调整每个应变范围内的样本数量，可以生成不同大小的应力-应变数据库。图5展示了不同规模的数据库在应变和应力空间中的分布（岩石数据库规模分别为Mr=10^3、30^3和50^3），为数据驱动计算提供了必要的数据基础。

下载：下载高分辨率图像（1MB）
下载：下载全尺寸图像
图5. 不同规模的岩石数据库。(a) 应变数据（Mr = 10^3），(b) 应变数据（Mr = 30^3），(c) 应变数据（Mr = 50^3），(d) 应力数据（Mr = 10^3），(e) 应力数据（Mr = 30^3），以及 (f) 应力数据（Mr = 50^3）。

(2) 关于节理法向和切向应力-位移数据库的构建方法：基于节理的法向和切向本构模型（如方程(1)所示），并考虑适当的法向位移和切向位移范围（例如 wn∈[-3×10-4, 3×10-4]；ws∈[-1.5×10-3, 1.5×10-3]），计算出相应的应力分量，从而获得节理的合成法向应力-法向位移和切向应力-切向位移数据库。对于二维问题，节理法向数据库包括法向位移wn及其对应的法向应力σn，而节理切向数据库包含切向位移ws及其对应的切向应力σs。图6展示了不同规模的节理法向和切向数据库在力学空间中的分布（节理法向数据库规模Mn和节理切向数据库规模Ms分别设置为20、100和2000），为数据驱动计算提供了必要的数据支持。

下载：下载高分辨率图像（665KB）
下载：下载全尺寸图像
图6. 不同规模的节理数据库。(a) 不同规模的节理法向应力-法向位移数据库，以及 (b) 不同规模的节理切向应力-切向位移数据库。

3.2. 穿透性断裂岩体

假设在断裂岩体中存在两组穿透性节理面。第一组节理面与水平方向成60°角，平行间距为1米；第二组节理面与水平方向成30°角，平行间距为3米。断裂岩体的尺寸为长度10米，高度10米，如图7a中的数值模型所示。岩石采用各向同性线性弹性模型表示，弹性模量为E = 10.08 GPa，泊松比为v = 0.26。节理采用线性弹性模型表示，法向刚度kn = 50 GPa/m，切向刚度ks = 5 GPa/m。y方向位移在底部边界受到约束，而x方向和y方向的位移在左下边界受到约束。在顶部边界施加均匀压缩应力5 MPa。将DDCM解与JFEM和等效有限元方法（EFEM）进行比较，其中等效刚度张量是通过解析解获得的（详见附录A）。DDCM和JFEM模型使用相同的网格，如图7b所示，包含1887个固体元素和314个节理元素（岩石由三节点三角形元素表示，节理由Goodman节理元素建模），而EFEM使用的网格包含638个固体元素，如图7c所示。

下载：下载高分辨率图像（632KB）
下载：下载全尺寸图像
图7. 穿透性断裂岩体的数值模型和网格配置。(a) 模型示意图，(b) JFEM和DDCM的网格，以及 (c) EFEM的网格。

在断裂岩体变形的EFEM分析中，将附录A中计算的等效刚度矩阵纳入计算模型。在JFEM分析中，应用了岩石和节理的本构模型。在DDCM方法中，岩石数据库（Mr = 503）、节理法向数据库（Mn = 2000）和节理切向数据库（Ms = 2000）替代了传统的岩石和节理本构模型。具体来说，岩石数据库中的应变范围如下：εx∈[0, 1×10-3]；εy∈[-2×10-3, 0]；γxy∈[-3×10-4, 3×10-4]；节理法向数据库中的法向位移范围为：wn∈[-3×10-4, 3×10-4]；节理切向数据库中的切向位移范围为：ws∈[-1.5×10-3, 1.5×10-3]。需要注意的是，为了确保数据库充分准确地反映岩体的变形特性，必须提前评估所解决问题的应力和应变水平，从而确保数据库包含给定情况下的解。图8显示了在5 MPa正常压力下，使用JFEM、DDCM和EFEM方法获得的水平和垂直位移的位移等值线图。比较表明，三种方法生成的位移等值线图一致，在x方向和y方向的位移上显示出相似的趋势。这表明这三种方法都能有效捕捉断裂岩体在宏观尺度上的变形行为。然而，在局部尺度上，DDCM和JFEM方法更好地捕捉了岩体内部断裂的不连续性，提供了更准确的变形特性描述。此外，由于DDCM方法依赖于真实的实验数据，没有引入额外的假设或近似，其结果更加准确可靠。

下载：下载高分辨率图像（1MB）
下载：下载全尺寸图像
图8. 穿透性断裂岩体的位移场比较。(a) JFEM-DIS-X，(b) DDCM-DIS-X，(c) EFEM-DIS-X，(d) JFEM-DIS-Y，(e) DDCM-DIS-Y，以及 (f) EFEM-DIS-Y。

在DDCM方法中，为了评估数据库对解空间的覆盖程度并确保结果的准确性和稳定性，本研究通过可视化分析了数据库中数据驱动解的分布。引入了最小相对距离（MRD）作为定量指标。MRD的定义如下：首先将每个解相对于相应数据库边界归一化到[0, 1]范围内；然后计算每个维度到边界（0或1）的最小距离；最后取所有维度中的最小值作为MRD。较大的MRD表示解位于数据库内部，反映了足够的覆盖范围。假设MRD低于某个阈值（例如0.05），则认为解位于数据库边界或外推区域，会立即触发警报。该区域被标记为“数据稀疏区”，并记录其应力/应变或应力/位移坐标。建议用户谨慎解释结果。图9a和b显示了穿透性断裂岩体中岩石数据驱动解在应变和应力数据库中的分布，图9c展示了相应的应变和应力MRD值。图9d和e分别展示了节理法向应力-位移和切向应力-位移数据库中数据驱动解的分布，图9f显示了相应的法向和剪切MRD值。如图所示，所有解点都位于数据库内部，MRD分布进一步确认没有点低于规定的阈值。这表明对于当前情况，数据库充分覆盖了解空间，计算结果非常可靠和稳定。

下载：下载高分辨率图像（675KB）
下载：下载全尺寸图像
图9. 穿透性断裂岩体的数据驱动解和MRD分布。(a) 岩石应变解，(b) 岩石应力解，(c) 岩石MRD分布，(d) 节理法向解，(e) 节理切向解，以及 (f) 节理MRD分布。

为了系统评估不同规模数据库对DDCM计算精度的影响，本研究采用均方根误差（RMS）（Kirchdoerfer和Ortiz，2016）作为评估指标，将DDCM解与JFEM参考解进行比较。此外，本研究还探讨了不同规模数据库的计算效率和收敛行为。值得注意的是，节理法向和切向数据库的规模固定为Mn = Ms = 2000，随后对不同规模数据库的引用特指岩石数据库，因为节理数据库规模的变化对计算结果的影响相对较小。图10a展示了不同规模数据库的应变和应力RMS误差以及所需的迭代次数。观察到，随着数据库规模的增加，应变和应力的RMS误差逐渐减小。这主要是由于数据库的扩展，扩大了围绕准确解的可行解集，从而引入了新的全局最优解。然而，随着数据库规模的增加，数据驱动迭代所需的时间显著增加。计算成本主要集中在子问题PC中，该问题涉及在每个积分点全局搜索最近的力学和材料状态。例如，当数据库规模从10^3增加到50^3时，计算时间从大约26秒增加到373秒。图10b显示了不同规模数据库所需的迭代次数。从图中可以看出，随着迭代次数的增加，全局目标惩罚函数G逐渐减小，并最终收敛到唯一的稳定解。特别是对于较大的数据库规模，惩罚函数的最终收敛值更小，表明数据驱动解的精度更高。此外，随着数据库规模的增加，迭代次数通常也会增加。

下载：下载高分辨率图像（392KB）
下载：下载全尺寸图像
图10. 不同规模数据库的RMS误差、计算时间和迭代次数。(a) RMS误差和计算时间，以及 (b) 迭代次数和收敛性（注：'pts'是'points'的缩写）。

3.3. 随机断裂岩体

根据Baecher节理网络模型（Baecher等人，1977）生成了两组节理。该模型考虑了节理的随机分布和几何特性，包括节理长度、方向和密度等参数（具体参数列在表1中）。随机断裂岩体的数值模型如图11a所示。岩石采用各向同性线性弹性模型表示，弹性模量为E = 10.08 GPa，泊松比为v = 0.26。节理采用线性弹性模型表示，法向刚度kn = 50 GPa/m，切向刚度ks = 5 GPa/m。y方向的位移在底部边界受到约束，而x方向和y方向的位移在左下边界受到约束。在顶部边界施加均匀压缩应力5 MPa。将DDCM解与JFEM进行比较，如图11b所示。网格包含1903个固体元素和160个节理元素。

表1. 节理参数的概率分布。
| 节理 | 节理强度 | 节理长度 | 倾斜/倾斜方向 | 节理数量/面积 | 平均值（m） | 标准差 |
| ---- | -------- | -------- | ------------ | -------- | -------- |
| J1 | 0.1 | Normal | 41 | 150 | 10 | 0.1 |
| J2 | 0.15 | Normal | 30 | 7 | 0.7 |

下载：下载高分辨率图像（498KB）
下载：下载全尺寸图像
图11. 随机断裂岩体的数值模型和网格配置。(a) 模型示意图，以及 (b) JFEM和DDCM的网格。

在DDCM方法中，使用岩石数据库（Mr = 503）、节理法向数据库（Mn = 2000）和节理切向数据库（Ms = 2000）替代岩石和节理的本构模型。岩石数据库的应变范围如下：εx∈[-5×10-4, 1×10-3]；εy∈[-2×10-3, 0]；γxy∈[-1×10-3, 1×10-3]。节理法向数据库中的法向位移范围为wn∈[-3×10-4, 3×10-4]，节理切向数据库中的切向位移范围为ws∈[-1.5×10-3, 1.5×10-3]。图12显示了在5 MPa正常压力下，使用JFEM和DDCM得到的水平和垂直位移等值线图。两种方法的比较表明，位移等值线图总体上是一致的，这证实了数据驱动计算方法在宏观尺度上的有效性和可靠性。在局部尺度上，由于裂缝对局部应力和位移场有显著影响，通常可以在裂缝处观察到位移不连续性。然而，由于众多裂缝的累积效应，在10米模型中可以观察到相对均匀的变形梯度趋势。下载：下载高分辨率图像（851KB）下载：下载全尺寸图像图12. 通过裂缝岩石体的位移场比较。(a) JFEM-DIS-X，(b) DDCM-DIS-X，(c) JFEM-DIS-Y，(d) DDCM-DIS-Y。图13a和b显示了随机裂缝岩石体中应变和应力数据库空间内的数据驱动解点分布，而图13c展示了相应的MRD值。图13d和e分别展示了接缝数据驱动解点在法向和切向数据库空间内的分布，图13f提供了相应的MRD值。结果表明，解点也主要集中在数据库的内部区域，MRD值超过了规定的阈值，从而避免了潜在的外推风险。这表明，即使在更复杂的裂缝分布下，构建的数据库仍然能够充分覆盖计算所需的状态空间，确保了结果的可靠性。下载：下载高分辨率图像（715KB）下载：下载全尺寸图像图13. 随机裂缝岩石体的数据驱动解和MRD分布。(a) 岩石应变解，(b) 岩石应力解，(c) 岩石MRD分布，(d) 接缝法向解，(e) 接缝切向解，(f) 接缝MRD分布。图14a展示了不同尺度数据库的应变和应力的RMS误差以及所需的迭代次数。随着数据库规模的增加，应变和应力的RMS误差逐渐减小，而数据驱动迭代解所需的时间显著增加。图14b显示了不同尺度数据库所需的迭代次数。可以观察到，随着迭代次数的增加，全局目标惩罚函数G逐渐减小并最终收敛到一个唯一的稳定解；此外，数据库规模越大，最终收敛的G值越小。下载：下载高分辨率图像（410KB）下载：下载全尺寸图像图14. 不同尺度随机裂缝岩石体的RMS误差、计算时间和迭代次数。(a) RMS误差和计算时间，(b) 迭代次数和收敛性（注：'pts'是'points'的缩写）。3.4. 工程尺度裂缝岩石体工程经验表明，裂缝在岩石边坡的表层最为发育，那里明显的卸载和张开现象显著影响了岩石体的力学性质。相比之下，内部岩石体相对完整，裂缝结合良好，力学性质良好。因此，可以采用一种简化的建模方法（Wang等人，2024年），该方法采用分层划分：在表层裂缝区构建详细的离散裂缝网络，而在更深层，小裂缝被表示为等效的均匀单元（Jiang等人，2017年），只有足够大的裂缝被明确建模。图15a展示了Lawa水电站右岸进水口D-D横截面的二维数值分析模型和开挖体。岩石使用各向同性线性弹性模型表示，其力学参数如表2所示。接缝使用线性弹性模型描述，法向刚度kn = 50 GPa/m，切向刚度ks = 5 GPa/m。模型的水平范围为625米，垂直范围为520米，x方向的左右边界以及x和y方向的底部边界都施加了约束。根据风化程度，将边坡岩石体分为三类：V类、IV类和III类。裂缝网络是根据地质调查得出的各种岩石体中结构表面的几何和空间分布特征生成的（见表3和4）。下载：下载高分辨率图像（3MB）下载：下载全尺寸图像图15. 工程岩石体的数值模型和网格配置。(a) 模型示意图，(b) JFEM的网格，(c) DDCM的网格。表2. 岩石力学参数值。岩土层自然密度，ρ (g/cm3)弹性模量，E (GPa)泊松比，v高度风化层（V）2.31.50.35中度风化层（IV）2.530.32轻微风化层（III）2.810.080.26表3. 并联确定性接缝模型接缝参数的概率分布。组裂缝集合方向（°）长度（m）间距（m）持久性V, IVJ1702050.6J2401050.5表4. Baecher接缝模型接缝参数的概率分布。组定义方法接缝强度接缝长度倾角/倾角方向接缝/面积分布平均值（m）平均倾角（°）平均倾角方向（°）IIIFisher (K=1)0.02指数3036135为了验证数据驱动计算方法的可行性，分别使用JFEM和DDCM对模型在自重和开挖条件下的变形进行了分析。裂缝岩石体的JFEM模型由67,346个固体元素和20,410个接缝元素组成，如图15b所示。在DDCM计算模型中，表现出显著力学变化的V类和IV类岩石体被视为由岩石和接缝组成的二元结构，计算使用岩石和接缝数据库进行。表现出更稳定力学性质的III类岩石体被视为由岩石和接缝组成的复合材料，基于等效连续介质力学的概念。使用代表性基本体积（REV）的概念创建了等效裂缝岩石体的数据库（详见3.4.1确定REV大小和3.4.2数据库生成部分）。最终的DDCM模型包括10,282个固体元素和1767个接缝元素，如图15c所示。3.4.1. 确定REV大小在为III类岩石体构建裂缝岩石体网络模型后，进行了数值实验以确定裂缝岩石体的REV大小。需要注意的是，REV研究旨在确定后续数值测试中使用的表征单元的适当大小，以生成等效的裂缝岩石体数据库。通过分析以同一点为中心的不同尺寸的岩石体，研究了裂缝岩石体样本的尺寸效应。随着计算岩石体尺寸的增加，当等效参数在所有方向上稳定时，认为达到了代表性元素的大小（Yang等人，2014年）。计算示意图和等效变形模量如图16所示。该图展示了在0°、30°、60°、90°、120°和150°方向上，边长为10米、20米、30米、40米、50米、60米、70米和80米的岩石体在轴向加载下的等效变形模量。玫瑰图表明，当研究区域的尺寸达到60米×60米或更大时，轴向等效变形模量在所有加载方向上稳定，表明已经达到了代表性元素的大小。因此，在以下部分中，使用60米×60米的代表性元素大小来生成数据库。下载：下载高分辨率图像（722KB）下载：下载全尺寸图像图16. 确定III类岩石体REV大小的方法。(a) 不同尺寸和方向的计算条件示意图，(b) 不同方向下岩石变形模量与岩石大小之间的关系。3.4.2. 数据库生成通过在REV上施加均质化位移/应变边界条件，获得了III类岩石体的等效应力-应变数据库。具体来说，在裂缝岩石体的REV尺度上，使用岩石和接缝数据库在均质化应变边界条件下进行计算，其中记录了每个加载步骤施加的应变和相应的计算等效应力作为数据点。这样，构建了REV尺度上的裂缝岩石体应力-应变数据库，随后用于工程尺度上的数值模拟。通过对REV重复施加不同的均质化应变边界条件，共生成了203个应力-应变数据点。这些数据点在每个分量方向上的应变分布如下：εx∈[-1×10-3, 1×10-3]，等间隔加载(Δεx=1×10-?)；εy∈[-2×10-3, 5×10-?]，等间隔加载(Δεy=1.25×10-?)；γxy∈[-2×10-3, 2×10-3]，等间隔加载(Δγxy=2×10-?)。V类岩石体（Mr = 203）、IV类岩石体（Mr = 203）、接缝法向数据库（Mn = 2000）和接缝切向数据库（Ms = 2000）的数据库是根据以下构成的方程结合适当的应变或位移范围生成的：εx∈[-1×10-3, 1×10-3]；εy∈[-2×10-3, 5×10-?]；γxy∈[-2×10-3, 2×10-3]；wn∈[-1×10-3, 1×10-3]；ws∈[-3×10-3, 3×10-3]。3.4.3. 计算解使用生成的工程岩石体数据库，通过数据驱动算法计算岩石体在自重和开挖作用下的变形场。在自重作用下，所有元素的应力和应变状态被随机初始化，并通过迭代过程逐渐收敛。在开挖过程中，开挖元素的应力在初始化时被设置为零，而它们的应变状态从数据库中随机采样。在PM映射中，开挖元素的平衡因子矩阵乘以一个软化因子（通常取为10-?），允许应变自由调整，同时保持应力为零。在PC映射中，移除了开挖元素的整数约束，使得应变可以直接继承，而应力保持为零。最后，通过交替迭代程序减去重力作用下的变形和开挖作用下的变形，可以获得由开挖引起的净变形。在数据驱动的计算解中，通过从数据库中的应力-应变数据构建冗余方程并使用最小二乘法求解来确定平衡因子矩阵C。(16)其中Ne和Nje分别是岩石和接缝数据库的索引；edata和jedata分别是岩石和接缝数据库中的数据点总数；σcl和εcl是岩石的应力-应变数据对；σcl和wcl是接缝的应力-位移数据对；Ce和Cje分别是岩石和接缝的平衡因子矩阵。图17和18显示了在自重和开挖条件下JFEM和DDCM的位移等值线图。结果表明，两种方法在大小和分布方面产生了高度一致的变形模式。在自重作用下，岩石体通常表现出向下运动，最大垂直位移发生在坡顶附近，因为重力应力集中，而最大水平位移发生在坡脚附近。在开挖过程中，由于应力的重新分布，会发生垂直回弹。同时，水平位移在开挖线上方向外，在开挖线下方向内，反映了岩石体内的局部应力释放和支持效应。在V类和IV类岩石体中，接缝和裂缝强烈影响变形和应力分布，使得传统的连续介质模型不够准确。DDCM使用岩石和接缝数据库而不是本构方程，准确地捕捉了这些效应。在III类岩石体中，由于不连续性较小，DDCM采用基于等效应力-应变数据库的等效连续介质方法，有效地再现了宏观力学行为。这种方法简化了模型，降低了计算成本，并保持了准确性，展示了其在复杂地质条件下的工程设计潜力。下载：下载高分辨率图像（710KB）下载：下载全尺寸图像图17. 自重条件下的位移场分布。(a) JFEM-DIS-X，(b) DDCM-DIS-X，(c) JFEM-DIS-Y，(d) DDCM-DIS-Y。下载：下载高分辨率图像（473KB）下载：下载全尺寸图像图18. 开挖条件下的位移场分布。(a) JFEM-DIS-X，(b) DDCM-DIS-X，(c) JFEM-DIS-Y，(d) DDCM-DIS-Y。4. 非线性问题的实际应用本节使用数据驱动计算力学来分析接缝的法向和切向变形，利用来自接缝法向加载试验和接缝切向直剪试验的实际数据。结果与传统的接缝本构模型进行了比较，全面突出了无模型数据驱动方法在分析非线性接缝变形方面的显著优势。上述介绍的算法是距离最小化数据驱动计算力学（DMDDCM），它使用空间距离作为度量标准，并采用每个单元的单个数据点信息来提取材料属性。材料数据库中的样本变化、测量误差、环境波动和人为操作等因素可能会引入数据噪声，从而可能影响传统DMDDCM方法的计算精度。为了解决这个问题，本研究提出了局部凸性数据驱动计算力学（LCDDCM）方法（He和Chen，2020年）。该方法利用局部凸性构建的概念有效减轻数据噪声的影响，从而提高数据驱动计算中的鲁棒性和收敛稳定性（详见附录B）。

4.1. 节缝正交加载试验示例
参考文献（Malama和Kulatilake，2003年）对亚利桑那州地表露头中的花岗岩-闪长岩样本进行了正交闭合试验。使用了每个包含一个垂直于加载方向的裂缝的矩形块。样品（Grd-1、Grd-2、Grd-3）的裂缝表面变化如图19所示。使用CT-500伺服控制系统施加了25-50 MPa的正交压缩载荷，并用LVDT（精度±1 μm）测量了正交位移。通过从总变形中减去完整样品的变形来获得裂缝的正交闭合。实验结果（图20）显示了在正交加载条件下的典型非线性行为。

下载：下载高分辨率图像（529KB）
下载：下载全尺寸图像
图19. 各个断裂试样底部表面的空间变化（Malama和Kulatilake，2003年）。
下载：下载高分辨率图像（685KB）
下载：下载全尺寸图像
图20. 节缝正交加载试验结果及不同本构模型的拟合。（a）Grd-1，（b）Grd-2，（c）Grd-3。
从裂缝的正交加载试验中获得了机械参数（表5），并系统地比较了不同的本构模型。如图20所示，使用BB双曲线模型（kn0, dmax）、经典指数模型（kn0, dmax）和幂函数模型（kn0, dmax）进行的数值模拟表明，在中等应力水平下，BB模型低估了正交变形，而指数模型则高估了它。在高应力下，BB模型达到最大闭合dmax的速度比实际观察到的慢，尽管幂函数模型在低应力下准确，但未能捕捉到渐近闭合行为。为了提高中等应力下的准确性，统一指数模型（Malama和Kulatilake，2003年）（kn0, dmax, n）和改进的双曲线模型（Yu等人，2008年）（kn0, dmax, m）引入了第三个参数来调整正交位移的演变，从而弥合了传统模型之间的差异。然而，这些模型仍然存在局限性：（1）初始刚度kn0在接近零应力时缺乏明确的物理意义；（2）额外的参数增加了敏感性，其物理意义尚不清楚，需要进一步研究，这限制了其更广泛的工程应用。

表5. 实验室测试中岩石样本的参数值（Malama和Kulatilake，2003年）。
岩石样本
边长（mm×mm×mm）
dmax（mm）
kn0（MPa/mm）
Grd-1
93×106×93
0.50
54.14
20.15
40.4
Grd-2
87×99×95
0.48
14.94
90.14
80.38
Grd-3
103×104×89
0.68
33.38
00.22
20.44

本研究使用真实的实验数据构建了一个数据库，其中裂缝数据直接来自裂缝正交加载试验。该岩石数据库与前面提到的模拟研究一致。本研究关注裂缝在正交加载下的变形响应特性，因此岩石数据库对裂缝正交闭合的机械响应影响可以忽略不计。使用LCDDCM方法（参数k=5）计算了在不同正交载荷下的变形响应。图20中的红色曲线展示了结果。比较分析表明，LCDDCM方法比BB模型、经典指数模型和幂函数模型更能匹配实验数据，从而显著提高了计算精度。

4.2. 节缝切向直剪试验示例
参考文献（Asadollahi等人，2010年）对红砂岩样本在各种正交应力下进行了系统的直剪试验，使用了人工制备的粗糙裂缝样本。在1 MPa的应力下对完整样本施加剪切直至破坏，然后移除正交应力并重置剪切执行器。试验使用了RDS-300伺服液压系统，该系统具有独立的正交（500 kN）和剪切（300 kN）载荷执行器，并用LVDT（±0.025 mm）测量了切向位移。图21显示了离散数据点的非线性剪切滑动行为。

下载：下载高分辨率图像（584KB）
下载：下载全尺寸图像
图21. 直剪试验结果及不同本构模型的拟合。（a）σn =0.5 MPa，（b）σn =1.0 MPa，（c）σn =1.5 MPa
切向应力-位移曲线在峰值前后都表现出显著的非线性，峰值前的增长模式各不相同。Barton-Bandis（BB）模型能够准确预测峰值剪切强度，但在估计峰值切向位移、峰值后行为、膨胀和表面退化方面存在局限性。为了解决这个问题，Asadollahi等人（2010年）提出了改进的BB模型。图21将模型拟合结果与实验数据进行了比较，表明由于裂缝表面异质性和实验复杂性，两种模型仍然难以完全捕捉峰值前的切向应力-位移行为。

本研究使用真实的实验数据构建了一个数据库，其中裂缝数据直接来自直剪试验结果，岩石数据与前面讨论的模拟研究一致。使用LCDDCM方法（k=5）迭代求解了在不同正交应力条件下断裂岩体的剪切变形响应，结果以图21中的红色曲线显示。如图所示，LCDDCM方法准确捕捉了在不同正交应力条件下断裂岩体的非线性力学行为。

5. 讨论
本研究提出的用于断裂岩体变形分析的无模型数据驱动计算力学是基于原始DDCM框架（Kirchdoerfer和Ortiz，2016年）开发和改进的。虽然传统的DDCM受到固体元素离散化策略的限制，因此仅适用于连续介质问题，但本工作将断裂岩体建模为块-缝二元系统，通过固体元素和Goodman缝元素进行离散化，从而明确捕捉了它们的不连续变形特性。该方法不需要本构方程，从而避免了由于主观选择本构模型和参数反演而产生的不确定性。相反，它直接使用实验材料数据进行计算，从而提高了结果的客观性和可靠性，同时保持了实验数据的真实性。值得注意的是，即使在存在高度复杂的材料行为（例如非线性和各向异性）的情况下，核心解决过程也不需要额外的假设。

然而，仍存在一些局限性。首先，随着数据库规模或积分点数量的增加（特别是在三维情况下），线性方程求解器（子问题PM）和全局最近邻搜索（子问题PC）的计算成本显著增加。其次，当前方法主要适用于非线性弹性阶段，要求岩石和裂缝的应力-应变关系一一对应；因此，尚未考虑岩石的塑性行为和裂缝的峰值后阶段。第三，由于数据驱动方法本质上依赖于现有的数据库进行变形预测，数据中的不确定性可能会直接影响结果的准确性和鲁棒性，超出数据库覆盖范围的推断可能会引入预测误差。

为了解决这些问题，未来的研究可以朝三个方向扩展：首先，可以采用GPU并行化（Chen和Gagan，2012年）来加速矩阵求解，并通过优化数据库构建和索引策略（如使用k-d树（Eggersmann等人，2021年）或分层算法（Hui等人，2022年）来提高最近邻搜索的效率。其次，可以引入历史状态变量或累积的塑性变量（Eggersmann等人，2019年；Ciftci和Hackl，2022年）为每个应变历史构建单独的数据库，并在求解过程中动态更新材料数据库，从而扩展该方法以捕捉岩石的塑性和裂缝的峰值后行为。第三，可以结合概率方法和不确定性量化技术（Liu和Dai，2018年；Liu和Dai，2021年）来系统评估数据不确定性对模拟结果的影响，提高方法在复杂材料行为建模中的鲁棒性和可靠性。

6. 结论
本研究提出了一种用于断裂岩体变形分析的无模型数据驱动计算力学方法。该方法绕过了构建本构模型的传统过程，直接利用实验获得的岩石应力-应变、裂缝正交应力-位移和裂缝切向应力-位移数据，实现了断裂岩体机械响应的精确解。通过几个典型示例和真实实验数据验证了该方法的有效性和计算精度。研究的主要结论如下：
（1）通过将数据驱动求解器集成到传统的断裂岩体有限元显式分析框架中，可以使用岩石和裂缝数据库准确再现断裂岩体的变形响应，从而消除了传统方法中对岩石和裂缝本构模型的依赖。Goodman缝元素的集成促进了详细的建模和断裂岩体的高效数值计算。
（2）该方法精确捕捉了裂缝网络的不连续变形特性，从而更准确地表示了岩体的变形行为。在涉及贯穿裂缝、随机裂缝和工程规模问题的数值实验中，该方法展示了优越的计算精度。
（3）采用分层和分区建模策略来分析工程规模问题：对于V类和IV类岩体，明确表示裂缝以准确捕捉不连续性的影响；而对于III类岩体，等效连续体方法可以有效表征宏观力学行为。这种策略提高了计算效率，同时保持了精度。
（4）通过使用局部凸集数据驱动计算力学，该方法已成功应用于噪声较大的非线性实验数据集，准确再现了给定数据中裂缝的非线性弹性变形行为，从而证明了其在处理非线性力学响应问题方面的显著优势。

CRediT作者贡献声明
刘瑶蕾：撰写 – 审稿与编辑。
姜庆辉：撰写 – 审稿与编辑，监督。
杨殿森：撰写 – 审稿与编辑，监督，资金获取。
冯帆定：撰写 – 原始草稿，软件，方法论，概念化。
刘东：撰写 – 审稿与编辑。

未引用的参考文献
Li等人，2022年；Li等人，2022年；Li等人，2024年。