CN103776559A - 断层成像激光剪切干涉三维温度测量装置及温度测量方法 - Google Patents

Abstract

一种断层成像激光剪切干涉三维温度测量装置及温度测量方法，属于工业自动化仪表领域。它由HeNe激光器(1)所发出的激光束经过两个透镜（2）和（3）等光学扩束后，光线经过第一个消色差透镜（6）使之成为平行光。平行光在透过温度测试空腔（7）后被第二个消色差透镜（8）会聚，并经过在会聚焦点位置的Wollaston棱镜（9）后成为干涉图，由CCD摄像仪（10）记录干涉图像。根据测试容器的特点，利用角度旋转装置（11）选取4个以上不同角度拍摄干涉图像。然后根据代数重建算法对不同角度拍摄的干涉图像进行光学断层成像技术处理，从而获得测量容器内的三维温度场。该测温装置不仅具有无干涉流体、系统简单、成本低等特点，而且具有高精度三维测量温度场的显著优点。

Description

断层成像激光剪切干涉三维温度测量装置及温度测量方法

所属技术领域

本发明涉及一种新型三维温度测量技术，属于工业自动化仪表领域。

背景技术

空腔内液体三维温度无干涉测量技术一直是工业界最需要解决的测量手段。在流体力学尤其是对流现象的研究中，人们通常使用光学干涉的方法观察和测量流体的温度分布。通过光学干涉可以得到空腔内流场的无限宽干涉条纹，这些条纹携带了光束路径（z方向）上的折射率在某一方向（可能是x或y方向）的变化率叠加的信息，这些信息是利用代数重建算法（ART-Algebraic ReconstructionTechnique）进行温度场重建所需要的投影数据，根据这些数据应用ART可以重建空腔内对流流场的三维温度分布。

折射率是介质本身的属性并和入射光的波长相关，其值等于光线入射角的正弦值与折射角的正弦值之比：

比值n₁₂称为介质2相对于介质1的折射率，这就是人们常说的斯涅耳定律。一种介质相对于真空的折射率，称作该介质的绝对折射率，也称为折射率。对于大多数的流体而言（极少数除外），流体的折射率与其密度具有直接的关系。在非等温流体系统中，流体的密度随着温度的变化而减小。通常情况下，流体密度的变化与温度呈简单的线性关系。由于非等温流体系统内流体折射率的非均匀性，当激光穿过这一流体系统就会被干扰，与未受干扰的光线相比会出现以下情形：

(1)光线偏离原传播方向；

(2)光线的相位发生变化。

利用这两种现象可以开展流动的可视化研究和流体密度的定量测量。定量测量必须知道密度和折射率之间精确的变化关系。

本世纪二、三十年代起干涉测温法就得到了广泛应用,它通过干涉法测量由温度变化引起的热膨胀和折射率变化来测定样品温度,可以非接触、无损地测量样品温度,测量响应速度快,但它只能测量固体表面的温度。

激光全息干涉法温度测量装置也是利用干涉法测量由温度变化引起的折射率变化的干涉图来进行温度测量。其原理是测量激光通过测量容器的全息图并与不经过测量容器的参考光源的干涉图相比较，得出测量容器一个二维平面上的温度场。其结构如下图3（见杨天亮等，全息实时干涉法测量二维温度场的实验研究，[J].西安交通大学学报,1987(5).）。

由于该方法是通过比较激光通过测量容器的全息图与不经过测量容器的参考光源的干涉图而得出温度场的，因此其缺点是光路中需要分出一参考光源，使系统复杂，且只能测量二维的温度场。轴对称和三维温度场虽然可测，但非常复杂。

发明内容

本发明目的在于提供一种不仅具有无干涉流体、系统简单、成本低等特点，而且具有高精度三维测量温度场的显著优点的断层成像激光剪切干涉三维温度测量装置及其温度测量方法。

一种断层成像激光剪切干涉三维温度测量装置，其特征在于：包括角度旋转装置和安装于角度旋转装置上的光学器件；所述的光学器件依次由激光器、光学扩束装置、第一消色差透镜、测试空腔、第二消色差透镜、Wollaston棱镜和CCD摄像仪组成；所述的光学扩束装置依次由第一透镜第一孔径光阑、第二透镜、λ/4波片、第二孔径光阑组成。

利用所述断层成像激光剪切干涉三维温度测量装置的温度测量方法，其特征在于包括以下过程：

一、预处理工作：

1-1、所述的断层成像激光剪切干涉三维温度测量装置在光线经过测量容器后会聚焦点位置上安装有Wollaston棱镜（9）；然后调整光；路中各元件的位置，确保在未加测试空腔（7）时，CCD摄像仪（10）可以拍摄到以有限宽方式显示的预加干涉条纹；

1-2、选取4个以上的角度对容器进行拍摄并获得清晰的干涉图，其中至少有1个角度所在的方向与其他角度的方向相互垂直；

1-3、对所获得清晰的干涉图做二值化处理，把图像转化成只含有黑白两种颜色的图片；用Matlab读各个角度的二值图像，得到每个像素的颜色值，此时各个角度干涉图的有效像素的数量都是48×291，每幅二值图像的颜色值分别存入一个颜色值数组；

1-4、按从下到上从左到右的顺序处理数据原则，计算每幅图像的每个像素点的像素条纹级数值，其中设每一列最下面的像素点的像素条纹级数值为0，同一列中若当前像素点的像素颜色与前一像素点的像素颜色产生黑到白或白到黑的变化则当前像素点的像素条纹级数值为前一像素点的像素条纹级数值+1，若当前像素点的像素颜色与前一像素点的像素颜色未产生变化则当前像素点的像素条纹级数值等于前一像素点的像素条纹级数值，每幅图像的像素条纹级数值分别存入一个条纹级数值数组；

1-5、对每幅二值图像所得的条纹级数值数组进行插值化处理，得到条纹级插值化数值数组；插值化处理的具体方法为：按从下到上的顺序处理数据原则，从每一列的第一个像素开始，寻找与其像素条纹级数值相差为1的像素；然后计算这两个像素之间间隔的像素数k，并计算得到△k=1/(k+1)；对于两个像素之间的间隔像素，用△k依次乘以间隔像素序号并加到其像素条纹级数值上，其中间隔像素序号从1开始计数；然后从这一列的下一个像素继续开始上述步骤，直到该列的第48个像素；所有的列都计算完毕后，得到条纹级插值化数值数组；

1-6、将重建区域划分为一个n×n的正方形网格，也称为像素，N＝n×n为像素总数；在每个像素内f(x,y)是常量，f_j代表第j个像素内的常数值；光束被定义为有一固定宽度的直线，把光束一一编号，编号用p的下标i表示；w_ji是第i个光束通过第j个像素的面积占整个第j个像素面积的比例，称做权重因子；

1-7、用每条光线经过的所有像素各自的权重因子w_ji和f_j的乘积表示每个像素对每条经过其本身的光线的贡献，一条光线经过的所有像素的这些值的和称作一个投影，这样，每个像素的函数值f_j和每条光线的投影值p_i的关系可以用下式表示：

$p_i=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_{\mathrm{ij}}{f}_j,i=\mathrm{1,2,3},...M$

其中，M为投影总数。

二、迭代过程

对于上面获得的投影值和权重因子之间的方程组，利用迭代法求解；

先要把方程(1)改写成以下形式：

$\begin{array}{c}{w}_{11}{f}_1+w_{12}{f}_2+w_{13}{f}_3+...+w_{1N}{f}_N=p_1\\ w_{21}{f}_1+w_{22}{f}_2+w_{23}{f}_3+...+w_{21N}{f}_N=p_2\\ .\\ .\\ .\\ w_{M1}{f}_1+w_{M2}{f}_2+w_{N3}{f}_3+...+w_{\mathrm{MN}}{f}_N=p_M\end{array}---\left(2\right)$

N个像素代表了一个图像的N个自由度，所以，(f₁,f₂,…f_N)可以被认为表示N维空间的一个点，在方程组(2)中的每一个方程，都可以被看作是这个空间的一个超平面。当这个方程组存在唯一解，就意味着这M个超平面相交于一点，该点即为唯一解表示的点；如果方程组存在唯一解，迭代就会收敛于唯一解表示的点；因此，通过迭代方法可以求得f_j的值，f_j即第j个像素的折射率，求解每个像素的折射率；

三、温度场的反演计算

利用上述求解获得的折射率的分布，通过折射率与流体密度之间的关系式n-1＝Kρ获得流体密度场分布数据，其中K为Gladstone-Dale常数，它具有1/ρ的量纲，并且与气体的特性及所用的光的频率或者波长有关，(n-1)称之为折射率；进一步通过流体密度与温度之间的线性变化关系式ρ＝f(T)，反演推导出温度场的分布；其中ρ＝f(T)这一公式根据具体的流体而定。

本发明主要在于提出一种利用激光剪切干涉的方法测量不同角度下因流体密度变化而引起流体折射率变化的干涉图像，通过代数重建算法对不同角度拍摄的干涉图像进行光学断层成像技术处理，从而获得测量容器内三维温度场。因此，该新型激光干涉测温装置不仅具有无干涉流体、系统简单、成本低等特点，而且具有高精度三维测量温度场的显著优点。

附图说明

图1是一种新型断层成像激光剪切干涉三维温度测量装置；

图2是光路系统安装角度底座图；

图3是二维激光全息测温度场的装置简图；

图4是光线投影示意图；

图中标号名称：1.HeNe激光器,2.第一透镜,3.第一孔径光阑,4.第二透镜,5.第二孔径光阑,6.第一消色差透镜,7.测试空腔,8.第二消色差透镜,9.Wollaston棱镜,10.CCD摄像仪，11.支架,12.偏振片,13.角度刻度盘指针,14.角度刻度盘底座，15.λ/4波片。

f₁、f₂一个二维平面内的两个常数值，H假定的起始点位置，G H的第一次投影位置，p₁、p₂光线1和2上的投影值，w_ji第i个光束通过第j个像素的面积占整个第j个像素面积的比例，称做权重因子；i光束的编号（1，2）,j像素的编号（1，2）H点的位置矢量，

G点的位置矢量，G点在光线2上投影的位置矢量。

具体实施方式

下面结合图1进一步说明本发明的新型断层成像激光剪切干涉三维温度测量装置。由HeNe激光器(1)所发出的激光束经过两个透镜（2和3）及两个孔径光阑（4和5）组成光学扩束装置，然后经过消色差透镜6形成单色平行光，进入测试空腔。经过测试空腔内流体的折射而透过的光携带有光线经过路程上流体折射率的信息致之和，再经过消色差透镜8汇聚到其焦点上的Wollaston棱镜9形成一干涉图象，由CCD摄像仪10记录其干涉图像。

通过角度旋转装置11旋转安装在其上的光路测量系统，测量任意一个角度下的激光通过容器后的干涉图像。这个角度可以精确地按照需要任意选取。

根据CCD摄像仪拍摄的不同角度上的干涉图像，根据代数重建算法软件进行光学断层成像技术处理，从而获得测量容器内三维温度场。

上述光学元件按照测量容器的大小选择。Wollaston棱镜在光路的位置对最后拍摄到条纹质量有决定性的影响，理论上，Wollaston棱镜要严格放置在第二个消色差透镜的焦点处。

新型断层成像激光剪切干涉三维温度测量方法实施的具体步骤如下：

一、预处理工作：

1-1、调整光路中各元件到可以拍摄到未加测试段的以有限宽方式显示的预加干涉条纹。

1-2、根据算法和实验容器的特点，选取4个以上角度对容器进行拍摄并获得清晰的干涉图片，其中至少有一个角度所在的方向与其他角度的方向相互垂直。

1-3、对干涉图做二值化处理，把图像转化成只含有黑白两种颜色的图片。

用Matlab读二值图像，得到每个像素的颜色值，此时各个角度干涉图的有效像素的数量都是48×291，每幅图像的颜色值分别存入一个数组。

1-4、接下来的过程由C++程序实现：程序按角度的顺序读入上面存有图像颜色的数组，读入之后，按照角度的顺序做以下工作：每一角度的图分为291列48行，按从下到上从左到右的顺序处理数据，建立一个新的数组，条纹级数值数组，数组大小为48×291，用来存储干涉图每个像素点的条纹级数值，设每一列最下面的像素条纹级数值为0，然后一个像素一个像素地向上数过去，见到相邻两个像素颜色变化（黑变白或和白变黑）就在变化后的像素的条纹级数值上加1，然后继续往上数，见到颜色变化再加1，这样直到本列的末尾（也就是第48个像素）。按这样的规则数完一列像素，最后存储在新建立数组中的条纹级数应该是一系列离散的量，如：0000011112222233333333344444444……，条纹数量越多，这个数值越大。

1-5、用上述得到的离散的量做图并不能得到平稳过渡的折射率变化量/温度等值线，对新建立的条纹级数值数组中的条纹级数进行插值，尽量得到平稳过渡的比较连续的折射率变化量/温度等值线。插值的具体方法为：在每一列的第一个像素开始，寻找相差为1的两个像素，然后计算这两个像素之间间隔的像素数，用1除以这个像素数然后依次乘以序号加到两个像素中间的各个像素中，如上面举例的条纹级数：00000111122222233333333334444444……中,从第一个像素开始，数到第六个像素找到差为1的像素，用1除以某个数，这个数为它们之间的像素数+1，用1除以这个数得到0.2，然后在第二个像素上加0.2，第三个像素上加0.2×2，第四个加0.2×3，第五个加0.2×4，加后的结果，条纹级数值变为0，0.2，0.4，0.6，0.8，1，1，1，1，2，2，2，2，2，3，……这样就得到比刚开始的条纹级数过渡更加平滑的一系列值。然后从这一列的下一个像素继续开始这个步骤，直到最后的第48个像素。所有的列都计算完毕后，这个角度的干涉图就计算完毕，最后得到的条纹级数的值有两个用途，一是用来做这个角度的折射率变化量/温度等值线，二是用于后面计算的迭代过程，作为ART算法中的投影数据。

1-6、将重建区域划分为一个n×n的正方形网格，也称为像素，在每个像素内f(x,y)是常量。f_j代表第j个像素内的常数值，N＝n×n为像素总数。光束被定义为有固定宽度的直线，把光束一一编号，编号用p的下标i表示。w_ji是第i个光束通过第j个像素的面积占整个第j个像素面积的比例，称做权重因子；

1-7、用每条光线经过的所有像素各自的权重因子w_ji和f_j的乘积表示每个像素对每条经过其本身的光线的贡献，一条光线经过的所有像素的这些值的和称作一个投影，这样，每个像素的函数值f_j和每条光线的投影值p_i的关系可以用下式表示：

$p_i=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_{\mathrm{ij}}{f}_j,i=\mathrm{1,2,3},...M$

其中，M为投影总数。

二、迭代过程

对于上面获得的投影值和权重因子之间的方程组，如果M，N都是很小的数，我们可以应用矩阵理论的方法来解这个方程组。但是，实际上N通常是65000(256×256像素)这个数量级，同时M也是这个数量级的量。这样对于方程(1)来说，矩阵的规模是65000×65000，这个规模排除了任何直接使用矩阵理论求解方程组的可能。当投影数据中存在噪声并且M＜N，即使N非常小，也不能使用矩阵转换的方法求解方程组。

迭代法是解决上述求解大型稀疏矩阵难题的可行方法。迭代法基于Kaczmarz^[40]首先提出的投影方法，Tanabe^[41]后来对这一方法进行了详细的阐述。为了解释迭代法的计算步骤，先要把方程(1)改写成以下形式：

$\begin{array}{c}{w}_{11}{f}_1+w_{12}{f}_2+w_{13}{f}_3+...+w_{1N}{f}_N=p_1\\ w_{21}{f}_1+w_{22}{f}_2+w_{23}{f}_3+...+w_{21N}{f}_N=p_2\\ .\\ .\\ .\\ w_{M1}{f}_1+w_{M2}{f}_2+w_{N3}{f}_3+...+w_{\mathrm{MN}}{f}_N=p_M\end{array}---\left(2\right)$

N个像素代表了一个图像的N个自由度，所以，(f₁,f₂,…f_N)可以被认为表示N维空间的一个点，在方程组(2)中的每一个方程，都可以被看作是这个空间的一个超平面。当这个方程组存在唯一解，就意味着这M个超平面相交于一点(唯一解表示的点)。为了说明，假设一个只有两个变量f₁，f₂的方程组如下：

$\begin{array}{c}{w}_{11}{f}_1+w_{12}{f}_2=p_1\\ w_{12}{f}_1+w_{22}{f}_2=p_2\end{array}---\left(3\right)$

这个方程组的计算步骤为：先假设一个任意的初始值H(见图4)，把这个初始值投影到第一条光线1，然后投影到第二条光线2，然后再向第一条光线投影，然后再向第二条，如此继续下去。如果方程组存在唯一解，迭代就会收敛于唯一解表示的点。

在编写此解法的计算程序时，开始要假设一个初始值，标记为

用N维空间的向量表示就是

在大多数情况下，

初始值都被设置为零，然后这个初始值被投影到在方程组(2)中p₁所代表的超平面上得到一组新值

接着

被投影到方程组(2)中p₂所代表的超平面上得到

以此方式继续下去。

投影到第i个方程表示的超平面得到的过程可以用以下方程描述：

在方程(4)中，

为与自身的点积。

当一条光线的N与ind[]数组都计算完毕后，就可以带入迭代公式来计算Image[v][ind[i+1]]，迭代语句为：

for(i=0;i<N;i++)

Image[v][ind[i+1]]+=(data[g-2][t-1][v]-Q)/N;

式中data[g-2][t-1][v]为前面计算好的投影数据。由这个迭代公式可以看出，每次对图像向量进行修正只是对那些与投影射线相交的像素进行修正，而没有与投影射线相交的像素不处理。

由于投影数据中噪声的原因，迭代次数不能过多，但至今有关最优迭代次数并没有明确的说法，所以本文采用4次迭代。

三、温度场的反演计算

3-1、流体折射率与密度之间的关系

定量测量的一个必要条件是精确了解折射率与流体密度之间的关系。在德国学者W.Merzkirch的著作《流动显示》中，推导了一个用流体的分子常数和性质以及所用光束的频率υ表示的流体折射率n与流体密度ρ之间的关系式：

式中

为Loschmidt数，e为单个电子的电荷量，m_e为单个电子的质量，M是流体分子量，f_i为假定的每个分子中若干个可以位移的电子的振子强度，υi为这些电子的共振频率，这个式子成立需要这样的假设：流体的共振频率与所用光的频率相差很远。这个式子被称为Clausius-Mosotti关系式，适用于气体和液体。说明流体折射率和流体密度之间存在着非常精确的对应关系。大多数气体的折射率很接近于1，所以对气体而言，可令(n²-1)≌2(n-1)，(n²+2)≌3,由此将Clausius-Mosotti关系式加以简化，得到Gladstone-Dale关系式：

其简化形式为：

n-1＝Kρ   (7)此式定义了Gladstone-Dale常数K，它具有1/ρ的量纲，并且与气体的某些特性及所用的光的频率或者波长有关。差值(n-1)也称之为折射率，液体的(n-1)比气体的要大几个数量级，因此，液体中极小的密度差异在光学方法中会产生很强的信号，获得有关的分子方面的数据，折射率和水(或任何其他液体)密度之间的关系可以利用上面的Clausius-Mosotti关系式得到，且这一关系需要通过标定来确定。

3-2、流体密度与温度之间的关系

通常而言，流体密度与温度之间成线性变化的关系ρ＝f(T)（这一公式根据具体的流体而定），从利用这一关系式反演推导出温度场的分布。

Claims (2)

1.一种断层成像激光剪切干涉三维温度测量装置，其特征在于：包括角度旋转装置（11）和安装于角度旋转装置（11）上的光学器件；所述的光学器件依次由激光器(1)、光学扩束装置、第一消色差透镜（6）、测试空腔（7）、第二消色差透镜（8）、Wollaston棱镜（9）和CCD摄像仪（10）组成；所述的光学扩束装置依次由第一透镜（2）第一孔径光阑（4）、第二透镜（3）、λ/4波片（15）、第二孔径光阑（5）组成。

2.利用权利要求1所述断层成像激光剪切干涉三维温度测量装置的温度测量方法，其特征在于包括以下过程：

一、预处理工作：

1-1、所述的断层成像激光剪切干涉三维温度测量装置在光线经过测量容器后所会聚的焦点位置上安装有Wollaston棱镜（9）；调整光路中各元件的位置，确保在未加测试空腔（7）时，CCD摄像仪（10）可以拍摄到以有限宽方式显示的预加干涉条纹；

1-2、选取4个以上的角度对容器进行拍摄并获得清晰的干涉图，其中至少有2个角度所在的方向相互垂直；

1-3、对所获得清晰的干涉图做二值化处理，把图像转化成只含有黑白两种颜色的二值图像；用Matlab读各个角度的二值图像，得到每个像素的颜色值，此时各个角度干涉图的有效像素的数量处理为48×291，每幅二值图像的颜色值分别存入一个颜色值数组；

1-4、按从下到上从左到右的顺序处理数据原则，计算每幅二值图像的每个像素点的像素条纹级数值，其中设每一列最下面的像素点的像素条纹级数值为0，同一列中若当前像素点的像素颜色与前一像素点的像素颜色产生黑到白或白到黑的变化则当前像素点的像素条纹级数值等于前一像素点的像素条纹级数值+1；若当前像素点的像素颜色与前一像素点的像素颜色未产生变化则当前像素点的像素条纹级数值等于前一像素点的像素条纹级数值，每幅图像的像素条纹级数值分别存入一个条纹级数值数组；

1-5、对每幅二值图像所得的条纹级数值数组进行插值化处理，得到条纹级插值化数值数组；插值化处理的具体方法为：按从下到上的顺序处理数据原则，从每一列的第一个像素开始，寻找与其像素条纹级数值相差为1的像素；然后计算这两个像素之间间隔的像素数k，并计算得到△k=1/(k+1)；对于两个像素之间的间隔像素，用△k依次乘以间隔像素序号并加到其像素条纹级数值上，其中间隔像素序号从1开始计数；然后从这一列的下一个像素继续开始上述步骤，直到该列的第48个像素；所有的列都计算完毕后，得到条纹级插值化数值数组；

1-6、将重建区域划分为一个n×n的正方形网格，也称为像素，N＝n×n为像素总数；在每个像素内f(x,y)是常量，f_j代表第j个像素内的常数值；光束被定义为有一固定宽度的直线，把光束一一编号，编号用p的下标i表示；w_ji是第i个光束通过第j个像素的面积占整个第j个像素面积的比例，称做权重因子；

1-7、用每条光线经过的所有像素各自的权重因子w_ji和f_j的乘积表示每个像素对每条经过其本身的光线的贡献，一条光线经过的所有像素的这些值的和称作一个投影，这样，每个像素的函数值f_j和每条光线的投影值p_i的关系可以用下式表示：

$p_i=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{w}_{\mathrm{ij}}{f}_j,i=\mathrm{1,2,3},...M$

其中，M为投影总数；

二、迭代过程

对于上面获得的投影值和权重因子之间的方程组，利用迭代法求解；

先要把方程(1)改写成以下形式：

$\begin{array}{c}{w}_{11}{f}_1+w_{12}{f}_2+w_{13}{f}_3+...+w_{1N}{f}_N=p_1\\ w_{21}{f}_1+w_{22}{f}_2+w_{23}{f}_3+...+w_{21N}{f}_N=p_2\\ .\\ .\\ .\\ w_{M1}{f}_1+w_{M2}{f}_2+w_{N3}{f}_3+...+w_{\mathrm{MN}}{f}_N=p_M\end{array}---\left(2\right)$

N个像素代表了一个图像的N个自由度，所以，(f₁,f₂,…f_N)可以被认为表示N维空间的一个点，在方程组(2)中的每一个方程，都可以被看作是这个空间的一个超平面；当这个方程组存在唯一解，就意味着这M个超平面相交于一点，该点即为唯一解表示的点；如果方程组存在唯一解，迭代就会收敛于唯一解表示的点；因此，通过迭代方法可以求得f_j的值，f_j即第j个像素的折射率，求解每个像素的折射率；

三、温度场的反演计算

利用上述求解获得的折射率的分布，通过折射率与流体密度之间的关系式n-1＝Kρ获得流体密度场分布数据，其中K为Gladstone-Dale常数，它具有1/ρ的量纲，并且与气体的特性及所用的光的频率或者波长有关，(n-1)称之为折射率；进一步通过流体密度与温度之间的线性变化关系式ρ＝f(T)，反演推导出温度场的分布；其中ρ＝f(T)这一公式根据具体的流体而定。

Images