CN109448122B - 对包含结构面的非连续体进行变形量测的方法及应用 - Google Patents

Abstract

本发明涉及对包含结构面的非连续体进行变形量测的方法，方法步骤如下：(1)对给定的含有结构面的非连续体，对其拍照并输入有限元软件中，分别对非连续体的不同结构面单独生成三角形单元，并输出对应的节点文件和单元文件；(2)生成内点；(3)求解内点位移；(4)求解三角块体顶点变形函数；(5)对顶点位移换算得到物理位移；(6)对位移进行可视化操作。本发明的方法能专门针对非连续体的变形测量，更好的反映非连续体的位移情况，从而使该方法能够更好的应用于含有节理或裂纹面的岩石试样的变形量测方面。

Description

对包含结构面的非连续体进行变形量测的方法及应用

技术领域

本发明属于非连续体变形处理技术领域，尤其涉及对包含结构面的非连续体进行变形量测的方法及其应用。

背景技术

岩石是由一种或多种矿物组成的固态集合体，作为地质作用的产物，往往具有非均质性、各向异性和非连续性的特点。其中，非连续性表现在岩石中具有不规则的结构面。此外，岩石的破坏本质上也是一种从连续到非连续的过程，岩石内部的微结构面逐渐增长、串联，逐渐变为宏观显著的裂纹面，最终使岩石发生破坏。为更深入理解岩石的变形破坏机理，需要对岩石的非连续变形及力学特性进行定量的分析。

数字图像相关法(Digital Image Correlation,DIC)最先由日本的I Yamaguchi和美国南卡罗来纳大学的Peters等科研人员分别独立提出。不同的文献给这个方法做了不同的名称，例如：数字散斑相关方法(digital speckle correlation method，DSCM)、纹理的相关性(texture correlation)、计算机辅助散斑干涉法(computer-aided speckleinterferometry，CASI)、电子散斑摄像(electronic speckle Photography，ESP)。20世纪初针对DIC，主要是基础性的研究，21世纪后，随着数码相机等数字图像采集设备的大量使用和计算机硬件性能的大幅提高，DIC作为一种基于计算机视觉技术的图像处理方法，凭借其具有全场和局部变形测量、非接触测量、对场地无特殊需求等优点，得到了迅速发展。

自DIC提出以来，国内外不少学者在此做了大量的研究。国际上，Konagi等使用激光辅助层析照相技术(LAT)分析了浸没在同样折射率液体中的碎玻璃颗粒的运动图像和全场变形；White等结合粒子图像测速(PVI)和DIC提出了一种测试土壤变形非接触测量的新方法；Gutberlet等通过PVI识别了被动土压力下的两种土壤剪切带的形成；Houda等通过DIC分析了软土地基中桩在竖向循环荷载作用下的沉降机理。如果测试对象是曲面，或者测试对象在加载后发生了三维变形，则二维的DIC方法不再适用，P.F.Luo等提出一种使用两个CCD相机获得对象表面位移的3D-DIC。在国内，芮嘉白等为提高DIC的处理速度和测量精度，提出了一种新的改进的数字散斑相关方法(十字搜索法)，由于散斑分布的随机性，该方法存在误判的可能。

对于非连续体，如图1，假设非连续体Ω存在一个非连续面τ，定义Ω₁处非连续面为τ⁺，Ω₂处非连续面为τ-，当使用数字图像相关法(DIC)进行变形分析时，对于非连续面处的像素点，将同时包含两个子块体(Ω₁和Ω₂)的变形信息，两个子块体的位移往往不同，使得DIC分析的非连续面的位移无法反映真实的情况，因此，针对这些问题，提供一种扩展DIC在求解非连续体变形方面的能力，能够对含有节理的岩体、含有裂纹面的试样等的变形进行变形量测的方法，具有重要的现实意义。

发明内容

本发明的目的在于克服现有技术的不足，提供一种扩展DIC在求解非连续体变形方面的能力，能够对含有节理的岩体、含有裂纹面的试样等的变形进行变形量测的方法。

本发明解决其技术问题是采取以下技术方案实现的：

对包含结构面的非连续体进行变形量测的方法，所述方法步骤如下：

(1)生成有限元网格

对给定的含有结构面的非连续体，对其拍照得到照片，将图片输入电脑，提取图片尺寸数据及结构面的位置信息，并输入有限元软件中，在软件中标明结构面的位置，分别对非连续体的不同结构面单独生成三角形单元，并输出对应的节点文件和单元文件；

(2)生成内点

根据节点文件在单元文件中选取三个节点组成三角块体，并计算所有三角块体内点坐标，如公式(1)所示：

其中，(x_in,y_in)为三角形内心坐标，ω为单元内点权重，(x_i,y_i)为三角形内点坐标，(X_i,Y_i)为内点对应顶点坐标；

(3)求解内点位移

在非连续体发生变形前采集的图像上取内点坐标作为变形前的像素点，且以内点坐标为中心选定标志区域；在变形后的图像上以内点坐标为中心选定搜索区域，其中，标志区域、搜索区域分别通过公式(2)、(3)确定：

t＝2×η+1 (2)

s＝2×(η+ξ)+1 (3)

其中，t为标志区域的尺寸，s为搜索区域的尺寸，η为标志区域半边大小，ξ为搜索区域半边大小，且公式(2)、(3)的单位均是像素单位；

在搜索区域中选一像素点生成标志区域等大小的矩阵，搜索区域中生成的矩阵为目标图像矩阵，计算标志区域和目标图像两个矩阵的差异度，差异度计算公式如(4)、(5)所示

S＝P-Q (4)

其中，D为差异值，P为目标图像矩阵(t×t)，Q为标志区域矩阵(t×t)，S为差异度矩阵，s_ij为P、Q中对应元素的差值(p_ij-q_ij)，其中，矩阵的元素代表数字图像每个像素的颜色值，此处采用RGB颜色空间；

计算完搜索区域中所有像素点的差异度，采用遍历搜索法求解最小差异度，差异度最小的偏移量就是内点对应的位移；

(4)求解三角块体顶点变形函数

假设三角块体位移满足线性函数，通过对内点的位移值求解系数，并将系数带入对应线性位移函数可得三角块体顶点位移矢量函数；若同一个顶点只有一个位移值，则三角块体顶点位移即为该唯一位移值；若同一个顶点有多个位移值，则三角块体顶点位移为多个位移值的算术平均值，如公式(6)；

其中，

为三角块体顶点i的位移，n为共享该顶点的三角块体数，

为对该顶点包含的n个三角块体顶点的位移矢量求和；

(5)对顶点位移换算得到物理位移

对步骤(4)中得到的像素位移

进行换算，得到物理位移，换算公式如(7)：

其中，

是物理位移，

是对应的像素位移，

是换算系数，单位是米/像素；

(6)对位移进行可视化操作

将位移值附在对应节点上生成位移矩阵，通过软件以云图的形式显示出来，最后输出位移云图和物理位移值。

所述步骤(4)中三角块体顶点位移矢量函数的求解步骤如下：

a.三角块体三个内点的位移函数如式(8)，其中，u、v为x方向位移和y方向位移，(x₁,y₁)、(x₂,y₂)、(x₃,y₃)为三角块体三个内点的坐标，a、b、c、e、f、g为系数，通过式(8)变形至式(9)、(10)，并求解系数a、b、c、e、f、g；

根据解得的系数，三角块体顶点对应的位移可表示为如式(11)所示：

式中，U、V为x方向位移和y方向位移，(X₁,Y₁)、(X₂,Y₂)、(X₃,Y₃)为三角形块体三个与内点对应顶点的坐标，a、b、c、e、f、g为通过公式(9)、(10)解得的系数。

所述步骤(3)中的RGB颜色空间可以采用灰度空间的灰度值，转换公式如(12)所示：

g(x,y)＝0.299×f_r(x,y)+0.587×f_g(x,y)+0.114×f_b(x,y) (12)

式中，g(x,y)为像素点的灰度值，f_r(x,y)、f_g(x,y)、f_b(x,y)为像素点的RGB颜色值。

所述步骤(1)中采用的有限元软件为ANSYS或ABAQUS。

所述步骤(3)中在确定标志区域和搜索区域时，为使内点可以代表所在三角块体的位移，应使区域内所有像素点位于内点所在三角块体内部。

所述步骤(2)中生成的内点坐标若不是整数，则需四舍五入取整数。

对包含结构面的非连续体进行变形量测的方法应用于含有节理或裂纹面的岩石试样的变形量测。

本发明的优点和积极效果是：

1、本发明通过对非连续体进行图像处理，通过有限元软件进行三角网格的划分，对于不同的结构面，在划分网格时做具有针对性的非连续处理，使得到的非连续变形场更加精确，所得非连续变形位移场更加利于后期分析非连续体的变形；而传统的图像处理方法假设非连续体表面为不含结构面的整体，通过计算得到的位移值与实际发生的位移有明显的误差；

2、本发明是以非连续体对象为基础，建立对应的非连续三角网格模型，每个结构面被划分成有限元网格，非连续体的复杂变形场由三角网格的分片变形来表征，相比于只能针对图像做整体变形分析的传统图像处理方法，本发明的方法能专门针对非连续体的变形测量，更好的反映非连续体的位移情况。

附图说明

以下将结合附图和实施例来对本发明的技术方案作进一步的详细描述，但是应当知道，这些附图仅是为解释目的而设计的，因此不作为本发明范围的限定。此外，除非特别指出，这些附图仅意在概念性地说明此处描述的结构构造，而不必要依比例进行绘制。

图1为本发明背景技术中提供的示意图；

图2为本实施例提供的对包含结构面的非连续体进行变形量测的流程图；

图3为基于数字图像的非连续模型生成流程图；

图4为非连续模型生成内点流程图；

图5为求解内点位移函数流程图；

图6为求解三角块体顶点变形函数流程图；

具体实施方式

首先，需要说明的是，以下将以示例方式来具体说明本发明的具体结构、特点和优点等，然而所有的描述仅是用来进行说明的，而不应将其理解为对本发明形成任何限制。此外，在本文所提及各实施例中予以描述或隐含的任意单个技术特征，仍然可在这些技术特征(或其等同物)之间继续进行任意组合或删减，从而获得可能未在本文中直接提及的本发明的更多其他实施例。需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、系统、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。

需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。

下面就结合图2至图6来具体说明本发明。

实施例1

图2为本实施例提供的对包含结构面的非连续体进行变形量测的流程图；图3为基于数字图像的非连续模型生成流程图；图4为非连续模型生成内点流程图；图5为求解内点位移函数流程图；图6为求解三角块体顶点变形函数流程图；如图2～6所示，本实施例提供的对包含结构面的非连续体进行变形量测的方法，所述方法步骤如下：

(1)生成有限元网格

对给定的含有结构面的非连续体，对其拍照得到图片，将图片输入电脑，提取图片尺寸数据及结构面的位置信息，并输入有限元软件中，在软件中标明非连续结构面的位置，分别对非连续体的不同结构面单独生成三角形单元网格，并输出对应的节点文件和单元文件，节点文件包括节点的编号和节点的坐标，单元文件包括单元的编号和组成该单元的节点编号；在此步骤中，采用的有限元软件可以为ANSYS或ABAQUS；

(2)生成内点

根据节点文件在单元文件中选取三个节点组成三角块体，计算三角块体内心坐标(x_in、y_in)，其中三角形的内心坐标公式在教科书上即可查到，属于现有技术，在此基础上，计算三角块体内点坐标，内点坐标如公式(1)所示：

其中，(x_in、y_in)为三角形内心坐标，ω为单元内点权重，(x_i,y_i)为三角形内点坐标，(X_i,Y_i)为内点对应顶点(即组成三角块体的节点)坐标；需要注意的是，在此实施例中，若生成的内点坐标若不是整数，需四舍五入取整数，这是因为后续计算内点位移时存在颜色值的像素点坐标均是整数，非整数的部分不存在颜色值，若此处不取整数，则会由于没有对应的颜色值而无法计算内点的位移；

经过计算，得到所有三角块体的内点坐标；需要注意的是，若没有得到所有三角块体的内点坐标，则选取未计算内点坐标的三角块体再次计算，直至得到所有三角块体的内点坐标；

在此需要说明的是，一个三角形块体有三个顶点，为使顶点含有的位移信息能够表征所在三角形块体的位移，取三角形内心和顶点连线上的一点，称为内点，需要注意的是：在本申请中，与三角形顶点对应的内点即为三角形内心和该顶点连线上的内点，因此，一个三角形块体的三个顶点对应着三个内点，所以，对应一个三角形块体，内点坐标(x_i,y_i)中的i分别取值1、2、3；权重为表征三角形内心、内点和对应顶点位置关系的参数，内点的坐标取决于权重，权重是用来表示选择的三角形内点与内心和对应节点间位置的一个描述量，在具体的实施例中，权重可为0-1之间的任意数字。

(3)求解内点位移

在非连续体发生变形前采集的图像上取内点坐标作为变形前的像素点，且以内点坐标为中心选定标志区域；在变形后的图像上以内点坐标为中心选定搜索区域，在确定标志区域和搜索区域时，为使内点可以代表所在三角块体的位移，应使区域内所有像素点位于内点所在三角块体内部；其中，标志区域、搜索区域分别通过公式(2)、(3)确定：

t＝2×η+1 (2)

s＝2×(η+ξ)+1 (3)

其中，t为标志区域的尺寸，s为搜索区域的尺寸，η为标志区域半边大小，ξ为搜索区域半边大小，且公式(2)、(3)的单位均是像素单位；

需要说明的是，η和ξ的值人为确定，即人为确定标志区域尺寸和搜索区域的尺寸，需要注意的是，区域尺寸过大会增加计算时间，尺寸过小会影响计算精度，因此，在人为确定时需要多方面考虑及尝试。

在搜索区域中选一像素点生成标志区域等大小的矩阵，生成的矩阵为目标图像矩阵，计算标志区域和目标图像两个矩阵的差异度，差异度计算公式如(4)、(5)所示

S＝P-Q (4)

其中，D为差异值，P为目标图像矩阵(t×t)，Q为标志区域矩阵(t×t)，S为差异度矩阵，s_ij为P、Q中对应元素的差值(p_ij-q_ij)，其中，矩阵的元素代表数字图像每个像素的颜色值，此处采用RGB颜色空间，根据不同的要求可以采用其他颜色空间，例如可以采用灰度空间的灰度值，转换公式如(12)所示：

g(x,y)＝0.299×f_r(x,y)+0.587×f_g(x,y)+0.114×f_b(x,y) (12)

式中，g(x,y)为像素点的灰度值，f_r(x,y)、f_g(x,y)、f_b(x,y)为像素点的RGB颜色值；

计算搜索区域中所有像素点的差异度，采用遍历搜索法(根据搜索区域像素点的坐标，按照先行后列的顺序)求解最小差异度，如果没有计算到最后一行，最后一列，则继续计算，所有像素点差异度计算完后，差异度最小像素点的偏移量就是内点对应的位移；

(4)求解三角块体顶点变形函数

假设三角块体位移满足线性函数，则可通过对内点的位移值求解系数，并将系数带入对应线性位移函数可得三角块体顶点位移矢量函数；若同一个顶点只有一个位移值，则三角块体顶点位移即为该唯一位移值；若同一个顶点有多个位移值，则三角块体顶点位移为多个位移值的算术平均值，如公式(6)；

其中，

为三角块体顶点i的位移，n为共享该顶点的三角块体个数，

为对该顶点包含的n个三角块体顶点的位移矢量求和；

按以上方法按照三角块体编号顺序求解三角块体顶点位移，如果没有计算完最后一个三角块体顶点位移，则继续求解。

(5)对顶点位移换算得到物理位移

对步骤(4)中得到的像素位移

进行换算，得到物理位移，换算公式如(7)：

其中，

是物理位移，

是对应的像素位移，

是换算系数，单位是米/像素，换算系数表示试件上两点间的实际长度与数字图像中相应两点之间包含像素点数量之比，计算公式如(13)

其中，

是试件上两点间的实际长度，

是数字图像中相应两点之间包含像素点数量；

(6)对位移进行可视化操作，具体的：将位移值附在对应节点上生成位移矩阵，使用相关软件，如MATLAB中片函数(patch)以云图的形式显示出来，称为位移云图，最后输出位移云图和物理位移值。

此外，需要说明的是，其中所述步骤(4)中三角块体顶点位移矢量函数的求解步骤如下：

a.三角块体三个内点的位移函数如式(8)，其中，u、v为x方向位移和y方向位移，(x₁,y₁)、(x₂,y₂)、(x₃,y₃)为三角块体三个内点的坐标，a、b、c、e、f、g为系数，通过式(8)变形至式(9)、(10)，并求解系数a、b、c、e、f、g；

根据解得的系数，三角块体顶点对应的位移可表示为如式(11)所示：

式中，U、V为x方向位移和y方向位移，(X₁,Y₁)、(X₂,Y₂)、(X₃,Y₃)为三角形三个与三个(x₁,y₁)、(x₂,y₂)、(x₃,y₃)内点对应的顶点的坐标，a、b、c、e、f、g为通过公式(9)、(10)解得的系数。

需要说明的是，在确定权重的具体值时，需要根据具体情况选择权重数值，这是因为权重值越大越靠近三角形内心越能代表该三角形块体的位移，但对三角形顶点的位移计算与真实值偏差较大；值越小，越靠近三角形顶点，越能准确反映三角形顶点的位移，后期通过数字图像相关法计算内点位移时，差异度最小点存在超出三角块体的可能，影响计算的准确性，因此，综合这些影响因素，在确定权重数值时，需要根据具体情况综合考虑确定。

本发明通过对非连续体进行图像处理，通过有限元软件进行三角网格的划分，对于不同的结构面，在划分网格时做具有针对性的非连续处理，使得到的非连续变形场更加精确，所得非连续变形位移场更加利于后期分析非连续体的变形；而传统的图像处理方法假设非连续体表面为不含结构面的整体，通过计算得到的位移值与实际发生的位移有明显的误差；相比于只能针对图像做整体变形分析的传统图像处理方法，本发明的方法能专门针对非连续体的变形测量，更好的反映非连续体的位移情况，从而使该方法能够更好的应用于含有节理或裂纹面的岩石试样的变形量测方面。

以上实施例对本发明进行了详细说明，但所述内容仅为本发明的较佳实施例，不能被认为用于限定本发明的实施范围。凡依本发明申请范围所作的均等变化与改进等，均应仍归属于本发明的专利涵盖范围之内。

Claims (7)

1.对包含结构面的非连续体进行变形量测的方法，其特征在于：所述方法步骤如下：

(1)生成有限元网格

对给定的含有结构面的非连续体，对其拍照得到照片，将图片输入电脑，提取图片尺寸数据及结构面的位置信息，并输入有限元软件中，在软件中标明结构面的位置，分别对非连续体的不同结构面单独生成三角形单元，并输出对应的节点文件和单元文件；

(2)生成内点

根据节点文件在单元文件中选取三个节点组成三角块体，并计算所有三角块体内点坐标，如公式(1)所示：

其中，(x_in,y_in)为三角形内心坐标，ω为单元内点权重，(x_i,y_i)为三角形内点坐标，(X_i,Y_i)为内点对应顶点坐标；

(3)求解内点位移

在非连续体发生变形前采集的图像上取内点坐标作为变形前的像素点，且以内点坐标为中心选定标志区域；在变形后的图像上以内点坐标为中心选定搜索区域，其中，标志区域、搜索区域分别通过公式(2)、(3)确定：

t＝2×η+1 (2)

s＝2×(η+ξ)+1 (3)

其中，t为标志区域的尺寸，s为搜索区域的尺寸，η为标志区域半边大小，ξ为搜索区域半边大小，且公式(2)、(3)的单位均是像素单位；

在搜索区域中选一像素点生成标志区域等大小的矩阵，搜索区域中生成的矩阵为目标图像矩阵，计算标志区域和目标图像两个矩阵的差异度，差异度计算公式如(4)、(5)所示

S＝P-Q (4)

其中，D为差异值，P为目标图像矩阵(t×t)，Q为标志区域矩阵(t×t)，S为差异度矩阵，s_ij为P、Q中对应元素的差值(p_ij-q_ij)，其中，矩阵的元素代表数字图像每个像素的颜色值，此处采用RGB颜色空间；

计算完搜索区域中所有像素点的差异度，采用遍历搜索法求解最小差异度，差异度最小的偏移量就是内点对应的位移；

(4)求解三角块体顶点变形函数

假设三角块体位移满足线性函数，通过对内点的位移值求解系数，并将系数带入对应线性位移函数可得三角块体顶点位移矢量函数；若同一个顶点只有一个位移值，则三角块体顶点位移即为唯一位移值；若同一个顶点有多个位移值，则三角块体顶点位移为多个位移值的算术平均值，如公式(6)；

其中，

为三角块体顶点i的位移，n为共享该顶点的三角块体数，

为对该顶点包含的n个三角块体顶点的位移矢量求和；

(5)对顶点位移换算得到物理位移

对步骤(4)中得到的像素位移

进行换算，得到物理位移，换算公式如(7)：

其中，

是物理位移，

是对应的像素位移，

是换算系数，单位是米/像素；

(6)对位移进行可视化操作

将位移值附在对应节点上生成位移矩阵，通过软件以云图的形式显示出来，最后输出位移云图和物理位移值。

2.根据权利要求1所述的对包含结构面的非连续体进行变形量测的方法，其特征在于：所述步骤(4)中三角块体顶点位移矢量函数的求解步骤如下：

a.三角块体三个内点的位移函数如式(8)，其中，u、v为x方向位移和y方向位移，(x₁,y₁)、(x₂,y₂)、(x₃,y₃)为三角块体三个内点的坐标，a、b、c、e、f、g为系数，通过式(8)变形至式(9)、(10)，并求解系数a、b、c、e、f、g；

根据解得的系数，三角块体顶点对应的位移可表示为如式(11)所示：

式中，U、V为x方向位移和y方向位移，(X₁,Y₁)、(X₂,Y₂)、(X₃,Y₃)为三角形块体三个与内点对应顶点的坐标，a、b、c、e、f、g为通过公式(9)、(10)解得的系数。

3.根据权利要求1所述的对包含结构面的非连续体进行变形量测的方法，其特征在于：所述步骤(3)中的RGB颜色空间可以采用灰度空间的灰度值，转换公式如(12)所示：

g(x,y)＝0.299×f_r(x,y)+0.587×f_g(x,y)+0.114×f_b(x,y) (12)

式中，g(x,y)为像素点的灰度值，f_r(x,y)、f_g(x,y)、f_b(x,y)为像素点的RGB颜色值。

4.根据权利要求1所述的对包含结构面的非连续体进行变形量测的方法，其特征在于：所述步骤(1)中采用的有限元软件为ANSYS或ABAQUS。

5.根据权利要求1所述的对包含结构面的非连续体进行变形量测的方法，其特征在于：所述步骤(3)中在确定标志区域和搜索区域时，为使内点可以代表所在三角块体的位移，应使区域内所有像素点位于内点所在三角块体内部。

6.根据权利要求1所述的对包含结构面的非连续体进行变形量测的方法，其特征在于：所述步骤(2)中生成的内点坐标若不是整数，则需四舍五入取整数。

7.权利要求1-6任一项所述的对包含结构面的非连续体进行变形量测的方法应用于含有节理或裂纹面的岩石试样的变形量测。

Images