基于射影几何的流体光强三维重建方法

摘要

本发明公开了一种基于射影几何的流体光强三维重建方法，具体按照以下步骤进行：首先，采用多台相机围绕流体拍摄图像；其次，确定相机的内参数与外参数，最后，利用流体投影到图像的射影关系，建立流体亮度与图像亮度的方程，采用优化算法求解方程得到流体亮度的三维重建结果。本发明无需采用机械定位装置对相机进行摆放，所需的相机数量少，是一种快捷、准确流体三维光强重建方法，为进一步的流体温度和密度的三维重建提供了基础。

说明书

技术领域

本发明属于火焰、等离子体流体和计算机视觉技术领域，特别是涉及一种基于射影 几何的流体光强三维重建方法。

背景技术

随着航空、航天、制造业及能源工程等的快速发展，关于火焰、等离子体流场的三维诊断技术研究备受关注，三维诊断结果为飞行器相关理论的发展以及参数优化提供了数据基础。

基于光学的流场诊断方法可以实现对流场的非接触、无扰动、瞬态、定量检测，能够在不干扰流场的前提下实现动态瞬时测量和空间立体测量。

光学层析技术(OCT)是常见的基于光学的流场测量和诊断技术，其优点表现在：能够 实现全场三维定量测量，突破了空间单点或单平面检测的局限性；通过同时采集不同类型的光学投影信息，通过反演算法，能够重建多种被测流场的关键参数，例如温度、电 子密度的三维空间分布。

过去，光学层析技术得到了一定的发展，A new method for density fieldgradient measurements in compressible flows.J Stricker,O Kafri.Aiaa Jouranl，1982，20(6):820-823。 利用莫尔偏折法进行流场三维重建，其中使用逆Abel变换法，但其只能用于轴对称流场 的测量。

Computed Tomography of Chemiluminescence(CTC):In-stantaneous 3Dmeasurements and Phantom studies of a turbulent opposed jet flame.J Floyd，PGeipei.Combustion&Flame， 2011，158(2):376-391。将五个CCD相机间隔36°放置于圆形实验平台上，并利用反射 镜同时形成两个方向的投影成像于每个相机上，因而该实验装置理论上可获取10个方向 的流场投影信息，并基于简化的平行投影模型重建了湍流冲射火焰的三维结构。

Implementation of tridirectioanal large lateral shearing displacementinterferometry in temperature measurement of a diffused ethylene flame.HCZhou，JR Zhu，W Lv.Applied Optics，2011，50(21):3924-36。采用光学干涉法，建立了较大横向剪切距离的三方向剪 切光干涉CT装置，并实现了火焰的三维温度场重建。

基于图像的启发式三维火焰重建算法.沈亮，朱登明，周果，王兆其.高技术通讯，2017，27(2):122-130。通过迭代优化求解能量约束模型，实现对火焰三维温度场的重建。该方法基于火焰的物理特性，没有使用相机的投影模型，且对相机数目要求较高，需要 12至16台相机。

上述流场的三维诊断方法主要存在如下几个问题：

1.假设流场是稳定的、特殊的轴对称流场。而实际是剧变、任意空间形状的。

2.假设流场与相机之间满足平行投影关系，即流场空间各层之间无相关性，但是在实 际中，只有当流场离相机足够远时才可以近似地认为流场发出的是平行光，因此此类方 法不满足基本的相机成像原理。

3.光学干涉层析技术由于测量灵敏度较高，在测试环境比较恶劣的情况下，容易受到 振动等因素的影响，测量误差将会增加。当待测场变化剧烈时，还会产生干涉条纹断裂 等更为严重的情况，使测量难以进行。另外使用的干涉仪器设备较为昂贵。

4.数值模拟具有高效率和低成本的特点，主要是用数值方法将流场的控制方程组离散 到一系列网格节点上，然后求出其离散数值解，并且整个过程遵循质量守恒定律、动量 守恒定律和能量守恒定律。由于复杂流场中通常会包含非平衡流、两相流、非定常流效应等，并且可能会伴有剧烈而复杂的化学反应，所以很难找到一个能够完整描述被测复 杂流场分布、反映其内在本质的理论模型。综上所述，尽管数值模拟方法在近年来取得 了很大的进展，是流体力学研究的重要方法之一，但是在复杂流场的研究中仍然受到了 制约。

针对光学层析方法的不足，近年来，学者从相机成像原理方面进行了三维重建的研 究:

基于灰度值的非对称等离子体弧形貌的三维重建.郑庆余，刘新，孙晶，孔红领.机电工程.2012，29(7):750-753。将18个方向的CCD图像进行立体匹配与插值，从亮度重 建出温度相关的非对称等离子体弧三维形貌，但精度低，机械装置复杂。

基于立体计算机视觉的火焰三维重建算法的研究.任李懋.华北电力大学.2015。利用 火焰的图像特征，提出了一种基于特征点匹配的双目视觉重建方法，但是由于火焰各个 区域存在相似性，且各相机观察视角差异大，特征点的提取与匹配困难，容易出现大量误匹配。

Direct three dimensio-nal tomography of flames using maximization ofentropy technique. Anurag Goyal，Combustion and Flame.2014，161(1):173-183。利用最大熵方法进行火焰的 三维重建，其将一个像素的成像光路视为有一定粗细度的射线，计算射线与空间三维立 方体格子的相交部分的体积，将此体积与立方体格子体积的比值作为该格子对这个像素 的亮度贡献的权重。这种方法需要计算三维直线方程以及射线形成的锥体与立方体表面 的交点，空间相交关系复杂，计算困难。

Investigation into the effect of Ar ration in He-Ar mixture sideassisting gas on plasma properties during CO₂laser>

发明内容

本发明的目的在于提供一种基于射影几何的流体光强三维重建方法，解决了现有技 术中光学层析模型不符合相机的射影几何投影关系；另外，现有方法需要相机或者视角数量较多，相机的机械定位、安装要求高的问题。

本发明所采用的技术方案是，基于射影几何的流体光强三维重建方法，具体按照以 下步骤进行：

步骤1、将多台相机围绕流体进行摆放并垂直于喷射方向，在任意一个180°范围内分布，各相机采用同步方式触发同时拍摄，在各相机所拍摄的图像中，流体所在区域只 有流体和暗黑背景；

确定流体与相机的关系，流体满足光学薄特性，即相机拍摄图像上像素的亮度I为该像素和相机光心的连线所通穿过的流体各部分光强D(l)的积分，即

步骤2、确定相机的内参数与外参数，并建立标准世界坐标系；

步骤3、根据步骤1流体和相机的关系，以及步骤2确定的相机的内参数、外参数以及建立的标准世界坐标系，进行流体三维重建。

进一步的，所述步骤2具体按照以下步骤进行：

步骤2.1、使用相机对棋盘格进行不同角度、不同距离的多次拍摄，利用拍摄得到的 棋盘格图片，使用棋盘格标定板对各台相机进行相机内参数标定，使用相机标定原理进行计算，得到各台相机的内参数K和畸变系数；

步骤2.2、拍摄三个或三个以上半径已知的圆球，其中三个圆球处于空间三个不同位 置；

步骤2.3、对拍摄的圆球图像利用畸变系数进行去畸变处理，然后提取圆球图像边缘 二次曲线c，利用圆球图像边缘二次曲线c、相机内参数K和圆球半径s计算每个圆球球心在图像对应相机坐标系下的非齐次坐标；

步骤2.4、利用步骤2.3提取得到的三个圆球球心的相机坐标系下的非齐次坐标，以 三个圆球球心所在的平面为XOY平面，建立标准世界坐标系即XYZ坐标系，求解相机 在标准世界坐标系中的旋转矩阵和平移向量，即相机的外参数。

进一步的，所述步骤2.3具体按照以下步骤进行：

步骤2.3.1、对每台相机拍摄的圆球图像进行去畸变处理；

步骤2.3.2、提取圆球投影的边缘，拟合得到圆球图像边缘二次曲线c；

步骤2.3.3、根据圆球射影几何关系，X_w、Y_w、Z_w坐标系是世界坐标系，X_c、Y_c、Z_c坐标系是相机坐标系，平面Ⅱ₁为相机成像的图像平面，空间平面Ⅱ₂是Z_w轴穿过圆球球>w轴的平面，平面Ⅱ₂不经过圆球球心的图像平面，Q_c为正圆锥，以相>w，圆C_w以及它的圆心O_c位于平面Ⅱ₂上，O_c到>w的半径为r，圆球半径为s，位于平面Ⅱ₂上的点X的齐次坐>T，点X在平面Ⅱ₁上的投影点x为：

其中，K为相机内参数，R为世界坐标系相对于相机坐标系的旋转矩阵，P是在世界坐标系下相机的投影矩阵，t是世界坐标系相对于相机坐标系的平移向量；

圆球图像边缘二次曲线c与C_w之间满足平面Ⅱ₁到平面Ⅱ₂的单应关系：

c＝H^-TC_wH^-1＝K^-TR^-Tdiag(1,1,-(r/h)²)R^-1K^-1；

其中，h为O_c到相机光心的距离，H^-T为H的转置逆，H^-1为H的逆，K^-T为相机内>-T为R的转置逆，K^-1为K的逆，diag为构造对角矩阵的函数；

H为平面II₂和平面Ⅱ₁之间的单应矩阵；

通过奇异值分解可得到常数h/r的值和矩阵R，球心到相机光心的距离为d，则

在相机坐标系下，球心位置的非齐次坐标为：依次计算出三个圆球球心在相机坐标系下的非齐次坐标其中，i＝1,2,3。

进一步的，所述步骤2.4具体按照以下步骤进行：

步骤2.4.1、以三个圆球球心所在的平面为XOY平面，设其中一个拍摄圆球球心O1为标准世界坐标系原点O，球心O1与另一个球心O2连线记为标准世界坐标系的OX轴， 利用两个矢量O1O2与O1O3进行叉乘，得到标准世界坐标系的OZ轴，由OX轴与OZ 轴的两个坐标轴矢量进行叉乘，可得到三个圆球球心同平面的标准世界坐标系的OY轴， 从而利用三个圆球球心建立了以三个球心同平面的标准世界坐标系，进一步求解球心O1 在标准世界坐标系下对应的世界坐标为球心O2在标准世界坐标系下对应的世界坐 标为球心O3在标准世界坐标系下对应的世界坐标为

步骤2.4.2、利用三个圆球球心在相机坐标系下的非齐次坐标其中，i＝1,2,3，标准世界坐标系下对应的世界坐标其中，i＝1,2,3，其中， 建立两个坐标系之间的旋转平移关系：

其中R_w为标准世界坐标系相对于相机坐标系的旋转矩阵，t_w为标准世界坐标系相对于相机坐标系的平移向量。

进一步的，所述步骤3具体按照以下步骤进行：

步骤3.1、利用步骤2计算得到的相机内参数和外参数以及步骤1流体在各个相机拍 摄图像的位置，确定流体空间范围；

步骤3.2、将流体所在空间划分成空间立方体格，流体在一个空间立方格内的光强均 匀，根据流体的光学薄特性，相机拍摄图像上每一个像素的亮度是由空间立方体格沿投影方向投影到该像素上的光强加权叠加而成；

对流体所在的空间范围进行立方体格划分，立方体格总数目为m，m＝M×N×Q，M、N、Q分别为标准世界坐标系的X、Y、Z轴三个方向上的立方体格数目；

步骤3.3、空间立方体格的投影并不会刚好覆盖某个像素，需计算空间立方体格光强 对图像上每个像素亮度贡献的权重，按照相机内外参数所确定的射影几何关系，将空间立方体格的八个顶点投影到图像上，提取空间立方体格的顶点投影的凸多边形，计算凸 多边形占据图像各个像素的面积比例作为空间立方体格对像素亮度的权重系数；

步骤3.4、按照步骤3.3，对每个空间立方体格都进行操作，建立权重系数矩阵W；

步骤3.5、对各个相机所拍摄的图像，建立像素亮度与空间立方体格光强之间的约束 方程，

将各个视角所拍摄图像上每一个像素的亮度值记为x＝[x₁，x₂，…，x_n]^T，空间立方体>1，D₂，…，D_m]^T，满足

其中，n为像素数目，m为立方体格数目，约束方程即x＝WD，W是权重系数矩阵；

采用最小二乘优化准则迭代求解约束方程，在约束条件为x-WD≤0,D≥0，损失函数Loss(D)取最小值时，D为最优解，得到各个空间立方体格的光强，从而得到流体 的三维光强重建；

进一步的，所述步骤3.3具体按照以下步骤进行：

步骤3.3.1、每个空间立方体格具有八个顶点，按照步骤2.4.2中的标准世界坐标系与 相机坐标系的关系，空间立方体格的三维顶点(X_w，Y_w，Z_w，1)^T投影到图像上，满足

其中，Z_c为图像景深即图像在相机坐标系下的Z坐标，(X_w，Y_w，Z_w，1)^T为三维顶点>T为三维顶点投影到图像上的点>

步骤3.3.2、通过获取八个顶点在图像上的投影点，黑点为空间立方体格的顶点在图 像上的投影，提取空间立方体格的顶点投影的凸多边形；

步骤3.3.3、计算凸多边形占每个像素的面积比例，凸多边形全覆盖的像素，权重为 1，凸多边形未覆盖的像素，权重为0，部分覆盖的像素，以覆盖部分占像素面积的比值作为权重系数，权重系数大于0小于1。

进一步的，所述步骤2.3.2中提取圆球投影的边缘采用亚像素边缘检测算法。

进一步的，所述步骤3.3中提取空间立方体格的顶点投影的凸多边形采用三角剖分算 法和凸包提取算法。

进一步的，所述步骤3.5中所采用的最小二乘优化准则替换为最大均匀性准则、平滑 准则、最大熵准则、交叉熵准则、贝叶斯准则中的任意一种。

本发明的有益效果是，与现有技术相比，基于射影几何的流体光强三维重建方法，相机摆放无需机械定位装置，拍摄视角少，是一种快捷、准确流体三维光强重建方法， 并为流体的温度和密度的三维重建提供基础。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例 或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的 附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造 性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是光强分布和像素亮度的积分关系图；

图2是相机从多个方向向内拍摄示意图；

图3是圆球及其投影图；

图4是圆球组成的标准世界坐标系图；

图5是流体空间立方体格光强叠加图；

图6是流体空间立方体格投影到相机拍摄图像的凸多边形示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、 完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的 实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前 提下所获得的所有其他实施例，都属于本发明保护的范围。

基于射影几何的流体光强三维重建方法，具体按照以下步骤进行：

步骤1、将多台相机围绕流体进行摆放并垂直于喷射方向，在任意一个180°范围内分布，各相机采用同步方式触发同时拍摄，在各相机所拍摄的图像中，流体所在区域只 有流体和暗黑背景，如图2所示，采用180°范围而不是360°，如果采用360°透明流 体后面的相机会被摄到，影响光强积分；并且流体透明，180°即得到流体各部分的光强； 相机的台数根据精度的要求进行设置；

确定流体与相机的关系，流体满足光学薄特性，即相机拍摄图像上像素的亮度I为该像素和相机光心的连线所通穿过的流体各部分光强D(l)的积分，即如图 1和图5所示，火焰和等离子体满足光学薄特性；

步骤2、确定相机的内参数与外参数，并建立标准世界坐标系，具体按照以下步骤进 行：

步骤2.1、使用相机对棋盘格进行不同角度、不同距离的多次拍摄，利用拍摄得到的 棋盘格图片，使用高精度加工的棋盘格标定板对各台相机进行相机内参数标定，使用张正友发表的A flexible new technique for camera calibration中的相机标定原理进行计算，得 到各台相机的内参数K和畸变系数；

步骤2.2、拍摄三个或三个以上半径已知的圆球，其中三个圆球处于空间三个不同位 置；

步骤2.3、对拍摄的圆球图像利用畸变系数进行去畸变处理，然后提取圆球图像边缘 二次曲线c，利用圆球图像边缘二次曲线c、相机内参数K和圆球半径s计算每个圆球球心在图像对应相机坐标系下的非齐次坐标，具体按照以下步骤进行：

步骤2.3.1、对每台相机拍摄的圆球图像进行去畸变处理；

步骤2.3.2、通过亚像素边缘检测算法，提取圆球投影的边缘，拟合得到圆球图像边 缘二次曲线c；

步骤2.3.3、根据圆球射影几何关系，如图3所示，X_w、Y_w、Z_w坐标系是世界坐标系，X_c、Y_c、Z_c坐标系是相机坐标系，平面Ⅱ₁为相机成像的图像平面，空间平面Ⅱ₂是Z_w轴>w轴的平面，平面Ⅱ₂不经过圆球球心，Q_c为正圆锥，以相>w，圆C_w以及它的圆心O_c位于平面Ⅱ₂上，O_c到>w的半径为r，圆球半径为s，位于平面Ⅱ₂上的点X的齐次坐>T，点X在平面Ⅱ₁上的投影点x为：

其中，K为相机内参数，R为世界坐标系相对于相机坐标系的旋转矩阵，P是在世界坐标系下相机的投影矩阵，t是世界坐标系相对于相机坐标系的平移向量；

圆球图像边缘二次曲线c即图3中平面Ⅱ₁上的曲线c与C_w之间满足平面Ⅱ₁到平面Ⅱ₂的单应关系：

c＝H^-TC_wH^-1＝K^-TR^-Tdiag(1,1,-(r/h)²)R^-1K^-1；

其中，h为O_c到相机光心的距离，H^-T为H的转置逆，H^-1为H的逆，K^-T为相机内>-T为R的转置逆，K^-1为K的逆，diag为构造对角矩阵的函数；

H为平面II₂和平面Ⅱ₁之间的单应矩阵；

通过SVD(奇异值)分解可得到常数h/r的值和矩阵R，圆球的半径s已知，球心 到相机光心的距离为d，则

球球心在相机坐标系下的非齐次坐标其中，i＝1,2,3；

步骤2.4、利用步骤2.3提取得到的三个圆球球心的相机坐标系下的非齐次坐标，以 三个圆球球心所在的平面为XOY平面，建立标准世界坐标系即XYZ坐标系，求解相机 在标准世界坐标系中的旋转矩阵和平移向量，即相机的外参数，具体按照以下步骤进行；

步骤2.4.1、以三个圆球球心所在的平面为XOY平面，设其中一个拍摄圆球球心O1为标准世界坐标系原点O，球心O1与另一个球心O2连线记为标准世界坐标系的OX轴， 利用两个矢量O1O2与O1O3进行叉乘，得到标准世界坐标系的OZ轴，由OX轴与OZ 轴的两个坐标轴矢量进行叉乘，可得到三个圆球球心同平面的标准世界坐标系的OY轴， 从而利用三个圆球球心建立了以三个球心同平面的标准世界坐标系，如图4所示，进一 步求解球心O1在标准世界坐标系下对应的世界坐标为球心O2在标准世界坐标系 下对应的世界坐标为球心O3在标准世界坐标系下对应的世界坐标为

步骤2.4.2、利用三个圆球球心在相机坐标系下的非齐次坐标其中，i＝1,2,3，标准世界坐标系下对应的世界坐标其中，i＝1,2,3，其中， 建立两个坐标系之间的旋转平移关系：

其中R_w为标准世界坐标系相对于相机坐标系的旋转矩阵，t_w为标准世界坐标系相对于相机坐标系的平移向量；当圆球为3个时，R_w为3×3的旋转矩阵，t_w为3×1的平移>

步骤3、根据步骤1流体和相机的关系，以及步骤2确定的相机的内参数、外参数以及建立的标准世界坐标系，进行流体三维重建，具体按照以下步骤进行；

步骤3.1、利用步骤2计算得到的相机内参数和外参数以及步骤1流体在各个相机拍 摄图像的位置，确定流体空间范围；

步骤3.2、将流体所在空间划分成空间立方体格，流体在一个空间立方格内的光强均 匀，根据流体的光学薄特性，相机拍摄图像上每一个像素的亮度是由空间立方体格沿投影方向投影到该像素上的光强加权叠加而成，如图5所示；

对流体所在的空间范围进行立方体格划分，立方体格总数目为m，m＝M×N×Q，M、N、Q分别为标准世界坐标系的X、Y、Z轴三个方向上的立方体格数目；

步骤3.3、空间立方体格的投影并不会刚好覆盖某个像素，需计算空间立方体格光强 对图像上每个像素亮度贡献的权重，按照相机内外参数所确定的射影几何关系，将空间立方体格的八个顶点投影到图像上，采用三角剖分和凸包提取算法提取空间立方体格的顶点投影的凸多边形，如图6所示，计算凸多边形占据图像各个像素的面积比例作为空 间立方体格对像素亮度的权重系数，具体按照以下步骤进行；

步骤3.3.1、每个空间立方体格具有八个顶点，按照步骤2.4.2中的标准世界坐标系与 相机坐标系的关系，空间立方体格的三维顶点(X_w，Y_w，Z_w，1)^T投影到图像上，满足

其中，Z_c为图像景深即图像在相机坐标系下的Z坐标，(X_w，Y_w，Z_w，1)^T为三维顶点>T为三维顶点投影到图像上的点>

步骤3.3.2、通过获取八个顶点在图像上的投影点，图6中每一个格子表示一个像素， 黑点为空间立方体格的顶点在图像上的投影，利用Delaunay三角剖分算法和ConvexHull 凸包算法提取空间立方体格的顶点投影的凸多边形；

步骤3.3.3、计算凸多边形占每个像素的面积比例，凸多边形全覆盖的像素，权重为 1，凸多边形未覆盖的像素，权重为0，部分覆盖的像素，以覆盖部分占像素面积的比值作为权重系数，权重系数大于0小于1；

步骤3.4、按照步骤3.3，对每个空间立方体格都进行操作，建立权重系数矩阵W；

步骤3.5、对各个相机所拍摄的图像，建立像素亮度与空间立方体格光强之间的约束 方程，

图像上第i个像素的亮度为x_i，x_i与流体空间中第j个立方体格的亮度D_j的关系为x_i＝D_back+∑w_jD_j；其中，w_j为权重系数，D_back为背景光强，当D_back远小于D_j时，D_back为零；

将各个视角所拍摄图像上每一个像素的亮度值记为x＝[x₁，x₂，…，x_n]^T，空间立方体>1，D₂，…，D_m]^T，满足

其中，n为像素数目，m为立方体格数目，约束方程即x＝WD，W是权重系数矩阵； 此时各个视角拍摄时，几个相机是同一环境下统一触发拍摄的，背景光强均远小于每个 立方体格的亮度，所以背景光强为零；

采用最小二乘优化准则迭代求解约束方程，在约束条件为x-WD≤0,D≥0，损失函数Loss(D)取最小值时，D为最优解，得到各个空间立方体格的光强，从而得到流体 的三维光强重建；

步骤3.5中所采用的最小二乘优化准则替换为最大均匀性准则、平滑准则、最大熵准 则、交叉熵准则、贝叶斯准则中的任意一种。

本发明将在燃烧分析诊断方面有重大应用价值，火焰的三维重建可以帮助人们更加 直观地了解火焰的温度分布和燃烧的充分程度,对于研究火焰的传播机理和发展火焰模 型的动态仿真研究有重要意义，用于发动机尾焰温度分析、锅炉燃烧温度和充分度分析等领域。

在等离子体诊断领域，等离子体内部结构的三维光强分布与电子密度的关系密切相 关，可以配合静电探针诊断、激光干涉诊断、发射光谱诊断等方法，对等离子各个区域的电子密度进行分析，也能够用在激光焊接等离子体诊断领域中。

本说明书中的各个实施例均采用相关的方式描述，各个实施例之间相同相似的部分 互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于系统实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法 实施例的部分说明即可。

以上所述仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。 凡在本发明的精神和原则之内所作的任何修改、等同替换、改进等，均包含在本 发明的保护范围内。