CN103591904A - 一种两步三维傅里叶变换测量物体内部三维变形场的方法 - Google Patents

Abstract

本发明提供了一种利用两次三维傅里叶变换测量物体内部三维位移场及应变场的方法。利用微焦点计算机层析，微核磁共振或激光扫描共焦显微镜获得内部具有天然或人工制作的散斑结构的物体变形前后的三维图像。将变形前后的三维图像分别分割成若干相同尺寸的子体块，相对应的成对子块体逐一进行第一步三维傅里叶变换，在频谱域进行数值干涉得到合成谱，对合成谱进行第二步的三维傅里叶变换，得到含位移信息的扩展脉冲函数，通过求取峰值位置得到相对应子块体的位移矢量，由所有子块***移矢量获得物体内部三维变形场。通过插值获得位移的亚体素精度。由应变与变形间的关系方程，求出物体内部的应变场。

Description

一种两步三维傅里叶变换测量物体内部三维变形场的方法

技术领域

本发明涉及物体内部三维变形场的测量，属固体实验力学领域。

背景技术

数字散斑相关法或数字图像相关法(Digital Image Correlation，DIC)已广泛用于二维平面变形及应变的测量。数字体相关法(Digital Volume Correlation，DVC)是数字图像相关法向物体内部三维变形测量的拓展，该方法利用微焦点CT，微焦点核磁共振或激光扫描共焦显微镜获取得变形前后物体的三维立体图像，直接在变形后的体图像中搜索与变形前体图像中相对应的点，得出点的位移矢量，由于要在搜索区域中进行逐一搜索，致使计算时间消耗较大，尤其当物体产生大变形，搜索区域增大，耗费的时间将大为增加，计算效率降低。

为提高计算效率，本发明提供一种利用两步三维傅里叶变换测量物体内部三维变形场的方法。

发明内容

本发明的目的是提供一种高效地测量计算内部具有散斑结构物体的三维变形及应变场的方法。

本发明用于测量计算内部具有散斑结构物体的三维变形及应变场，物体可以是透明的或不透明的，散斑结构可以是物体内部的天然结构所形成的，如骨骼、木材、孔隙类岩石或含多种成份的岩石等，也可以是人为在待研究物体中加入某种在图像中易于辨识的微颗粒形成散斑结构。

本发明利用微焦点计算机层析(Micro focus computer tomography,Micro-CT)，微核磁共振成像(Micronuclear magnetic resonance imaging，Micro-MRI)，激光扫描共焦显微镜(Laser Scanning Confocal Microscope，LSCM)或其他可以获取反映物体内部结构的三维图像的方法获取具有散斑结构物体变形前后的三维立体图像，所得立体图像分别定义为参考体图像和变形体图像，将参考体图像与变形体图像分割成若干大小相等的子体块，参考体图像与变形体图像中相对应的一对子体块构成“斑对”，对构成“斑对”的两个子体块分别进行第一步三维傅里叶变换，得到相应的频谱函数，在频谱域中对所得到的频谱函数进行数值干涉，得到合成光谱函数，对该函数进行第二步傅里叶变换，得到一个峰值中心偏移的扩展脉冲函数，通过求取该扩展脉冲函数峰值的坐标，得到参考子体块的位移矢量。对参考体图像与变形体图像分割后所有相对应的“斑对”逐一进行上述两步三维傅里叶变换，就可以得到物体内部的变形场。通过应变与变形的关系，可以求解物体内部应变场。

本发明的具体步骤是：

步骤一：利用微焦点计算机层析(Micro focus computer tomography,Micro-CT)，微核磁共振成像(Micronuclear magnetic resonance imaging，Micro-MRI)，激光扫描共焦显微镜(Laser Scanning Confocal Microscope，LSCM)或其他可以获取反映物体内部结构的三维图像的方法获取具有散斑结构物体变形前后的三维立体图像，分别定义为参考体图像和变形体图像；

步骤二：将参考体图像与变形体图像分割成若干大小如32×32×32体素的子体块，相对应的参考子体块和变形子体块构成一个“斑对”；当子体块尺寸较小时，可以认为子体块的内各点的位移是相同的；

本步骤中子块体的尺寸可以是4×4×4体素，8×8×8体素，16×16×16体素，32×32×32体素，64×64×64体素，……，根据物体内部散斑结构的大小，适当选取子块体的大小，但每个子块体中散斑结构的数量应大于9。

本步骤中体图像的分割相邻子体块之间可以有二分之一子块体尺寸的重叠，提高计算的准确性。

步骤三：令h₁(x，y，z)和h₂(x，y，z)分别表示一个“斑对”中的参体子体块和变形子体块，h₁(x，y，z)与h₂(x，y，z)具有如下关系：

h₂(x，y，z)=h₁[x-u(x，y，z)，y-v(x，y，z)，z-w(x，y，z)]   (1)

式中，u，v，w分别为变形后参考子体块在x，y，z三个方向上的位移，对参考子体块和变形子体块分别进行第一步三维傅里叶变换，则有：

$\begin{array}{c}{H}_1(f_x,f_y,f_z)\\ =\underset{\mathrm{\Δ}}{\∫\∫\∫}{h}_1(x,y,z)\exp [-j2\mathrm{\π}({\mathrm{xf}}_x+{\mathrm{yf}}_y+{\mathrm{zf}}_z)]\mathrm{dxdydz}\\ =\left|H_1(f_x,f_y,f_z)\right|\exp \left[{\mathrm{j\φ}}_1(f_x,f_y,f_z)\right]\\ H_2(f_x,f_y,f_z)\\ =\underset{\mathrm{\Δ}}{\∫\∫\∫}{h}_1(x-u,y-v,z-w)\exp [-j2\mathrm{\π}({\mathrm{xf}}_x+{\mathrm{yf}}_y+{\mathrm{zf}}_z)]\mathrm{dxdydz}\\ =\left|H_2(f_x,f_y,f_z)\right|\exp \left[{\mathrm{j\θ}}_2(f_x,f_y,f_z)\right]\\ =\left|H_1(f_x,f_y,f_z)\right|\exp \left\{j\right[{\mathrm{\φ}}_1(f_x,f_y,f_z)-2\mathrm{\π}({\mathrm{uf}}_x+{\mathrm{vf}}_y+{\mathrm{wf}}_z)\left]\right\}\end{array}---\left(2\right)$

式是Δ表示子体块区域，(f_x，f_y，f_z)为第一次三维傅里叶变换后的频谱域坐标，H₁(f_x，f_y，f_z)为h₁(x，y，z)傅里叶变换后的函数，|H₁(f_x，f_y，f_z)|和φ₁(f_x，f_y，f_z)分别为H₁(f_x，f_y，f_z)的幅值与相位，H₂(f_x，g_y，f_z)为h₂(x，y，z)傅里叶变换后的函数，|H₂(f_x，f_y，f_z)|和φ₂(f_x，f_y，f_z)分别为H₂(f_x，f_y，f_z)的幅值与相位；

步骤四：对步骤三中得到了H₁(f_x，f_y，f_z)和H₂(f_x，f_y，f_z)通过下式进行数值干涉生成新的合成光谱

$F(f_x,f_y,f_z)=\frac{H_1(f_x,f_y,f_z)H_2^{*}(f_x,f_y,f_z)}{{\left|H_1(f_x,f_y,f_z)H_2(f_x,f_y,f_z)\right|}^{1-\mathrm{\α}}}---\left(3\right)$

也可以写成

F(f_x，f_y，f_z)=|H₁(fx，f_y，f_z)H₂(f_x，f_y，f_z)|^αexp{j[φ₁(f_x，f_y，f_z)-φ₂(f_x，f_y，f_z)]}

                                                                                  (4)

式中0≤α≤1。

步骤五：对步骤四中所得的合成光谱进行第二步傅里叶变换。

将式(2)代入式(3)中可得

F(f_x，f_y，f_z)≈|H₁(f_x，f_y，f_z)|^2αexp[j2π(uf_x+vf_y+wf_z)]   (5)

$\begin{array}{c}G(\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})=\underset{{\mathrm{\Δ}}_f}{\∫\∫\∫}F(f_x,f_y,f_z)\exp [-j2\mathrm{\π}(f_x\mathrm{\ξ}+f_y\mathrm{\η}+f_z\mathrm{\ζ})]{\mathrm{df}}_x{\mathrm{df}}_y{\mathrm{df}}_z\\ \≈\underset{{\mathrm{\Δ}}_f}{\∫\∫\∫}{\left|H_1(f_x,f_y,f_z)\right|}^{2\mathrm{\α}}\exp \{-j2\mathrm{\π}[f_x(\mathrm{\ξ}-u)+f_y(\mathrm{\η}-v)+f_z(\mathrm{\ζ}-w)\left]\right\}{\mathrm{df}}_x{\mathrm{df}}_y{\mathrm{df}}_z\\ =\stackrel{\‾}{G}(\mathrm{\ξ}-u,\mathrm{\η}-v,\mathrm{\ζ}-w)\end{array}---\left(6\right)$

式中Δ_f表示子体块经第一步三维傅里叶变换后的频谱区域，(ξ，η，ζ)表示第二步三维傅里叶变换后的频谱域的坐标。

是一个峰值点在(u，v，w)处的扩展脉冲函数。

a的取值会影响

扩展脉冲函数峰值的锐利程度及峰值信号对噪声信号的抗干扰性，当a=0.5时，扩展脉冲函数峰值具有较好的锐利程度和抗干扰性。在本发明中，a取值为0.5。

步骤六：搜索求取步骤五中扩展脉冲函数峰值点的坐标值(u，v，w)，从而得到参考子体块的位移矢量；

步骤七：在数字图像中，步骤六中所求得的三个方向的位移均为整体素，为了提高测量精度，通过插值获得位移矢量的亚体素精度。在

扩展脉冲函数中，围绕着步骤六中所求得的峰值点(其坐标为整体素)，选取3×3×3或是5×5×5的矩阵，可采用已有的三维线性(trilinear)，三维多义线(tri-spline)或三维三次多项式插值(tri-cubic)等插值方法来求取亚体素精度的位移；

步骤八：对参考体图像与变形体图像分割成的若干子体块形成的“斑对”，逐一进行上述步骤三至步骤七的操作，最终获得整个物体内部的变形场。

步骤九：计算物体内部的应变场，由步骤八得到的物体内部变形场，通过应变与变形的关系方程，可以求解物体内部的应变场。本发明中采用逐点最小二乘法来求解应变。由步骤八所得物体内部的变形场，以待求点为中心，围绕其选取(2N+1)×(2N+1)×(2N+1)立方区域，当所选区域足够小时，其内部各点的位移可以看作是坐标的线性函数，因此有：

u(x，y，z)=a₀+a₁x+a₂y+a₃z

v(x，y，z)＝b₀+b₁x+b₂y+b₃z

w(x，y，z)=c₀+c₁x+c₂y+c₃z

式中，x,y,z取值范围为[-N,N]，且为整数。u(x，y，z)，v(x，y，z)和w(x，y，z)为步骤八中所求得的物体内部变形。a_{i＝0，1，2，3}，b_{i＝0，1，2，3}和c_{i＝0，1，2，3}是待定的多项式系数，可以通过已有的最小二乘法或多元拟合方法来求取。这些系数确定后，待求点的Cauchy应变值为：

${\mathrm{\ϵ}}_x=\frac{\∂u}{\∂x}=a_{1,}{\mathrm{\ϵ}}_{\mathrm{xy}}=\frac{1}{2}(\frac{\∂v}{\∂x}+\frac{\∂u}{\∂y})=\frac{1}{2}(b_1+a_2)$

${\mathrm{\ϵ}}_y=\frac{\∂v}{\∂y}=b_{1,}{\mathrm{\ϵ}}_{\mathrm{yz}}=\frac{1}{2}(\frac{\∂w}{\∂y}+\frac{\∂v}{\∂z})=\frac{1}{2}(c_2+b_3)$

${\mathrm{\ϵ}}_z=\frac{\∂w}{\∂z}=c_{1,}{\mathrm{\ϵ}}_{\mathrm{xz}}=\frac{1}{2}(\frac{\∂u}{\∂z}+\frac{\∂w}{\∂x})=\frac{1}{2}(a_3+c_1)$

附图说明

图1为本发明的原理示意图；

其中1.参考体图像，2.变形体图像，3.分割后的参考子体图像矩阵，4.分割后的变形子体图像矩阵，5.参考子体块，6.变形子体块，7.扩展脉冲函数的等值分布

图2为中砂岩变形后的内部三维位移分量，x方向的位移u场截面图；

图3为中砂岩变形后的内部三维位移分量，y方向的位移v场截面图；

图4为中砂岩变形后的内部三维位移分量，z方向的位移w场截面图；

具体实施方式

下面结合附图与实施例对本发明作进一步说明

步骤一：利用微焦点计算机层析(Micro focus computer tomography,Micro-CT)扫描获得中砂岩试件变形前后的三维立体图像，分别定义为参考体图像1和变形体图像2；

步骤二：将参考体图像1和变形体图像2分别分割成若干大小如32×32×32体素的子体块矩阵3和4，相对应的每个子块体如5和6构成一个“斑对”，令h₁(x，y，z)和h₂(x，y，z)分别表示变形前参考子体块5及变形后相应的子体块6，h₁(x，y，z)与h₂(x，y，z)具有如下关系：

h₂(x，y，z)=h₁[x-u(x，y，z)，y-v(x，y，z)，z-w(x，y，z)]   (1)

式中，u，v，w分别为子块体5在x，y，z三个方向上的位移，

步骤三：将构成“斑对”的变形前后的子块体5和6分别进行第一步三维傅里叶变换，则有：

式中

表示傅里叶变换，H₁，H₂分别为h₁和h₂的傅里叶函数，f_x、f_y、f_z分别是第一步三维傅里叶变换后频谱空间的坐标。

步骤四：在频谱域中，对上述两个三维散斑体频谱H1和H2进行数值干涉，得到合成谱如下：

$\begin{array}{c}F(f_x,f_y,f_z)=\frac{H_1(f_x,f_y,f_z)H_2^{*}(f_x,f_y,f_z)}{\sqrt{\left|H_1(f_x,f_y,f_z)H_2(f_x,f_y,f_z)\right|}}\\ =\exp \left\{j\right[{\mathrm{\φ}}_1(f_x,f_y,f_z)-{\mathrm{\φ}}_2(f_x,f_y,f_z)\left]\right\}\end{array}---\left(3\right)$

式中*表示共轭，Φ₁(f_x，f_y，f_z)和Φ₂(f_x，f_y，f_z)分别为H₁(f_x，f_y，f_z)和H₂(f_x，f_y，f_z)的相位角，由相位差与光程差之间的关系则有：

φ₁(f_x，f_y，f_z)-φ₂(f_x，f_y，f_z)=2π(uf_x+vf_y+wf_z)   (4)

步骤五：对数值干涉得到的合成谱F(f_x，f_y，f_z)进行第二次三维傅里叶变换，有：

上式表示中心峰值点位于(u,v,w)的扩展脉冲函数，其三维分布如图1中7所示。

步骤六：通过搜索步骤五所得扩展脉冲函数的峰值点坐标，得到参考子体块5的位移矢量；

步骤七：为了提高测量精度，通过插值获得位移矢量的亚体素精度。在

扩展脉冲函数中，围绕着步骤六中所求得的峰值点坐标(为整体素)，选取3×3×3或是5×5×5的矩阵，可采用已有的三维线性(tri-linear)，三维多义线(tri-spline)或三维三次多项式插值(tri-cubic)等方法来求取亚体素精度的位移；

步骤八：对所有子块体形成的“斑对”进行上述步骤三到步骤七的运算，最终获得整个物体内部的变形场。图2为中砂岩变形后的内部三维位移分量，x方向的位移u场截面图；图3为中砂岩变形后的内部三维位移分量，y方向的位移v场截面图；图4为中砂岩变形后的内部三维位移分量，z方向的位移w场截面图；

步骤九：计算物体内部的应变场，由步骤八得到的物体内部的变形场，通过应变与变形的关系方程，可以求解物体内部的应变场。

Claims (10)

1.一种利用两步三维傅里叶变换测量物体内部三维变形的方法，其特征是：利用计算机层析技术，核磁共振成像技术或激光扫描共焦显微镜等方法获取具有散斑结构特征物体在变形前后的三维立体图像，将物体变形前后的三维立体图像分割成一系列子体块，相对应的子体块构成“斑对”，分别进行第一步三维傅里叶变换，在频谱域中进行数值干涉，对所得合成光谱进行第二次三维傅里叶变换，得到该子体块的位移矢量，对所有相对应的子体块进行上述两步三维傅里叶变换，就可以求出物体内各子体块的位移矢量。根据应变与变形的关系方程，可以获得物体内部应变场。 

2.根据权利要求1所述的一种两步三维傅里叶变换测量物体内部三维变形场的方法，该物体可以是透明的或是非透明的，但其内部应具有天然或人为制作的结构作为散斑结构。 

3.根据权利要求1所述的一种两步三维傅里叶变换测量物体内部三维变形场的方法，物体的三维立体图像可以由微焦点计算机层析(Micro focus computer tomography,Micro-CT)，微核磁共振成像(Micro nuclear magnetic resonance imaging，Micro-MRI)或激光扫描共焦显微镜(Laser Scanning Confocal Microscope，LSCM)或其他可以获取反映物体内部结构的三维图像的方法获取。 

4.根据权利要求1所述的一种两步三维傅里叶变换测量物体内部三维变形场的方法，在将参考体图像与变形体图像进行分割时，子块体的尺寸可以是4×4×4体素，8×8×8体素，16×16×16体素，32×32×32体素，64×64×64体素……，根据物体内部散斑结构的大小，适当选取子块体的大小，但每个子块体中散斑结构的数量应大于9。在将参考体图像与变形体图像进行分割时，相邻子体块之间可以有二分之一子块体尺寸的重叠，提高计算的准确性。 

5.根据权利要求1所述的一种两步三维傅里叶变换测量物体内部三维变形场的方法，参考体图像与变形体图像分割后相对应的子体块构成“斑对”，分别进行第一步三维傅里叶变换。 

令h₁(x，y，z)和h₂(x，y，z)分别表示变形前参考子体块及其变形后相应变形子体块，h₁(x，y，z)与h₂(x，y，z)具有如下关系： 

h₂(x，y，z)=h₁[x-u(x，y，z)，y-v(x，y，z)，z-w(x，y，z)] 

式中，u，v，w分别为参考子块体在x，y，z三个方向上的位移， 

将变形前后的子块体进行第一步三维傅里叶变换，有： 

式中

表示傅里叶变换，H₁，H₂分别为h₁和h₂的傅里叶函数，f_x、f_y、f_z分别是第一步傅里叶变换后频谱空间的坐标。 

6.根据权利要求1所述的一种两步三维傅里叶变换测量物体内部三维变形场的方法，在频谱域中，对两个三维散斑体频谱H₁和H₂进行数值干涉，得到合成光谱如下： 

式中*表示共轭，Φ₁(f_x，f_y，f_z)和Φ₂(f_x，f_y，f_z)分别为H₁(f_x，f_y，f_z)和H₂(f_x，f_y，f_z)的相位角，由相位差与光 

程差之间的关系可得： 

φ₁(f_x，f_y，f_z)-φ₂(f_x，f_y，f_z)=2π(uf_x+vf_y+wf_z)。 

7.根据权利要求1所述的一种两步三维傅里叶变换测量物体内部三维变形场的方法，对数值干涉得到的合成光谱F(f_x，f_y，f_z)进行第二次三维傅里叶变换，有： 

上式表示中心峰值点位于点(u，v，w)的扩展脉冲函数。通过搜索得到扩展脉冲函数的峰值点坐标，从而得到参考子体块的位移矢量。 

8.根据权利要求1所述的一种两步三维傅里叶变换测量物体内部三维变形场的方法，为了提高测量精度，通过插值获得位移矢量的亚体素精度。在

扩展脉冲函数中，围绕所求得的峰值点(其坐标为整体素)，选取3×3×3或是5×5×5的矩阵，可采用已有的三维线性(trilinear)，三维多义线(tri-spline)或三维三次多项式插值(tri-cubic)等插值方法来求取亚体素精度的位移。 

9.根据权利要求1所述的一种两步三维傅里叶变换测量物体内部三维变形场的方法，对所有的子块体形成的“斑对”进行两步三维傅里叶变换的运算，最终获得整个物体内部的变形场。 

10.根据权利要求1所述的一种两步三维傅里叶变换测量物体内部三维变形场的方法，通过应变与变形之间的关系方程，由得到的物体内部变形场可以求解物体内部的应变场。由所求得物体内部的变形场，以待求点为中心，围绕其选取(2N+1)×(2N+1)×(2N+1)立方区域，当所选区域足够小时，其内部各点的位移可以看作是坐标的线性函数，因此有： 

u(x，y，z)=a₀+a₁x+a₂y+a₃z 

v(x，y，z)=b₀+b₁x+b₂y+b₃z 

w(x，y，z)=c₀+c₁x+c₂y+c₃z 

式中，x，y，z取值范围为[-N,N]，且为整数。u(x，y，z)，v(x，y，z)和w(x，y，z)是由两步三维傅里叶变换所求得的物体内部变形。a_{i＝0，1，2，3}，b_{i＝0，1，2，3}和c_{i＝0，1，2，3}是待定的多项式系数，可以通过已有的最小二乘法或多元拟合方法来求取。这些系数确定后，待求点的Cauchy应变值为： 

。

Images