基于DSCM-FEM的花岗岩力学参数反演方法研究

凌小康， 张敬宗， 钱 鼎， 宋义敏

(北方工业大学 土木工程学院，北京 100144)

1 引 言

岩石是矿山工程和交通工程等广泛存在的一种地质材料，岩石材料的力学参数作为岩土工程理论分析和数值计算的数据基础，其准确性对岩土工程设计、施工和维护等具有重要的影响[1,2]。因此，研究岩石材料力学参数反演具有理论和现实意义。

目前材料力学参数反演的方法主要分为两类[3]。第一类是通过加载条件和荷载信息来对材料的力学参数进行反演分析。通过测量标准试件上某个区域的应变值，结合荷载信息来计算其相应的力学参数。吉建民等[4]应用数字散斑相关方法和电测法对比测试了三种不同航空复合材料的弹性常数，证明了数字散斑相关方法能够很好地反映航空复合材料的弹性常数；耿红霞等[5]利用三维光学应变测量系统采集不同荷载下的氧化锆陶瓷图像，通过应力-应变曲线拟合得到的斜率即为弹性模量，为陶瓷材料的弹性变形行为研究提供了一种可靠的方法。易亚楠等[6]通过数字散斑相关方法测得核石墨的杨氏模量和泊松比与电阻应变片法的测量结果具有很好的一致性，验证了数字散斑相关方法的可行性与准确性。第二类是基于数字图像相关技术和数值计算相结合进行反演分析。Cooreman等[7]使用数字散斑相关方法测量带孔金属的应变场，用有限元软件模拟相同条件，得到该金属的力学参数。贺体人等[8]将数字图像相关方法和有限元模型修正技术相结合，通过短梁剪切试验获得碳纤维增强环氧树脂沿厚度方向压缩的本构关系参数。李璐璐等[9]采用修正的Drucker-Prager Cap模型对金属粉末压制成形过程进行模拟，基于 ABAQUS -MATLAB 联合仿真平台，获取混合金属粉末本构模型参数。但是针对反演目标函数选择的研究较少，因此，本文基于数字散斑相关方法和有限单元法对岩石力学参数反演进行研究。

本文以花岗岩三点弯曲实验为基础，假设材料初始本构，采用有限元模型计算获得位移场仿真值，以数字散斑相关方法DSCM(Digital Speckle Correlation Method)计算的位移测量值和有限元模型的位移仿真值之差，建立不同的目标函数，采用鱼群优化算法实现力学参数的更新，最终获得花岗岩的力学参数。

2 基于DSCM-FEM的力学参数反演方法

2.1 力学参数反演方法

根据本构关系可知，岩石在受载时产生的位移场是多个力学参数作用的结果，基于反分析的思想，采用数值计算结合试验的方法，设计了三点弯曲试验反演花岗岩多个力学参数的反演方法。三点弯曲是一种典型的非均匀变形形式，该变形能够为力学参数测量提供充分的力学响应信息[10]。本文的反演方法通过鱼群优化算法不断地调整有限元模型中的力学参数，使反演区域内有限元模型的位移仿真值与实验的位移测量值之差最小化，以此时有限元模型中的力学参数表征花岗岩的力学参数。图1给出了力学参数反演的流程。

图1 力学参数反演流程

整个力学参数反演流程通过MATLAB-PYTHON-ABAQUS联合完成，通过Matlab协同驱动各个软件，设置合理的参数范围和优化条件就可以完成力学参数的反演。

2.2 数字散斑相关方法

数字散斑相关方法通过跟踪变形前后图像中同一点的移动来获得变形信息，如图2所示，其关键在于匹配参考图像与变形图像中对应点的位置。假设在变形过程中图像的灰度具有不变性[11,12]，就能以相关函数为度量，搜索变形前后图像中最相关的点。当两个子区完全一致时，相关系数为1；完全不一致时，相关系数为-1[13]。

图2 图像变形分析原理

根据位移场的连续性假设，图像子区内任意一点的水平位移和竖直位移[14]可表示为

(1)

式中Δx和Δy为点(x,y)到子区中心点的距离，∂u/∂x,∂u/∂y，∂v/∂x和∂v/∂y为位移的一阶导数。

2.3 不同目标函数反演程序验证

为了确定不同目标函数对反演结果的影响，本文通过三个目标函数进行目标函数反演精度验证，目标函数1～3分别如式(2～4)所示，目标函数1考虑水平位移和竖直位移；目标函数2仅考虑水平位移；目标函数3仅考虑竖直位移。虚拟数据是检验反演方法可行性的有效方法[15]，用已知材料力学参数数据，输入有限元模型，将有限元模拟位移场作为虚拟的实验数据，进行力学参数反演。

(2)

(3)

(4)

针对弹性模量和泊松比进行参数反演辨识，设置弹性模量为30000 MPa，泊松比为0.3作为力学参数的输入，建立有限元计算模型，按位移荷载 0.3958 mm 加载。优化算法选择鱼群算法，因为参数反演中涉及的反问题模型是高度非线性的，鱼群算法拥有较强的鲁棒性，且有较好的全局寻优能力[16]。设置鱼群优化算法的鱼群数为4，迭代次数为25，最多试探次数为2，感知距离为 0.2，拥挤度因子为0.618，步长为0.05，不同目标函数的鱼群初始值相同[17]。考虑到弹性模量和泊松比数据的离散性问题，对输入数据进行归一化处理[18]。

不同目标函数优化过程中最优坐标变化如 图3 所示，三个目标函数的最优坐标经过前期的波动之后，逐渐趋近于参考值；但是目标函数1在优化过程中最优坐标的变化波动最大，而且目标函数1的最优坐标在优化结束后离参考值差距还很大。这是因为目标函数1包含了更多的信息，所以目标函数1对弹性模量和泊松比的变化更加敏感，即弹性模量和泊松比发生相同变化，对目标函数1的值影响更大，最终影响反演精度。目标函数值在优化过程中的变化历程如图4所示，三个目标函数值呈现一直减小的趋势，目标函数1在经过10次迭代后基本趋于收敛，而目标函数2和目标函数3仅在经过5次迭代后就已经基本趋于收敛。这是由于目标函数1的初始目标函数值最大，所以影响了目标函数收敛的速度。

随着科学技术的不断进步，传统纳税申报方式逐渐会被纳税人网络自主申报取代，纳税申报的真实性、准确性责任由纳税人承担，税务机关侧重对纳税申报的后续管理是未来税收征管的趋势。申报纳税成为纳税人的自主行为，在主观因素和客观因素共同作用下，税收申报变得具有不确定性，存在税款流失的风险，无法保证税款及时足额入库。税务机关通过纳税评估，审查纳税人申报行为的准确性和全面性，涉及税款应缴未缴或存在税收违法违纪行为的，及时让纳税人进行补缴税款、采取违法处罚或移交稽查等方式，目的是促进纳税申报质量的提升，提高纳税人的税法遵从度。

图3 目标函数在优化过程中最优坐标的变化

coordinates during the optimization process

图4 目标函数值在优化过程中的变化历程

表1给出了不同目标函数的反演结果以及反演结果的误差。综合考虑反演的收敛速度和反演结果的准确性，本文采用目标函数2，即式(3)作为花岗岩三点弯曲实验反演中的目标函数。

表1 反归一化处理后E和μ的反演

3 测量实验

3.1 试件制备

三点弯曲试件的尺寸为400 mm×100 mm×50 mm，试件的有效跨度为360 mm，如图5所示。在试件表面进行人工制斑，首先在试件表面喷涂黑漆作为底色，再向表面喷涂白漆，使其呈随机分布的白色斑点，待喷漆完全晾干以后，试件制备完成。

图5 三点弯曲试件

3.2 实验测量

本次实验系统包括加载系统和图像采集系统两部分。加载系统为电子万能试验机加载，采用位移控制方式，试件上方压头以0.05 mm/min向下运动，试件下方保持位移为0。图像采集系统使用CCD相机采集实验全过程的试件表面图像，采集速率为2 帧/s，图像分辨率为1600×1200像素，物面分辨率为0.089 mm/pixel，测量精度为0.01 pixel，实际位移精度为0.00089 mm，实验系统布置如图6所示。

图6 实验系统布置

3.3 正问题模型分析

根据实验试样尺寸和边界条件，建立有限元模型。有限元网格的单元类型为C3D8，在试件表面中间区域(100 mm×100 mm)划分更密集的单元，作为感兴趣区域，感兴趣点是感兴趣区上所有的节点，其中感兴趣区域的单元尺寸为2.5 mm×4 mm，因此共有感兴趣点41×26个，如图7所示。模型边界条件为下方两个支座固支，本文选取t=475 s时刻的试件位移场来进行参数反演，通过试件上端像素点坐标的变化得到t=475 s时刻试件发生的位移为0.1460 mm，在有限元中施加 0.1460 mm 的位移控制。目标函数中仿真值为感兴趣点集在有限元模型中的计算结果；测量值的获得首先需要通过图像分辨率和物面分辨率找到感兴趣点集在散斑图像上的对应位置，然后在数字散斑位移场中通过插值得到感兴趣点集的测量值。

图7 有限元模型

4 反演结果及分析

4.1 花岗岩力学参数反演结果

以t=475 s时刻试件位移场作为待反演时刻，设置鱼群优化算法的鱼群数为8，迭代次数为20，最多试探次数为4，感知距离为0.2，拥挤度因子为0.618，步长为0.05。同样在反演前对力学参数E和μ的值进行归一化处理，在反演结束后再对数据进行反归一化处理。表2给出了各个参数的反演结果以及反演误差。

表2 反归一化后E和μ的反演

反演结果表明，基于DSCM-FEM的参数反演方法获得的各个力学参数具有较高的精度和较好的稳定性。对在较大搜索范围内的待反演力学参数仍具有较强的识别能力，反演的精度和速度皆比较理想。其中泊松比的反演误差为8.00%，弹性模量的反演误差仅为1.98%。但是，应当认识到弹性模量和泊松比的反演结果仅是在现有条件下最优的近似结果，这是因为不同的力学参数经过不同的变化可能产生相同的位移结果。

4.2 反演分析

选择t=0 s时刻散斑图像作为参考图像，选择t=475 s时刻散斑图像作为变形后的图像，设置子区为41 pixel×41 pixel，步长为5 pixel，采用数字散斑相关方法计算得到相应位移场。图8给出了试件表面水平方向位移场云图。将反演得到的力学参数作为数值模型相应参数输入，通过正演分析，得到模型表面的位移场云图。图9给出了试件表面水平方向的位移场云图。

图8 散斑水平位移场云图

图9 有限元模型水平位移场云图

可以看出，散斑位移结果与有限元位移结果位移云图分布特征基本相同，都是试件上端压缩，试件下端拉伸，水平位移呈中心对称分布，具有良好的复杂性。局部位移分布仍然存在缺陷，主要表现在散斑水平位移场在感兴趣区域中间部分过渡不光滑，这可能是有限元模型中位移结果被平滑了，导致部分细节信息隐藏，但是因为整个感兴趣区域都参与计算，可以降低局部误差对反演结果的影响。综合来看，反演得到的力学参数较为理想。

5 结 论

岩石力学参数反演是一个多参数组合的大空间搜索问题。本文提出一种基于数字散斑相关方法和有限单元法的力学参数反演方法，以数字散斑相关方法获得试件的全场位移作为测量值，以有限元仿真的正问题分析结果作为仿真值，考察了不同的目标函数组合，采用鱼群优化算法，求解岩石试件的弹性模量和泊松比。主要结论如下。

(1) 针对力学参数反演中目标函数对反演精度影响，本文考察了不同的目标函数，结果表明基于水平位移的目标函数具有更高的反演精度和更快的反演速度。

(2) 花岗岩力学参数反演结果表明，本文提出的基于数字散斑相关方法和有限单元法的力学参数反演方法，对在较大区间下的力学参数具有较强的识别能力，且反演精度可靠，可为工程设计与分析提供支持。