CN105737769A - 一种数字体积相关算法在边界问题上的处理方法 - Google Patents

Abstract

本发明公开了一种数字体积相关算法在边界问题上的新型处理方法，该方法包括如下步骤：根据三维图像数据的灰度和边界形状定义计算区域(VOI)，利用VOI区分有效体素点和无效体素点，VOI内的体素点为有效体素点，VOI外的体素点为无效体素点；将传统子区内所有有效的体素点组成改进的子区，利用改进的子区确定计算点的整体素位移；利用改进的子区计算亚体素位移；利用改进方法得到的位移场，且只选取VOI内有效计算点的非规则窗口计算应变场。本发明方法克服了传统数字体积相关方法在边界附近的计算错误问题，提高了在边界及边界附近计算点的位移和应变计算结果的准确性。

Description

一种数字体积相关算法在边界问题上的处理方法

技术领域

本发明涉及一种数字体积相关算法在边界问题上的处理方法，属于材料与结构的变形测量领域。

背景技术

数字体积相关方法(DVC)是目前唯一能直接测量材料内部三维变形场的实验力学方法。DVC在1999年由B.K.Bay等人(BayBK，SmithTS，FyherieDP，etal.Digitalvolumecorrelation：Three-dimensionalstrainmappingusingX-raytomography[J].ExperimentalMechanics，1999，39(3)：217～226)提出后，在生物材料、泡沫材料、复合材料等领域取得了广泛的应用。

DVC是依据变形前后三维数字体图像灰度的变化确定三维位移场，进而计算应变场的一种实验方法。DVC原理上基于连续介质力学，不但要求材料相对连续，也要求材料的变形场连续。而实际的材料和结构不可避免地存在不同的内外边界，比如复合材料中的增强相，粗细骨料、颗粒、纤维等，再比如很多材料(比如混凝土、陶瓷等)中均存在气孔、裂纹等典型的缺陷，还有规则或不规则材料形状也构成外边界。传统DVC方法在跨越这些边界时，所依赖的理论基础不再成立，会出现明显的计算错误。基于体素子区的DVC在计算变形时以围绕计算点的规则子区为基本单位。计算点在边界或者边界附近时，子区会跨过边界，由于边界内外图像灰度差别和材料变形的不连续性，导致边界上及边界附近的计算点的位移和应变结果均存在明显的错误。

除了基于体素子区的DVC，还有一种基于有限元的DVC，也叫全局DVC(RouxS，HildF，ViotP，etal.ThreedimensionalimagecorrelationfromX-Raycomputedtomographyofsolidfoam[J].CompositesPartAAppliedScience&Manufacturing，2008，39(8)：1253-1265.)。基于有限元的DVC同样存在边界和非均匀变形场测量问题。为了解决该问题，JRéthoré等人在扩展有限元的基础上提出了解决三维问题的方法(RéthoréJ，LimodinN，BuffièreJY，etal.Three-dimensionalAnalysisofFatigueCrackPropagationusingX-RayTomography，DigitalVolumeCorrelationandExtendedFiniteElementSimulations[J].ProcediaIutam，2012，4(9)：151-158.)，该方法的核心思想是通过改变有限元中的形状函数来实现对不均匀变形场的描述。基于体素子区的DVC和基于有限元的DVC本身具有巨大的差别，可视为两种不同的独立方法，前者是对数据局部利用进行相关性计算，后者是对数据全局进行有限元计算。JRéthoré等人的核心思路是改变扩展有限元中的形函数来准确描述边界附近的变形场，本发明的核心思路根本性地区别于基于有限元的DVC。

发明内容

发明目的：本发明所要解决的技术问题是提供一种数字体积相关算法在边界问题上的处理方法，本发明方法是通过在边界附近区分有效体素点和无效体素点，构建不规则形状的改进子区，从而准确计算边界附近的变形场，解决了传统方法存在边界计算错误的问题。

发明内容：为解决上述技术问题，本发明所采用的技术方案为：

一种数字体积相关算法在边界问题上的处理方法，包括如下步骤：

步骤1，根据三维图像数据的灰度和边界形状确定计算区域，根据计算区域区分三维图像数据的有效体素点和无效体素点，位于计算区域内的点为有效体素点，位于计算区域外的点为无效体素点；

步骤2，整体素位移场的计算：当计算点位于边界上或者位于边界附近时，将围绕计算点传统子区内所有有效的体素点组成改进的子区，保持计算点在传统子区的位置不变，利用改进的子区确定计算点整体素的位移场；

步骤3，亚体素位移场的计算：将围绕计算点传统子区内所有有效的体素点组成改进的子区，保持计算点在传统子区的位置不变，利用改进的子区确定计算点亚体素的位移场；

步骤4，应变场计算：利用步骤2和步骤3获得的位移场，对于位于边界及边界附近的点，只选取计算区域内有效体素点的非规则窗口计算应变。

本发明方法中所指的数字体积相关是指基于体素子区计算的数字体积相关。

其中，所述边界包括内边界和外边界，所述内边界包括材料中的骨料、颗粒、纤维、气孔或裂纹或者为材料或结构的预制内边界，所述外边界包含规则形状样品的外边界和不规则形状样品的外边界。

其中，步骤1中，所述三维图像数据是指采用X射线断层照相、伽马射线断层照相、核磁共振断层照相、中子断层照相、聚焦离子束-电子束-双束显微镜或共聚焦显微镜进行三维成像得到的三维图像数据。

其中，所述X射线断层照相是指医用X射线断层照相、工业X射线断层照相、显微X射线断层照相、纳米X射线断层照相或同步辐射X射线断层照相中的任意一种。

有益效果：本发明方法不仅使测量复合材料内部夹杂和缺陷附近的准确变形场成为可能，也使测量具有复杂形状样品的位移场和应变场成为了可能；同时，本发明方法还能确定材料内部由于变形而产生的裂纹；本发明方法有效解决了传统数字体积相关方法关于边界问题存在的计算错误问题，提高了边界附近计算点在位移和应变上计算结果的准确性；最后本发明方法为数字体积相关方法在实验力学和材料领域的广泛应用起到了推动作用。

附图说明

图1为含内边界的子区示意图；

图2为将传统子区变为改进子区的示意图；

图3为用改进子区跟踪参考图像和目标图像中相同计算点(B点)的位移量；

图4为实施例中I型裂纹的示意图；

图5为实施例中样品通过X射线断层照相法获得的第一次CT扫描二维数据切片；

图6为实施例中样品通过X射线断层照相法获得的第二次CT扫描二维数据切片；

图7为对第二次CT扫描的数据施加I型裂纹位移场得到的二维数据切片；

图8为在裂纹左右两侧分别定义两个计算区域的示意图；

图9为实施例采用不同方法计算得到的样品变形后沿X方向和Y方向的位移场；其中，(a)和(d)为理论的位移场，(b)和(e)为本发明方法得到的位移场，(c)和(f)为现有DVC方法得到的位移场；

图10为实施例中采用不同方法计算得到的样品变形后Y方向的应变场；其中，(a)为理论的应变场，(b)为本发明方法得到的应变场，(c)为现有DVC方法得到的应变场。

具体实施方式

下面结合附图和实施例，对本发明的技术方案进行详细说明，本发明实施例采用写入的I型裂纹位移场来进行说明。但本发明要求保护的范围不局限于实施例表示的范围，故凡依本发明专利申请范围所述的构造、特征及原理所做的等效变化或修饰，均包括于本发明专利申请范围内。

一种数字体积相关算法在边界问题上的处理方法，包括如下步骤：

步骤1，两次原位CT扫描的三维图像数据分别为G₁(x，y，z)和G₂(x，y，z)；对第二次扫描得到的三维数据G₂(x，y，z)施加I型裂纹的变形场，裂纹示意图如图4所示，裂纹的变形场为：

$u(r,\mathrm{\θ})=\frac{K_I}{2E}\sqrt{\frac{r}{2\mathrm{\π}}}(1+v)[(2\mathrm{\κ}-1)\cos \left(\frac{\mathrm{\θ}}{2}\right)-\cos \left(\frac{3\mathrm{\θ}}{2}\right)]$

$v(r,\mathrm{\θ})=\frac{K_I}{2E}\sqrt{\frac{r}{2\mathrm{\π}}}(1+v)[(2\mathrm{\κ}-1)\sin \left(\frac{\mathrm{\θ}}{2}\right)-\sin \left(\frac{3\mathrm{\θ}}{2}\right)];$

w(r，θ)＝0

其中，(r，θ)为极坐标，将u(r，θ)、v(r，θ)和w(r，θ)转换成直角坐标下的u(x，y，z)、v(x，y，z)和w(x，y，z)，将三维数据G₂(x，y，z)的每个体素点加上裂纹场相应位置的位移，通过插值，得到变形后三维数据G₃(x，y，z)，其典型二维切片如图7所示；G₁(x，y，z)和G₃(x，y，z)将被用于数字体积相关计算；

步骤2，根据三维图像数据的灰度和边界形状定义计算区域(VOI)，根据计算区域区分三维图像数据的有效体素点和无效体素点，本例根据裂纹的位置在裂纹左右两侧分别定义两个VOI，如图8所示为具有两个VOI的二维示意图，当计算点在第一个VOI内时，第一个VOI内的体素点均为有效体素点，而第二个VOI内的体素点均为无效体素点；相反，若计算点在第二个VOI内，第二个VOI内的体素点均为有效体素点，而第一个VOI内的体素点均为无效体素点；

步骤3，整体素位移计算：将传统子区内所有有效的体素点组成改进的子区，如图2所示，并且计算点在传统子区的位置不变，利用改进的子区确定计算点的整体素位移u₀，v₀，w₀，即得到改进子区的零均值化相关函数(C_iZNCC)：

$C_{\mathrm{iZNCC}}=\frac{\underset{(x_i,y_i,z_i)\∈V}{\mathrm{\Σ}}[f(x_i,y_i,z_i)-{\stackrel{\‾}{f}}_m][g(x_i+u_0,y_i+v_0,z_i+w_0)-{\stackrel{\‾}{g}}_m]}{\sqrt{\underset{(x_i,y_i,z_i)\∈V}{\mathrm{\Σ}}{[f(x_i,y_i,z_i)-{\stackrel{\‾}{f}}_m]}^2\sqrt{\underset{(x_i,y_i,z_i)\∈V}{\mathrm{\Σ}}{[g(x_i+u_0,y_i+v_0,z_i+w_0)-{\stackrel{\‾}{g}}_m]}^2}}}---\left(1\right);$

其中， ${\stackrel{\‾}{f}}_m=\frac{\underset{(x_i,y_i,z_i)\∈V}{\mathrm{\Σ}}f(x_i,y_i,z_i)}{N};{\stackrel{\‾}{g}}_m=\frac{\underset{(x_i,y_i,z_i)\∈V}{\mathrm{\Σ}}g(x_i+u_0,y_i+v_0,z_i+w_0)}{N};$

式(1)中，V表示改进子区内部的所有有效体素点的集合，N是有效体素点的数目，f(x_i，y_i，z_i)和g(x_i+u₀，y_i+v₀，z_i+w₀)是变形前后在参考子区和目标子区中同一点的灰度值，和分别是参考子区和目标子区中有效体素点的灰度平均值；

步骤4，亚体素位移计算：采用最小二乘法计算亚体素位移，并利用如下的位移映射函数来描述同一计算点在变形前后的坐标对应关系：

$\left\{\begin{array}{c}{x}_i^{\′}=x_i+u_0+\mathrm{\Δu}+u_x\mathrm{\Δx}+u_y\mathrm{\Δy}+u_z\mathrm{\Δz}\\ y_i^{\′}=y_i+v_0+\mathrm{\Δv}+v_x\mathrm{\Δx}+v_y\mathrm{\Δy}+v_z\mathrm{\Δz}\\ z_i^{\′}=z_i+w_0+\mathrm{\Δw}+w_x\mathrm{\Δx}+w_y\mathrm{\Δy}+w_z\mathrm{\Δz}\end{array}\right.---\left(2\right);$

式(2)中(x_i，y_i，z_i)和(x_i’，y_i’，z_i’)是改进子区内同一点变形前后的坐标，Δx、Δy、Δz是该点和计算点(x₀，y₀，z₀)的距离，u₀、v₀、w₀和Δu、Δv、Δw分别为整体素位移和亚体素位移，u_x，u_y，u_z，v_x，v_y，v_z，w_x，w_y，w_z是位移梯度，同计算整体素位移时对传统子区的处理一样，利用改进的子区确定计算点的亚体素位移；

最终，理论的位移场、本发明方法得到的位移场(整体素位移加上亚体素位移)和传统DVC方法得到的位移场如图9所示；

步骤5，利用步骤3和步骤4获得的位移场，选取5×5×5的立方体计算应变场，对于所有的边界上和边界附近的点，选取只含有有效体素点的非规则窗口计算应变，即在

5×5×5的立方体内仅利用有效体素点按下式进行最小二乘拟合：

$\begin{array}{c}u(x,y,z)=a_0+a_{1}x+a_{2}y+a_{3}z,\\ v(x,y,z)=b_0+b_{1}x+b_{2}y+b_{3}z,\\ w(x,y,z)=c_0+c_{1}x+c_{2}y+c_{3}z,\end{array}---\left(3\right);$

式(3)中a₀，a₁，a₂，a₃，b₀，b₁，b₂，b₃，c₀，c₁，c₂，c₃是需要求解的系数，应变和系数的对应关系为：

$\begin{array}{c}{\mathrm{\ϵ}}_{\mathrm{xx}}=a_1,{\mathrm{\ϵ}}_{\mathrm{yy}}=b_2,{\mathrm{\ϵ}}_{\mathrm{zz}}=c_3,\\ {\mathrm{\ϵ}}_{\mathrm{xy}}=(b_1+a_2)/2,{\mathrm{\ϵ}}_{\mathrm{yz}}=(c_2+b_3)/2,{\mathrm{\ϵ}}_{\mathrm{xz}}=(a_3+c_1)/2.\end{array}---\left(4\right);$

理论的应变场、本发明方法得到的应变场和传统DVC方法得到的应变场如图10所示。

从图9和图10可明显看出，本发明方法能准确地得到理论写入的位移场和应变场，而传统DVC方法在边界附近存在巨大的计算错误问题，本发明方法有效地提高了边界附近计算点的位移场和应变场计算结果的准确性。

显然，上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无需也无法对所有的实施方式予以穷举，这些引伸出的变化或变动也处于本发明的保护范围之中。

Claims (4)

1.一种数字体积相关算法在边界问题上的处理方法，其特征在于，包括如下步骤：

步骤1，根据三维图像数据的灰度和边界形状确定计算区域，根据计算区域区分三维图像数据的有效体素点和无效体素点，位于计算区域内的点为有效体素点，位于计算区域外的点为无效体素点；

步骤2，整体素位移场的计算：当计算点位于边界上或者边界附近时，将围绕计算点传统子区内所有有效的体素点组成改进的子区，保持计算点在传统子区的位置不变，利用改进的子区确定计算点整体素的位移场；

步骤3，亚体素位移场的计算：将围绕计算点传统子区内所有有效的体素点组成改进的子区，保持计算点在传统子区的位置不变，利用改进的子区确定计算点亚体素的位移场；

步骤4，应变场计算：利用步骤2和步骤3获得的位移场，对于位于边界及边界附近的点，只选取计算区域内有效体素点的非规则窗口计算应变。

2.根据权利要求1所述数字体积相关算法在边界问题上的处理方法，其特征在于：所述边界包括内边界和外边界，所述内边界包括材料中的骨料、颗粒、纤维、气孔或裂纹或者为材料或结构的预制内边界，所述外边界包含规则形状样品的外边界和不规则形状样品的外边界。

3.根据权利要求1所述数字体积相关算法在边界问题上的处理方法，其特征在于：步骤1中，所述三维图像数据是指采用X射线断层照相、伽马射线断层照相、核磁共振断层照相、中子断层照相、聚焦离子束-电子束-双束显微镜或共聚焦显微镜进行三维成像得到的三维图像数据。

4.根据权利要求3所述数字体积相关算法在边界问题上的处理方法，其特征在于：所述X射线断层照相是指医用X射线断层照相、工业X射线断层照相、显微X射线断层照相、纳米X射线断层照相或同步辐射X射线断层照相中的任意一种。