使用PROS方法近似实验数据，在单轴压缩的预峰值阶段估算花岗岩的弹性模量 根纳季·科列斯尼科夫（Gennady Kolesnikov）和维塔利·舍科夫（Vitali Shekov）

《Applied Sciences》：Estimation of the Elastic Modulus of Granite Under Uniaxial Compression at the Pre-Peak Stage Using the PROS Method for Approximating Experimental Data Gennady Kolesnikov and Vitali Shekov

　　摘要
本研究探讨了2021年文献中提出的纯随机正交搜索（PROS）方法在近似岩石试样测试所得的力和位移测量数据中的应用，以花岗岩作为案例研究。主要目的是简化数据近似过程，并通过在预定义协议内减少主观因素的影响来提高实验数据分析的准确性。研究重点在于确定花岗岩试样在单轴压缩前峰变形阶段的弹性模量最大值。研究采用了岩石力学行为的数学建模和实验数据分析方法。为了近似实验数据，提出了一种改进的二维S曲线方程。通过PROS方法确定最优参数值，该方法将问题简化为求解二维目标函数最小化任务。这种优化问题的维度与实验数据点的数量无关，从而提高了计算效率。研究开发了一套系统的计算程序，用于自动化计算近似方程的参数并确定弹性模量的最大值。针对使用传统测试设备准确测量位移时遇到的挑战，提出并实施了一种数值校正程序来考虑加载系统的柔顺性。研究结果与文献报道的实验数据以及本研究中获得的数据一致。该方法论和发现可以适用于分析混凝土作为天然岩石材料的人工类似物的性质。

1. 引言
1.1. 研究问题
本文提出了一种近似实验数据的方法，该方法可以减少在确定岩石单轴压缩弹性模量时主观因素的影响。表示有无侧压力情况下岩石轴压试验结果的传统方法是处理测量结果后得到的应力-应变曲线，这些结果最初表示一组载荷-位移坐标中的点。应力-应变曲线结合其他依赖关系，对于建模和分析岩石行为所需的关键力学属性是必要的[1,2,3]。岩石的起源条件及其后续历史使其包含微裂纹、中等裂纹和宏观裂纹，同时保持一定的强度。然而，如果自然或人为的影响增强并达到某些阈值，会出现新的裂纹，从而降低岩石的强度。因此，岩石作为承重材料的自然或人工结构发生失效的风险会增加。因此，预测裂纹起始的阈值应力问题就出现了[4,5]。这个阈值对应于应力-应变曲线上的一个特定点。为了确定足够准确的裂纹起始阈值应力，有必要考虑岩石的变形和失效过程。
许多关于岩石试样单轴压缩的测试（包括考虑侧压力的测试）表明，其变形和破坏过程可以分为几个阶段[6,7,8]。这些阶段和阈值应力对应于应力-应变曲线上的以下部分和点（图1）：(1) 裂缝闭合和裂缝闭合应力；(2) 弹性区域和裂纹起始应力；(3) 稳定裂纹生长和裂纹损伤应力；(4) 不稳定裂纹生长和峰值应力；(5) 失效和后峰行为。图1. 变形阶段和应力阈值。确定峰值应力相对直接，因为相应的点（图1）通常可以从实验数据中明确识别。然而，识别其他阈值应力需要更仔细的分析。
Eberhardt等人[7]确定裂纹损伤应力对应于相对体积应变达到最大值的应力水平。使用这一标准有效减少了在确定应力阈值水平时主观因素的影响。Zhao等人[9]报告的声发射测量进一步支持了这一方法的有效性。
确定裂纹起始（CI）应力是一个更为复杂的挑战。为了提高估计的准确性并深入理解岩石的变形和失效机制，建议将声发射方法与其他技术结合使用。
为了能够不同类型岩石的裂纹起始应力进行一致比较，许多研究考察了单轴压缩下裂纹起始应力与峰值应力的比率（CI/UCS）。Nicksiar和Martin[10]应用统计方法分析了376个火成岩、沉积岩和变质岩试样的测试结果。他们发现：
在无约束单轴压缩下，CI/UCS的范围为0.42到0.47，不受材料属性影响；
在受约束压缩下，这一比率增加到0.50–0.54。
Cai等人[11,12]为地下作业附近的脆性岩体开发了裂纹起始和失效应力的广义估计。
Yan等人[5]利用群论对岩石损伤进行建模，推导出脆性岩石中CI/UCS的理论表达式，范围为0.43–0.52。然而，作者指出，对于与其研究中考虑的属性显著不同的岩石，首先需要评估材料异质性以评估这一CI/UCS范围的适用性。这一观察结果突显了现有方法的不足，并强调了确定CI/UCS值需要更通用方法的需求。
针对这一问题的实际重要性，Dong等人[13]研究了裂纹成核和传播的应用方面。他们使用离散元方法[14]对水-冰相变过程中断裂岩石中的冻融诱导裂纹传播进行了数值模拟。
在相关研究中，Li等人[15]开发了一个模型，用于确定岩石单元格在单轴压缩下的强度和失效行为。他们的方法结合了现象学理论和断裂力学原理。
Lin等人[16]使用应力-应变曲线分析了断裂岩石中的裂纹闭合应力和裂纹起始应力。同样，Qiang等人[17]探索了通过最小化主观因素的影响来提高裂纹起始应力确定准确性的替代方法。
这两项研究[16,17]强调了在数据解释中减少主观性的重要性，特别是对于对实验数据的小偏差敏感的参数，如弹性模量和阈值应力。虽然使用多种方法[16]有助于减轻主观性的影响，但并不能完全消除它。
例如，在使用声发射方法研究岩石变形阶段时，可靠的数据采集需要考虑外部噪声干扰、岩石的结构和矿物学特性以及先前应力历史的潜在表现（如Kaiser效应所描述[18,19,20]）。这些因素导致了每个试样的固有独特性[21]。因此，虽然可能无法完全消除主观影响，但可以通过开发新的方法、技术和算法逐步减少其影响。本研究为这一持续的努力做出了贡献。

确定岩石弹性模量的方法与这一挑战密切相关。Bri?evac等人[22]证明，如果已知特征点的应力和应变值（见图1），可以使用方程（1）计算弹性模量E：
$$E = \frac{\sigma_{cc}\cdot\delta_{cc}\cdot\delta_{i}\cdot\delta_{ud}}{1 - \delta_{cc}\cdot\delta_{i}}$$ （1）
对由岩石试样单轴压缩测试结果构建的应力-应变曲线进行分析，可以确定关键力学参数，包括裂纹闭合应力（$\sigma_{cc}\cdot\delta_{ud}$）、裂纹起始应力（$\sigma_{ci}$）和弹性模量（E），如方程（1）所定义。然而，从测试设备获得的实验数据通常记录为一组离散的测量点，无论是应力-应变坐标系还是应力-时间坐标系。时间分辨率（即数据采集间隔）可能精细到0.01秒。因此，为了生成一个与图1中理想化表示相当的连续应力-应变曲线，需要应用近似程序。这涉及以下方法之一：
- 用适当的数学函数拟合离散数据集，以捕捉材料的微观力学行为。
- 构建一个平滑的近似曲线，准确表示岩石在整个变形范围内的物理响应。
这种近似对于可靠地识别曲线上的特征点以及以最小的主观性计算衍生力学属性是必不可少的。

1.2. 近似
近似是指开发一个能够最小偏差地再现岩石力学行为的数学模型的过程。这样的模型通常包括一个或多个用于模拟应力-应变关系的近似方程，例如在变形的前峰和后峰阶段。为了解决模型与实验观测之间的偏差问题，采用了优化方法（即寻找函数最小值的技术）。近似在岩土数据分析中提供了几个关键优势[23]，包括平滑实验数据中的噪声和测量误差；减少分析信息的量而不显著损失岩土评估的准确性；简化测试数据处理的计算；允许将数据外推到观测范围之外；识别应力-应变曲线上的特征点以及关键趋势和模式。
目前，最小二乘法是确定近似方程参数最广泛使用的方法之一[24,25,26]。然而，这种方法并不总是数值稳定的。例如，在[27]中报告了在压缩测试中解决应用问题时的不稳定性，分析产生了不可接受的结果——可能是由于数据集集中在狭窄的应力或应变范围内。
除了最小二乘法之外，当实验数据显示异常值、异常误差分布或异方差性时，还采用了其他技术进行参数估计[28,29]。然而，这些方法并不总是在算法复杂性、多功能性和数值结果的准确性之间取得可接受的平衡。通过引入参数无关的正交随机搜索（PROS）方法[30]，这一局限性得到了显著缓解，该方法在处理噪声或稀疏数据集时保持了鲁棒性。当与Blagojevi?模型[32]中推导出的新修改版本的S曲线方程[31]结合使用时，PROS方法作为近似工具表现出很高的有效性。本研究的目的是通过在预定义协议内减少主观因素的影响来简化数据近似过程，并提高实验数据分析的准确性。研究重点在于确定花岗岩试样在单轴压缩前峰变形阶段的弹性模量最大值。这需要足够准确的实验数据。在使用传统测试设备绘制载荷-位移曲线时准确确定位移所面临的问题背景下，本研究探讨了测量位移的数值校正可能性。

2. 方法论
2.1. 实验数据的近似
2.1.1. 近似方程
在本研究中，采用了改进的二维S曲线方程[31]（方程（2））来近似实验数据。这些参数的值使用PROS方法确定：
$$\sigma = \sigma_{cc}\cdot(\epsilon_{cc}\cdot e^{xp(\tau_n))$$ （2）
参数$\sigma_{cc}$和$\epsilon_{cc}$对应于峰值点（见图1），是通过实验确定的，并包含在初始数据集中。这些值在整个近似过程中保持不变。
以下部分详细介绍了使用PROS方法计算参数n和B的过程。
2.1.2. 近似程序
如果没有先前的实验数据（例如，来自文献资料），则将试样测试到失效。在测试过程中，定期记录位移和载荷（力）——例如，每0.01秒一次——然后分别转换为应变（ε）和应力（σ）。因此，测量的应变和应力值分别表示为$\bar{\epsilon}$和$\bar{\sigma}$，适当的下标用来指代特定的数据点。
在所考虑的问题中，初始数据必须包括以下内容：
- 峰值点的测量应力值（$\bar{\sigma}_{cc},\bar{\epsilon}$）；
- 在几个中间点$\bar{\epsilon}_{i} < \bar{\epsilon}_{cc}$，$i = 1, 2,\ldots, m, m_l < m \leq m_u$测量的应变值；例如，对于花岗岩，$m_l = 4$，$m_u = 10$，这些值确定的是$ \bar{\epsilon}_{l}, \bar{\epsilon}_{u}$的下限（$m_l$）和上限（$m_u$）。
这些数据表示为一组对：
$$\left(\bar{\sigma}_{0},\bar{\epsilon}_{0}; \bar{\sigma}_{1},\bar{\epsilon}_{1}; \ldots, \bar{\sigma}_{m},\bar{\epsilon}_{m}; \bar{\sigma}_{c},\bar{\epsilon}_{c}\right)$$ （3）
这里$\bar{\sigma}_{0} = 0$，$\bar{\epsilon}_{0} = 0$。上述中间点数量m是通过实验确定的；在本研究中，m=4。近似曲线段的曲率取决于具体岩石类型的力学特性。一般来说，某一部分的曲率越大，就越需要更多的点来准确表示它。需要单独进行研究，以针对不同岩石类型制定精确的推荐方案。

PROS算法确定了参数n和B的值，这些值可以将计算得到的应力值????（来自方程（2））与测量得到的应力值?????????之间的差异降到最小，其中i=1，…，m。近似的准确性是通过目标函数???????来评估的，其定义如下：
???????=??1+??2+…+????,0. （4）

其中
????=(??????????????)2,??=1，…，??
这里，????和?????????分别表示通过方程（2）计算得到的应力值和测试过程中测量的应力值。

PROS算法解决了???????的最小化问题，其中自变量是参数n和B。能够使???????的值最小化（理想情况下为零）的n和B的最优值，就是问题的解决方案。因此，确定近似方程（方程（2）中的参数n和B，归结为最小化一个二维目标函数，而与上述实验点的数量m无关。

参数n和B可能取值的上下界是通过实验建立的。对于花岗岩，典型的界限是1 ≤?? ≤50和1 ≤?? ≤50。

使用方程（2）和选定的?????（对于i = 1，…，m），计算预测的应力值????。点（?????，σi）位于近似曲线上。

计算残差：对每个i = 1，…，m计算???? =(??????????????)2。

评估目标函数：将????的值相加，得到???????（4）。

使用PROS进行优化：使用PROS算法[30]寻找???????的最小值。

随着实验数据集中指定的中间点数量m的增加，近似的准确性会达到一定的极限，这是数值方法的典型特点。可以通过试验计算来确定中间点的合理数量。这些点的均匀分布不是强制性的，但在自动选择时是首选。未来，一个集成到测量结果处理程序中的额外算法将能够自动选择指定的点，从而减少主观因素对岩土分析结果的影响。

无论中间点的数量如何，搜索都是在参数n和B的二维范围内进行的，其上下界分别为：??????? ≤?? ≤???????和??????? ≤?? ≤???????。这些界限是基于解决类似问题的经验和研究目标而的经验性确定的。

根据PROS算法[30]，在求解过程中只有一个参数（n或B）会发生变化。如果参数n的变化使目标函数（4）的值减小，则搜索将继续在同一方向n上继续。如果参数n的变化使目标函数的值增大，则切换到正交方向B进行搜索。如果目标函数的值增大，则切换回方向n进行搜索，依此类推。

迭代基于最大尝试次数或结果的准确性终止。以下示例说明了这种近似方法。

2.1.3 示例1：确定近似方程的参数
我们考虑花岗岩变形的预峰值阶段，其实验数据是从应力-应变曲线中获得的，这些数据以及其他结果是通过Aboayanah等人[33]的研究中的数字图像相关性（DIC）技术获得的。实验数据由初始点和峰值点的应力及应变值以及在四个中间点的应力及应变值表示（表1）。在这个例子中，m=4。

表1. 实验数据*。基于解决类似问题的经验，我们为参数n和B的可能值建立了近似界限：1 ≤?? ≤50；1 ≤?? ≤50。

生成一对随机数来初始化参数的初始值，例如，?? = 2.909和?? = 20.806。这些是参数n和B的起始值。

ε的值初始化为一个随机数（???0 ≤?? ≤?????）。应力使用方程（2）计算得出。

这样，就确定了初始值n、B和ε，这些值是使用PROS算法[30]开始计算所必需的。

图2显示了目标函数（4）的图形。图3显示了使用上述PROS算法达到目标的 geral 视图和轨迹细节。

图2. 目标函数的图形。最小点位于左侧的“峡谷”中。

图3. 根据PROS算法的运动轨迹。（a）算法从坐标为B = 20.807，n = 2.909的点开始。最小点（×）位于图形的左侧“峡谷”中（B = 1.3398，n = 20.4523）。（b）详细显示了最小点附近的轨迹。结果（图3）反映了到达目标的路径和PROS算法的准确性。在1.05 ≤?? ≤1.5范围内的轨迹细节（图3b）展示了随着接近目标，自动逐步减少的步长和优雅性。在这个例子中，PROS算法得出的参数值为n = 20.4523和B = 1.3398，这些值是使用方程（2）进行计算所必需的。

尽管进行的测试数量有限，但使用PROS算法确定参数n和B的方法具有以下优势，值得进一步研究：
- 参数仅使用测量结果（表1）和方程（2）确定，不需要任何额外的参数。
- 对于近似方程（2）中参数??和??的值的搜索，简化为解决二维目标函数???????（4）的最小化问题，而与实验点的数量无关。
- 确定参数n和B的过程可以自动执行，无需用户干预，从而减少主观因素对分析结果准确性的影响。

2.1.4 示例2：近似
这个例子基于与上述示例1相同的实验数据[33]。

使用找到的参数??和??的值（图3），以及实验值?????和?????????（表1），我们指定了方程（2），在所考虑的例子中方程（2）的形式如下（5）：
??=170.8?(??0.00399?e?x?p?(120.4523?(1?(??0.00399)20.4523)))1.3398。

图4显示了表1中的实验数据（圆圈）和拟合曲线（5）。

图4. 实验数据（圆圈）和拟合曲线（红色）。如果根据方程（5）计算应力时的实验应变值??与实验值??的偏差在+5%（海军蓝）和?5%（绿色）范围内，则显示了??的界限。图4表明方程（5）在变形的预峰值阶段相当准确地近似了实验数据；小的偏差是由于数字四舍五入不可避免的影响。高精确度的近似允许使用近似方程（5）通过微分来确定弹性模量，如下一节所示。

2.2 确定最大预峰值切线模量

2.2.1 必要的说明
从物理角度来看，弹性模量是材料刚度的特征。从几何角度来看，弹性模量在数值上等于应力-应变曲线上切线倾斜角的切线值。岩石的裂隙越少，其刚性越高。

在变形的第一阶段（图1），随着负荷的增加，外部力用于闭合裂隙的工作比例逐渐减少，相应地，用于变形矿物颗粒、它们的集合体以及它们之间连接的力的工作比例增加。少数最弱的颗粒和颗粒接触点可能会被破坏，这一点通过声发射测量得到证实[18]。裂隙面收敛并闭合，岩石变得更加致密，刚性增加，并在第二阶段的某个点达到最大值（图1）。在这一点上，曲率的符号发生变化；在应力-应变曲线上，向下的凸性消失（切线斜率停止增加），出现向上的凸性（切线斜率减小）。

从物理角度来看，曲率符号的变化意味着在达到其最大点后，切线模量减小，这只有在出现新的裂隙、最弱的（或承受最大负荷的）颗粒及其集合体被破坏时才可能发生。可以合理假设，在第二阶段的这一点（图1），所有裂隙都闭合了，因此岩石密度最大，切线模量将大致对应于完整岩石的切线模量。因此，岩石在这种状态下的机械行为类似于胡克定律描述的理想弹性材料的行为。严格来说，这种状态只对应于应力-应变曲线上的一个点，如果已知足够准确的近似方程，就可以对此进行正式的搜索（在所考虑的例子中，这是方程（5））。

从数学角度来看，弹性模量的切线值在数值上等于函数（2）（或所考虑例子中的函数（5）的一阶导数：
??=??????(??)????? （6）

使用方程（2），经过变换后得到：
??=????????(???????e?x?p?(??????))????????,???=(1?(??????)??) （7）

预峰值分支是一个S形曲线，有一个拐点（图1）。在这一点上，切线模量（7）达到最大值，然后在峰值点开始减小到零。可以通过解以下方程（8）找到切线模量（7）达到最大值的?? =????：
??????(??)?????=??2????(??)?????2=0 （8）

将方程（2）代入方程（8）后，我们可以确定弹性模量达到最大值的变形（9）[31]：
????=?????? ? ? ??1+???1?√4??????+(???1)22?????? ? ? ??1?? （9）

通过代入方程（9）中的?? =????，可以确定切线模量达到最大值的应力????。通过代入值?? =????（9），利用方程（7）可以确定最大预峰值切线模量。

2.2.2 示例3：预峰值切线模量
这个例子基于与上述示例1和2相同的实验数据[33]。

通过将上述（第2.1.3节）定义的数值?????????、??和??代入方程（7），可以计算预峰值阶段（0
（a）应力???(??)（红色）和切线模量???(??)（灰色）与轴向应变的依赖关系。（b）如果根据方程（5）计算应力??时的实验应变值??的偏差在+5%（深蓝色）和?5%（绿色）范围内，那么切线模量与归一化应力（???(??/????)）的依赖关系。

表2. 初始数据和计算结果。对于所考虑的例子，切线模量的值来自文献；它们是在应力-应变曲线的0.45????和0.55????边界范围内通过两种方法实验确定的。使用第一种方法（传感器）得到52.57 GPa的值，使用第二种方法（数字图像相关性）得到51.79 GPa的值[33]。

相比之下，所提出的方法在0.55σc的应力水平下得到的弹性模量值为51.84 GPa。切线模量在坐标为 ε = 0.00325，σ = 138.6 MPa 的点达到最大值 56.2 GPa（见表 2）。如上所述，弹性模量可以使用方程（1）计算，该方程与公式 ?? =?????????????????? 类似。需要注意的是，使用此公式计算出的模量可以称为弦弹性模量。这个术语通过与“切线弹性模量”和“割线弹性模量”的类比来反映问题的几何特性。为了便于比较，图 6 以图形形式展示了弦弹性模量的计算结果，假设应力-应变曲线被分成了 s 个段。黑色点对应于在应力水平 0.55???? 时的弦弹性模量值 51.79 GPa；绿色点对应于使用所提出的方法在同一应力水平下获得的切线弹性模量（51.84 GPa）。图 6 中的红色点表示切线弹性模量的最大值（56.2 GPa）。随着段数的增加，使用所考虑方法得到的模量值趋于一致，这证实了结果的可靠性。图 6. 使用所讨论方法得到的弹性模量值的收敛性（具体解释见上文）。

2.3. 岩石变形的线性弹性阶段
第 2.2 节的分析表明，在满足条件（8）的点，切线模量最大，这间接证实了裂纹闭合。因此，该点是可以确定裂纹闭合应力（???????，见图 1）的候选点。然而，如上文（2.2.1 节）详细描述的，在应力-应变曲线上经过该点之后，切线模量在峰值应力处降至零，并在峰值后阶段变为负值。这种情况只有在出现新裂纹、部分失效以及最弱颗粒及其团聚体的失活时才会发生（这种力学行为在单轴压缩下的岩石试样中有所体现）。因此，该点也是确定新裂纹起始应力的候选点。因此，从形式上讲，满足条件（8）的点定义了现有裂纹闭合阶段与新裂纹形成阶段之间的界限。然而，应力-应变曲线的某一段在视觉上可以被视为直线段，这是线性弹性变形阶段的特征，如图 1 通常所示，并在图 7 中更详细地说明了这一点。图 7：(a) 实验点（小圆圈）、近似曲线以及通过拐点 1 的切线（在点 1 处曲率改变符号，向下凸性消失，出现向上凸性）。点 2 和 3 是使用方程（5）在切线和近似曲线的纵坐标差等于点 1 处应力的 1%（??? =0.01?????）的条件下确定的。切线与横坐标轴的交点为 ?? = 0.0007873。(b) 当 ??? =0.01????? 时的图表片段。(c) 当 ??? =0.001????? 时的图表片段。图 7 显示，在满足条件（8）的点 1 处的切线几乎与第 2 节和第 3 节中的曲线重合。将切线和近似曲线的纵坐标差从点 1 处应力的 1% 减少到 0.1% 后，几乎线性变形的长度缩短了。显然，可以选择 2–3 区域中的任何点作为标准。这一特征是主观因素存在及其对岩土分析结果产生影响的原因。因此，使用自动化处理实验数据是提高计算可靠性和简化相应方法和算法的必要条件。本文没有详细讨论确定岩石变形线性弹性阶段精确边界的方法。然而，可以假设客观存在的点 1（其坐标 ???? 和 ???? 是自动确定的），在该点切线模量最大，这在某种程度上与岩石变形线性弹性阶段的边界有关（见图 7）。

3. 结果与讨论
3.1. 考虑加载系统刚度对测量结果的校正
试验机的加载系统刚度（LSS）显著影响在压缩加载下测量试样端位移的准确性[34,35]。本节介绍了一种校正方法，可更准确地确定岩石的变形特性，特别是在使用常规（或软）试验机时。为了最小化实际位移和测量位移之间的差异，通常采用几种方法：使用更刚性的加载系统[36]；直接将物理位移传感器附着在试样上[33]；实现基于数字图像相关（DIC）技术的虚拟位移传感器[33,37]。提高位移测量精度的另一种策略是在测试后对原始测量数据进行校正。使用未经校正的位移值可能会导致确定岩石材料的物理和力学性质时出现显著误差[33,37,38]。值得注意的是，校准方法的研究仍然非常相关，原因有以下几点：不存在用于构建加载系统的无限刚性材料；高刚性试验机并不总是可用或经济可行；即使对于高刚性机器，位移校正也可能有益——特别是在测试具有不同刚性特性的试样时。让我们考虑一种在校正常规试验机在单轴压缩下对试样施加载荷 P 时测量的位移的方法。刚度（K）定义为施加力（P）与由此产生的位移（Δ）的比值，即 K = P/Δ。如果已知 K 和 P 的值，位移可以按照以下公式计算：??=????（10）。顺应性（C）定义为刚度的倒数，即 ?? =???1 或等效地，?? =????。如果已知 C 和 P 的值，位移由以下公式给出：??=??·??（11）。设 ??1 和 ??2 分别为试验机和试样的刚度。在常规试验机上测量的位移记为 ?。测量的位移 ? 来自于加载系统和试样的变形。因此，测量位移等于 ? =?1+?2，其中位移 ?1 仅由加载系统的弹性变形引起，而 ?2 仅由试样的变形引起。重新排列方程后可得到试样位移的表达式：?2=???1（12）。将关系式（10）纳入后，方程（12）可以改写为方程（13）的形式：?2=????1??1（13）。如果忽略加载系统内部的摩擦，则作用在各个组件中的力——加载系统（??1）、试样（??2）和组合系统（P）——是相等的：??=??1 =??2。这个等式意味着，在方程（13）中，可以用试验机测量的力 P 替换力 ??1，因为在无摩擦系统的假设下所有组件都承受相同的载荷：?2=??????1（14）。从计量学的角度来看，方程（14）包括以下量：Δ 和 P——串联连接系统（加载系统 + 试样）中的测量位移和力；??1——试验机加载系统的刚度；Δ1 =????——反映加载系统变形的校正项；?2——仅由试样变形引起的校正位移。考虑到方程（11），方程（14）可以写成（15）：?2=????1???（15）。因此，方程（10）–（15）从理论上模拟了安装在试样上的虚拟伸长计的操作。这种校正程序有效地起到了虚拟伸长计的作用，提供了与直接在试样表面进行物理测量获得的位移数据相当的数据——而无需额外的硬件。为了进行计算（14），需要确定试验机加载系统的刚度（??1）。为此，考虑到方程（11）和上述定义的等式 ??=??1 =??2，我们将方程（14）重写为（16）：??2???=??????1（16）。如果试样的顺应性（??2）为零，则从方程（16）可以确定方程（12）中指定的校正值：?1=????1（17）。在实验室条件下，校正值（17）是根据顺应性为零的虚拟试样的测试结果确定的（??2 =0；在这种情况下 ??2 →∞）。这意味着力-位移曲线是在“空”的机器上绘制的（没有试样）。在这种情况下，上板和下板直接通过接触面以力 P 相互作用。如果位移 ?1 与 P 成线性关系，则加载系统的刚度是恒定的。然而，在实际条件下，力-位移曲线的初始部分是向下凸的，因为实际加载系统的刚度不是瞬间增加的，而是在短时间内逐渐增加的。加载系统中的每个额外环节都会增加其顺应性，从而降低其刚度。这意味着在测试之前必须确定加载系统的实际顺应性，以确定校正因子（17）的值。图 8 显示了本研究中使用的试验机的校正（17）与压缩力的关系图（见第 3.2 节）。曲线在初始区域向上凸，因为它是一个位移-载荷曲线，即一个反向的载荷-位移曲线。图 8. 实验曲线显示了作为测量力 P 函数的校正值（17）对测量位移的校正。这条曲线仅适用于生成该曲线的试验机及其附加设备（如果有的话）。对于由试验机及其附加设备组成的每个系统，都必须获得类似的曲线。通过近似上述顺应性为零的虚拟试样的实验数据来实际确定校正 ?1 的方程具有以下形式：如果 ??≤39?kN，则 ?1=?0.1755972×10?6???4+0.1937955×10?4???3?0.8121482×10?3???2+0.2274613×10?1???，??2=0.9998；如果 ??>39 kN，则 ?1=0.5959798×10?2???+0.1649199，??2=0.9999（18）。在方程（10）–（15）中，用 Δ 表示试验机测量的真实试样的位移。从理论到实践，我们将测量的位移重新定义为 ?m?e?a?s?u?r?e?d。在方程（12）–（15）中，校正后的位移表示为 ?2；为了方便使用，我们将校正后的位移重新定义为 ?c?o?r?r?e?c?t?e?d。然后，方程（12）变为以下形式（19）：?c?o?r?r?e?c?t?e?d =?m?e?a?s?u?r?e?d??1（19）。作用在试样上的力 ?? =??m?e?a?s?u?r?e?d 总是同时由机器测量与位移 ?m?e?a?s?u?r?e?d 测量的。这种方法用于处理下一节中介绍的 rapakivi 花岗岩测试的结果。

3.2. Rapakivi 花岗岩测试结果的近似
基于第 2 节介绍的方法，我们考虑了来自 Mustavara 矿床的五块 rapakivi 花岗岩圆柱试样的峰值前变形阶段（见表 3 和图 9）。表 3. Rapakivi 花岗岩（Mustavara 矿床，卡累利阿）的初始数据和计算结果。图 9. AGS-X 300 试验机以及测试前后的试样。试样的矿物组成如下：钾长石—58%；斜长石—7%；石英—32%；黑云母、角闪石和矿石矿物—剩余 3%。与 rapakivi 花岗岩相关的详细矿物学、地球化学和地质年代学方面在 [39,40] 中有讨论。所有试样的直径均为 42.6 毫米。按试样编号顺序列出的高度分别为 38.1、39.0、38.8、38.5 和 39.8 毫米。在 AGS-X 300 试验机上以位移控制模式进行单轴压缩测试，横梁移动速率为 0.2 毫米/分钟。施加载荷使用内部载荷传感器（容量 300 kN；精度 ±1% 全量程）测量，并以 100 Hz 的频率记录。每个试样进行了 50,000 到 60,000 次测量。位移使用前一节描述的方法和公式（19）确定。图 10 显示了 PROS 算法 [30] 为第一个试样生成的搜索轨迹。起始点用方块标记，无论初始条件如何，算法始终收敛到一个最终点。图 10. PROS 算法从不同起始点出发的搜索轨迹，用小方块表示。无论初始坐标如何，算法都收敛到相同的最终点，坐标为 n = 9.676 和 B = 2.418。对角线标记了最终点。利用已确定的参数n和B的值，根据第2节描述的方法论，并基于方程（2）进行了进一步计算。这一过程得到了五个试样的应力-应变曲线（图11）。图11显示了根据第2节描述的方法论构建的应力-应变曲线。实验数据点用黑色圆圈表示，这些数据点是通过校正后的测量结果获得的（19）。大的红色圆圈表示切向弹性模量?? =??????????达到最大值的点；在该点上，应力关于应变的二阶导数为零：??2????????2 =0。小的红色圆圈表示在视觉上被视为应力-应变曲线线性段的边界。灰色的斜线代表在其二阶导数为零处绘制的切线。该切线的斜率在数值上等于切向弹性模量的最大值。需要指出的是，方程（2）中的参数n和B是为每个试样分别确定的。这使得不仅可以通过对方程进行微分来分析应力-应变曲线上的任意点的切向弹性模量，还可以分析变形过程中的能量相关特性。通过后续研究中对提出的方程进行积分，可以进一步获得这些特性，从而实现更全面的机械特性描述。具体来说，给定应变水平ε下的切向弹性模量E可以通过对应变函数???(??)关于应变的一阶导数来获得。同样，通过积分（?? =∫??*0???(??)????）可以计算出在达到某个应变??*之前变形过程中积累的比应变能（单位体积的能量），这为了解材料的能量吸收能力和失效机制提供了有价值的见解。

从拉帕基维花岗岩试样的单轴压缩试验中获得的额外实验数据在表3中呈现。这些结果在图9中进行了说明，包括应力-应变响应的关键机械参数和特征点，进一步支持了所提出的分析方法的验证和应用。除了固有的地质特性外，外在的实验因素也可能导致观察到的差异。特别是，在机械测试中使用的加载系统的刚性（或硬度）会显著影响测量的弹性模量。刚性较低的测试系统可能会引入顺从效应，从而人为降低试样的表观硬度，导致较低的切向模量值。因此，为了确保对测试结果的可靠解释，建议在表征岩石力学性质时使用多种测试方法，比较不同实验设置获得的结果，并继续系统研究以进一步阐明地质历史、微观结构和宏观力学响应之间的关系[41,42,43]。

此外，本研究中讨论的方法在近似静态和动态加载条件下的实验数据方面也显示出应用潜力。这种多功能性增强了其在岩土工程中的实用性，因为在广泛的应变率和加载情景下理解材料行为对于准确建模和设计至关重要[44,45,46]，包括将这种方法适应于混凝土这种人造岩石的力学行为分析，并考虑新的方法论[47,48]。

3.3. 所提出的数据近似方法的范围和局限性
目前，所提出的方法适用于近似实验数据的范围仅限于测试花岗岩及其类似岩性。相关的应用领域还包括测试混凝土样品，因为混凝土可以作为天然岩石的人造类似物。本文中提供的实施示例和支持性论据展示了减少主观因素对单轴压缩测试中位移测量及相关应变计算准确性的影响的潜力。然而，应该强调的是，使用应力-应变曲线估计最大切向弹性模量时的主观影响并未被完全消除，而是得到了缓解。在以下方面仍存在一定的主观性：
- 应力-应变曲线函数的选择：在本研究中采用了具有两个无量纲参数（n和B）的S曲线；
- 识别应力-应变曲线上线性区域的标准，这对于准确计算弹性模量至关重要；
- 选择中间数据点（m ≥ 4）的程序，以及这些数据点在应变范围内的分布。

需要进一步的研究来开发能够进一步缓解这些限制的方法。这类研究可以集中在以下方面：
- 功能形式选择的具体标准；
- 用于识别线性区域的自动化算法；
- 选择中间点的优化策略，可能利用机器学习或高级统计技术，包括对参数值变化的敏感性分析。

4. 结论
为了提高对花岗岩及其类似地质材料的压缩测试结果近似的准确性，本研究提出了一种在峰值前阶段拟合实验数据的替代方法。所提出的方法结合了修改后的S曲线方程和纯随机正交搜索（PROS）优化方法[30]。本工作的主要发现和贡献总结如下：
- 应用了PROS进行参数识别。该研究证实了2021年文献中介绍的纯随机正交搜索（PROS）方法用于确定近似方程参数的可行性。这种方法通过减少定义实验协议内的主观因素影响，提高了实验数据分析和的客观性和准确性。
- 使用基于PROS的参数估计修改后的S曲线模型。通过使用PROS方法解决二维目标函数最小化问题来确定模型参数的最优值。值得注意的是，优化问题的维度与实验数据点的数量无关，从而提高了计算效率。
- 估计最大切向弹性模量的程序。开发并展示了一种系统程序，用于确定单轴压缩下试样变形峰值前阶段的最大切向弹性模量。获得的结果与现有研究中的发现一致，从而验证了所提出的方法。
- 根据虚拟传感器概念对测量位移进行数值校正。在利用软测试机进行单轴压缩试验时准确测量位移所面临的挑战背景下，引入了位于试样上的虚拟伸长计的概念。提出并实现了一种数值校正程序，以考虑机器的顺从性并提高测量位移的可靠性。

工程意义。在未来经过进一步研究后，可以使用PROS算法近似实验数据的方法可以在岩土测试机的软件中实现。此外，为了提高岩石测试中位移测量的准确性，本文中阐述的虚拟伸长计概念内的位移校正程序也可以在同一岩土测试机的软件中实现。这种校正程序有效地起到了虚拟伸长计的作用，提供了与直接在试样表面进行物理测量获得的位移数据相当的结果——而无需额外的硬件。

未来的研究方向。所提出的计算框架的潜在扩展包括将其适用于不同地质沉积物中岩石的比较岩土分析。进一步的发展还可能涉及将峰值后的变形行为纳入分析中，从而实现对整个应力-应变曲线范围内岩石力学响应的更全面表征。