基金项目:国家自然科学基金(11927802); 中国博士后科学基金(2023M730168)
作者简介:胡 炜(1994—),女,博士,副研究员,研究领域为背景纹影技术及其在燃烧测量的应用。
通信作者:李敬轩(1984—),男,教授,研究领域为不稳定燃烧的预测、测量和控制。
1.北京航空航天大学 宁波创新研究院,浙江 宁波 315800; 2.北京航空航天大学 宇航学院,北京 100191
1.Ningbo Institute of Technology, Beihang University, Ningbo 315800, China; 2.School of Astronautics, Beihang University, Beijing 100191, China
background-oriented schlieren; tomographic technique; three-dimensional reconstruction; incomplete angle; ill-posed inverse problem
DOI: 10.3969/j.issn.1672-9374.2024.06.001
实现对非定常流动的三维测量一直是航空航天、能源环境、燃烧诊断等领域的研究重点[1-2]。光学诊断技术因其非接触、无干扰的优势成为近年来的研究热点[3]。在众多光学诊断技术中[4-5],背景纹影技术(background-oriented schlieren, BOS)是近20年来兴起的一种新型流场测量技术[6-7]。该技术利用光在非均匀流场中会发生偏折的原理,通过比较有无流场前后的背景图像偏差和对沿光程积分信息的反解,仅需相机和背景就能动态测量非定常流场,极大地简化了测量装置的结构与光学布局。因此背景纹影技术具有大视场测量、对环境要求低、不干扰流场、成本低等优点,能更好地满足在线测量的工程需要。
经过20余年的发展,BOS已在人文[8]、健康[9]、能源[10]、建筑[11]、电力[12]、航空航天[13]等领域崭露头角,其应用场景涉及超声速流动、燃烧环境、等离子体流动等[6]。通过物理方程转换,BOS可测量内容包括折射率场[14]、密度场[15]、温度场[16]与速度场[17]等。结合层析技术与多相机技术,BOS可重构流场维度覆盖二维[18]与三维[19],可重构流场结构包含层流[20]与湍流[21]。随着背景板的不断扩大与自然化,小至毫米级的湍流火焰[22],大至数米级的直升机桨叶旋涡[23],BOS均证明了其全尺寸测量能力。配合高速相机,BOS还可实现流场的高频动态测量。美国NASA Bathel等[24]采用4个高速相机实现了10 kHz的氦气射流流场测量。
需要注意的是,在大视野、高频动态、低成本、好测量的优势背后,BOS仍面临高空间分辨率、高精度、高敏感度、高实时性等测量性能的挑战。例如,背景聚焦而待测流场失焦的测量方式使得拍摄到的背景斑点存在模糊圈[25],而图像模糊造成图像分辨率降低,进而影响系统的空间分辨率。其次,BOS中沿光程积分的反解用于重构流场信息,而边界条件对反解精度具有重要影响。尽管现有BOS采用狄利克雷边界条件或纽曼边界条件解决了一些二维边界情况[18],但其在三维重构时基本使用固定边界条件,即假设流场外部均匀到特定值[19]。这种简化假设与实际情况存在矛盾。例如对发动机腔体内燃烧场测量时或待测流场存在明显扩散的情况下,其边界条件均非零。再者,当BOS测量具有极端梯度(非常小或非常大的梯度)的流场时,其测量精度和灵敏度还受背景图案[26]、图像偏移算法[26]、近光轴假设[27]的限制。此外,计算成本一直是制约BOS在线测量的关键因素,特别是流场重构时对大规模线性方程的求解存在耗时长、占用内存大的问题。即使在GPU加速下,实时BOS目前也处于二维阶段,即实现单相机的图像偏移场测量[28]。同时,计算成本还会约束流场重构中所能处理的未知量数量和方程个数,从而限制空间分辨率和精度的上限。
为提升BOS性能,国内外多个知名团队以及美国NASA纷纷展开了BOS理论研究与应用。德国宇航中心Raffel[29]、美国宾夕法尼亚州立大学Settles和Hargather[30]、北京航空航天大学熊渊[6]发表了BOS代表性综述。德国Lavision、丹麦Dantec公司更是提出了商业化解决方案。按研究内容分,BOS方法的发展主要集中在硬件搭建、偏移场计算和流场重构3个方面。以瑞士苏黎世联邦理工学院Meier等[31]和澳大利亚蒙纳士大学Amjad等[32]为代表的学者们结合激光、全反射玻璃片、全息扩散片等光学器件,发展了二维/三维激光散斑BOS测量系统,即以激光散斑代替背景板,相机可以聚焦在任意位置,以解决自然光穿透性不强、时空分辨率不足等问题。千叶大学[33-34]、国防科技大学[35-36]结合远心光学结构设计了远心BOS测量系统,实现平行光测量、背景和待测流场双聚焦并获得更大的景深。上海交通大学[37]发展了光纤+单相机的BOS测量系统,实现单相机的三维流场动态测量。中国民航大学[38]设计了激光干涉BOS,将原先的斑点式背景替换成高灵敏度的条纹式背景。值得注意的是,随着BOS硬件部分的提升,系统的空间分辨率、测量精度和敏感度均有所改善,但常以损失视野范围或提高成本为代价。在保持原有BOS硬件布局的前提下,德国汉诺威莱布尼茨大学Goldhahn等[15]讨论了BOS硬件布局对敏感度、准确度和分辨率的影响。德国宇航中心Schwarz等[25]更是提出了能最大化测量敏感度的光路布局方法。
除了从光源选择、镜头选择、光路设计方面优化光学布局外,BOS软件部分也得到了发展,以进一步提高测量性能。在图像偏移场计算方面,各国学者研究了背景图案分布(如泊松圆盘、随机、小波、条纹)[26,39]、图案大小与密度[40-41]、图案颜色(如黑白、彩色)[42-43]对偏移场估计的影响,以期找到对流场梯度变化最敏感、偏移显示最明显的背景图案。为进一步提高有无流场前后的背景差,德国宇航中心Gardner等[44]和德国通快公司Reichenzer等[45]还提出了可移动式和可变式的新背景生成方案。在此基础上,学者们比较了互相关[46]、光流[47]、点追踪[48]、调制解调[49]、深度神经网络[50-52]方法下的图像偏移计算准确度,并探讨了各方法的适用条件。由费马定理可知,光线偏折角等于折射率梯度沿光程的积分。在获取图像偏移量并完成光线偏折角转换和矫正后,还需对积分方程反解,以重构流场信息[6]。因此重构方法是保证流场预估准确性的最后关键步。除了特定型流场(如轴对称型、轴向均匀型)可用二维泊松方程[53]、Abel变换[54]及其扩展重构外,对于一般型非均匀无规则流场,BOS需要结合层析技术利用多个视角的图像偏移信息进行重构。学者们分别从重构基函数、重构顺序、重构算法、重构优化等方面展开了研究。
需要注意的是,由于待测流场未知量个数的庞大性,BOS需要多相机密集布置,获得大样本积分方程,以避免求解陷入不适定。德国柏林工业大学Lang等[55]发现投影数量需超过10个才能保证重构结果的准确性。然而对于实际测量条件而言,受空间、成本限制,可用的测量光路并不多,因此BOS只能在有限角度内稀疏采集投影数据。投影数据的不完备与不可避免的投影噪声使得BOS处于高度不适定的重构过程中。显然,这种有限角度内稀疏测量的BOS系统(命名为不完全角度BOS)容易造成重构结果不准确与过拟合。如何提高不完全视角下层析BOS的流场重建质量,是燃烧诊断与先进流体测量中的研究要点和亟待解决的应用难点。
本文旨在对现阶段不完全角度层析重构技术在BOS测量中的发展与应用情况进行综述,从以下3个方面展开:首先回顾BOS成像原理、层析重构原理以及BOS层析测量系统,讨论BOS测量时面临的不完全角度问题; 其次对现阶段不完全角度BOS层析技术重构算法进行综述; 最后对BOS层析测量技术的发展提出展望。
图1是基于薄透镜相机模型的背景纹影成像系统,由背景图案、待测流场和成像系统3个部分组成。
当待测流场具有非均匀折射率时,产生的折射率梯度使得来自背景图案的光线以一定偏转角发生偏折,引起光线透过镜头入射到投影面上的位置产生偏移。由费马定理可得,光线在介质中传播的光线方程为
式中:p=[x,y,z]T为光传播路径ray上任意一点的位置; n(p)为在p点处的折射率; ds=为光传播路径的积分变量; SymbolQC@为空间梯度算子。
式(1)可分解为一阶常微分方程组,即
式中J为被p点处折射率缩放后的光线方向向量。
当除待测区域外的环境折射率n0恒定时,偏折角ε=[εx,εy,εz]T等于折射率梯度沿光传播路径的积分[21],即
式中:Jin和Jout分别为光线进入和离开待测区域时被缩放后的光线方向向量; pin和pout分别为光线进入和离开待测区域时的位置。若传播过程中介质折射率n(p)≈1,则偏折角还可简化为ε=。但若传播过程中除待测区域外的环境折射率既不恒定、又无法确定、还不可近似为1时,偏折角公式应为ε=。现有BOS测量都基于式(3)展开。式(3)又可表达为
式中εα为α方向上的实际偏折角。
假设光线实际传播路径接近光轴,即偏折角满足小角度时,实际发生偏折的光线离开待测区域的位置pout与初始不发生偏折的光线离开待测区域的位置poutO几乎一致。此时,偏折角可由背景斑点的偏移场推测得到,即
式中:ε'α为α方向上poutOpbk和poutOpbkO形成的近似偏折角,pbkO为背景斑点发生偏移前位置,pbk=pbkO+[Δx Δy Δz]T为背景斑点发生偏移后的位置; ZBO为背景板到镜头的距离; ZPO为流场中心到镜头的距离; Δα为背景斑点在α方向上的偏移量,其与投影面的图像偏移u和v成比例。
假设BOS测量时有K个投影,每个投影由I×J个像素组成。由相机标定模型[56]可知,偏折角与图像偏移间的关系为
式中:K(k)∈瘙綆3×3为第k个相机的内参; R(k)∈瘙綆3×3为统一相机坐标系下第k个相机的绝对旋转外参; Z(k)c为第k个相机的尺度系数; Z(k)BP=Z(k)BO-Z(k)PO为第k个背景到待测流场中心的距离。
以上这些参数均可通过相机标定获得[21]。该过程先将像素坐标系中第k个投影的图像偏移[记为u(i,j)(k)和v(i,j)(k)]转为第k个相机坐标系中的背景偏移,再将其转换到统一相机坐标系下,并利用式(5)转化为第k个背景发出光线对应的偏折角εα(i,j)(k),α∈{x,y,z}。而投影面上有无流场前后的图像偏移可借助知名开源软件包PIVlab[57]、OpenCV[58]等,通过互相关、光流等方法计算得到。
背景纹影测量的核心是通过图像处理与几何光学转换得到偏移角的观测值,而后逆求解式(4)以重构折射率场。为方便求解,一般在近光轴假设下,将光线偏折角εα(i,j)(k)近似看作是由沿初始未发生偏转的光程上折射率梯度的积分。此时,式(4)可近似线性离散成为
式中ε(k)α, α∈{x,y,z}包含第k个投影面上各光线沿α方向的偏折角。
将待测流场划分成L×W×H个体素,P(k)为第k个层析投影矩阵,由每条直线光线沿着每个体素的弦长s组成,如图2所示。
图2 BOS层析离散示意图[19]
Fig.2 Tomographic BOS discretization illustration[19]
▽n和n分别为待测流场各体素上的折射率梯度和折射率,本文统称为未知量矩阵x。而流场特定点处的未知量x可由各体素基函数及其对特定点处的权重和表示,即
式中基函数gi包括阶跃基函数、线性基函数、高斯基函数、Kaiser-Bessel基函数等。美国宾夕法尼亚州立大学[54]、南京理工大学[59]、加拿大英属哥伦比亚大学[60]分别研究了阶跃型、线型、高斯型、径向对称Kaiser-Bessel基函数对重构结果的影响。
对折射率场的重构顺序目前有4种,如图3所示。第一种(间接法1)是先重构三维折射率梯度分布SymbolQC@n,再求解三维泊松方程来获得折射率[60]。当折射率梯度在边界为零时,第二种重构顺序(间接法2)先重构各投影面折射率积分分布n^-,再求解三维折射率[61-62],其推导参见文献[20]。第三种为目前最为常见的直接法,由法国航空航天实验室Nicolas等[19]在2016年首次提出,即将折射率梯度的微分算子通过有限差分离散为矩阵Dα,直接求解折射率,从而避免折射率梯度/积分病态求解误差在下一个病态方程求解中的放大。当折射率梯度变化不大时,美国宾夕法尼亚州立大学Grauer等[14]结合光流法提出一种更为精炼的统一直接法,通过推导偏移算子U(k)op和V(k)op,实现直接从图像强度变化来重构三维折射率场,跳过图像偏移估计以避免该估计误差在下一步病态方程求解中的放大。需要注意的是,相较另两种方法,间接法2和统一直接法的使用需要更严苛的流场条件[20,63]。
BOS层析重构模型如图3所示。由图3可知,无论何种重构顺序,近光轴假设下的式(4)终会简化成Ax=B的线性方程形式。由于待测流场未知量个数的庞大性,以及多相机密集布置引起的方程个数庞大性,使BOS重构成为一个经典的大规模线性方程求解问题。巨大的方程规模、不可避免的观测值噪声B以及理想假设下系数矩阵A与实际的差异性对重构算法的高效性、正确性和鲁棒性提出了挑战。加之实际不完全角度的测量限制,使重构结果存在诸多问题。本文将在第2节详细论述相关不完全角度下的BOS层析重构算法。
需要注意的是,BOS层析技术重构的是一个相对的值,如何选择边界条件对重构量值的正确性至关重要。如图4(a)所示,在二维折射率场的重构过程中[18],一般在区域边界上给定折射率(狄利克雷边界条件)或折射率梯度的值(纽曼边界条件)。而在三维折射率场的重构过程中[19,64],一般在某些面上施加固定边界,即假设边界处及边界以外区域折射率一致,见图4(b)白面。同时对其他无法确定边界条件的面施加自由边界,见图4(b)灰面。美国奥本大学[65]采用了三维纽曼边界条件,假设流场空间最外两层的折射率梯度相等。美国北卡罗来纳州立大学Wahls等[66]在测量一个用石英玻璃圆筒包围着的预混轴对称火焰的温度分布时,用红外热成像测量了玻璃表面温度,并作为边界条件提供给BOS测量内部温度。测量结果表明,提供更准确的边界条件能使测量准确性提升1.5%~3.8%。除了边界条件,筛选有效的光线也是组建Ax=B线性方程组时所要着重考虑的[19,64]。若光线在待测区域外仍会发生偏折,则被判定为无效光线,应当在构建模型时去除,如图4(b)中的红色光线所示。笔者曾尝试将该类光线引入方程组,结果加重了求解结果的不适定性。此外,待测流场的最小体素由BOS测量中解析的最小尺度决定,过大的体素会过平滑流场并影响数值结果[55]。
图3 BOS层析重构模型[46,48-49]
Fig.3 Tomographic BOS reconstruction model[46,48-49]
图4 BOS二维和三维重构过程中的常见边界条件[18-19]
Fig.4 Common boundary conditions in BOS two-dimensional and three-dimensional reconstruction processes[18-19]
在获得三维折射率场后,可通过Gladstone-Dale方程转换为第j个体素处的密度场ρj,即
nj=1+ρjGj (9)
式中:nj为第j个体素处的折射率; Gj为第j个体素处的Gladstone-Dale系数。结合理想气体状态方程,可推得第j个体素处的温度场Tj为
式中:R为通用气体常数; pj为第j个体素处的压强; Gi、Yi、Mi分别为混合气体中第i种组分的Gladstone-Dale系数、摩尔分数、摩尔质量。当待测流场中气体组分浓度分布与压强分布未知时,一般对其做均值处理,即假设待测流场各体素间的压强和组分浓度分布相同[1,21]。此外,结合其他物理方程,如流体连续性方程[67]、Navier-Stokes方程[68]、Euler方程[69],流体密度还可转换为气体不同状态下的速度场。
对于非轴对称或非均匀的不规则流场,典型的BOS层析测量系统包括以下两种:①旋转流场实现相对于单个相机的不同位置投影; ②在流场周围放置多个相机(通常为10个或更多),以获取大范围的投影[20,42,70-71],如图5所示。第一种仅适用于流场稳定的测量场景。例如笔者通过旋转流场,采用180°范围内10个角度下的投影信息,重构了层流预混双火焰的三维温度场[20]; 日本千叶大学Ota团队[70]和法国-德国圣路易斯研究所Sourgen团队[42]设计的旋转喷嘴,通过90°范围内19个角度下的投影信息分别重构了超声速非对称锥形钝体和突刺钝体的三维密度场。该测量系统的优势在于,可以以较低的空间成本与设备成本获得非常多的角度投影,从而避免不完全角度层析重构问题。对于具有特定变化规律的流场,如具有旋转对称性的旋涡射流,德国柏林工业大学Lang等[55]利用单相机锁相技术,通过360°范围内42个角度的投影信息,重建了时间和相位平均的非等温旋流三维温度场,如图6所示。
图5 BOS层析典型测量系统[20-21]
Fig.5 Classical system of tomographic BOS[20-21]
图6 非等温旋流三维温度场[55]
Fig.6 Temperature field reconstruction in a non-isothermal swirling jet undergoing vortex breakdown[55]
第二种多相机系统需要充足的相机和足够大的空间环境才能完成测量,但其所能重构的流场类型可以更为复杂多变。例如,法国航空航天实验室Nicolas等[19]、南京理工大学Cai等[59]、澳大利亚蒙纳士大学Amjad等[32]、加拿大英属哥伦比亚大学Atcheson等[72]、美国宾夕法尼亚州立大学Grauer等[21]分别采用12~23个普通相机测量了扩散火焰、热气流、燃烧器螺旋羽流等不稳定非规则流体的瞬时三维流场。然而由于使用的是普通相机,测量的时空分辨率都是受限的,远无法满足非定常复杂流动(如湍流火焰)的高频动态测量需求。在不考虑成本的前提下,理想方案是同时使用多台高速相机。目前报道的BOS中使用高速相机数量最多为4台,其测量频率达到了10 kHz[24]。需要注意的是,随着测量频率的上升,相机曝光时间减少,会对光源的亮度和脉冲宽度提出更高的要求。此外,随着相机数量的上升,对相机标定的准确性要求也会变高。
在高速相机数量受限的情况下,上海交通大学Liu等[37]、日本大阪工业大学Ukai[22]通过多个光线内窥镜和单台高速相机配合的方式,实现了千赫兹时间分辨率的非层流火焰三维折射率场和密度场测量,如图7所示。尽管光纤内窥镜测量方式可以减少测量成本并提高测量时间分辨率,但以损失视野范围和空间分辨率为代价。西北工业大学[73]指出当采样频率高于10 kHz时,光纤束会造成极大的光信号损失,进而大大降低图像信噪比。
图7 利用光纤内窥镜的单相机BOS层析测量系统[22]
Fig.7 Single-camera system of tomographic BOS using optical fiber endoscopes[22]
另外,多重瞬时层析测量结合多种层析测量技术实现多流场类型的同时测量正成为新趋势,如图8所示。德国杜伊斯堡-埃森大学[74-76]首次结合BOS层析、光发射层析和盐发射层析这3种成像技术,采用43个相机实现了折射率(反映温度、压力和组分浓度)、CH*化学发光(反映火焰锋面)、SrOH发射(反映热喷雾火焰产物)以及热喷雾火焰产物温度4种流场三维分布的瞬时测量。多重瞬时层析测量提供了一种低成本的燃烧诊断替代方案,可评估火焰形态、区域限制、局部化学计量和温度等参数,对理解反应流和喷雾火焰合成机制具有重要意义。
图8 多重瞬时层析测量系统[74]
Fig.8 Experimental setup of TIMes[74]
BOS的层析重构早期算法代表为滤波反投影算法(filtered back projection, FBP)。该算法基于平行光束假定和中心切片定理:二维非均匀场f(x,y)投影P(r,θ)的一维傅里叶变换F[P(r,θ)],等于f(x,y)二维傅里叶变换F(u,v)在平面上沿着θ方向过原点直线上的值,如图9所示。因此,为了覆盖二维非均匀场的整个傅里叶空间以实现高精度重构,需要获得180°范围内的全部投影数据。受空间、成本的限制,能够获得完备投影数据的测量条件并不常见,一般只出现在用单相机测量稳定流场的情况[15,42]。当投影数量无法满足时,FBP算法的重构精度将受到严重影响。
图9 FBP算法原理[77]
Fig.9 FBP principles[77]
需要注意的是,由图2可知,BOS的光束模型实际上为锥形束。因此,在二维非均匀场的BOS测量中,扇形束假设比平行光束假设更贴近真实。此时中心切片定理不成立,需要获得360°范围内的全部投影数据才能实现高精度重构,进一步加大了对投影数量的需求。图 10(b)和图 10(c)比较了平行束逆Radon变换和扇形束逆Radon变换下,采用180°范围内的全部投影数据进行重建的非均匀场质量。可以看到,180°范围内的全部投影数据无法支撑扇形束逆Radon变换完成准确的流场重构。
图 10 平行束逆Radon变换和扇形束逆Radon变换下流场重构结果
Fig.10 Reconstructed results of the flow field under the inverse Radon transform of parallel beam and fan beam
本节将在扇形光束假设下讨论BOS所面临的不完全角度问题。不完全角度问题隶属于不完全数据问题的一种。根据数据缺失情形,不完全数据问题可分为3种[78],如图 11所示。
1)内问题。光线视场不能覆盖待测流场的全部空间,仅采集到流场内部区域全部角度下的投影数据。由于BOS重构依赖准确的边界条件,而内问题下流场内部区域的边界条件不好获得,外部区域又缺乏全部投影数据,造成内部和整体流场的求解都有困难。
2)外问题。待测流场内部存在光线(自然光、LED光、激光)无法穿透的物体,如金属火焰、碳烟等,该物体所在区域在每个角度下都无法得到投影数据。
3)不完全角度问题。某些角度下无法获得投影数据,又可分为稀疏角度问题和有限角度问题。稀疏角度问题指相邻扫描角度间存在间隔,通常情况下角度是等间距且覆盖360°区间。由于投影数据等角度间隔缺失,重建的流场会出现图 10(d)所示的条纹状伪影,且图像细节存在模糊。有限角度问题指连续扫描角度范围小于360°。由于在一段连续角度范围内损失投影信息,重建的流场会出现图 10(c)所示的严重畸变和信息丢失。在BOS实际测量时,稀疏角问题和有限角问题往往同时存在。
从数学角度看,BOS层析重构属于大规模线性方程求解问题。投影数据的缺失使得大规模方程变得病态和欠定,而不可避免的观测值噪声和理想假设下系数矩阵与实际的差异性,加重了对病态方程求解的不适定性。相较内问题和外问题,不完全角度问题在实测中更为普遍。
图 11 BOS测量中3类不完全数据问题[78]
Fig.11 Three types of incomplete data issues of BOS measurement[78]
对于BOS病态矩阵Ax=B求逆问题,最常用的处理方法是通过添加适当的约束条件,使不适定问题转变为适定问题。正则化方法通过给目标函数增加平滑惩罚项,使新目标函数中相对稳定且有意义的解从大量振荡无意义的解集中脱颖而出。新目标函数的表达式为
x*=argminx[‖Ax-B‖2+λL(x)] (11)
式中:‖Ax-B‖2为数据拟合项; λ≥0为正则化系数,控制着惩罚项和数据拟合之间的权衡。需要注意的是,λ的取值对重构结果有很大影响,过大的λ容易过平滑流场,而过小的λ不足以形成适定解。目前常用的λ优化方法为L曲线法和广义交叉验证法(generalized cross-validation, GCV)。L曲线法通过绘制对数-对数尺度下不同正则化系数的惩罚项log‖Ax*(λ)-B‖2和数据拟合项间log L[x*(λ)]的曲线,取曲率最大点作为最优正则化系数点[19-20,71,79]。GCV法[80-81]的基本原理为即使数据拟合项中的任意一条数据被移除时,仍能从剩余数据得到的正则化解中估计出移除的数据项。通过代数变换构造一个广义交叉GCV函数,即
式中:C(λ)=A(ATA+nλI)-1AT; n为数据拟合项个数; trace(·)表示矩阵的迹,即矩阵对角线上的元素和; I为单位矩阵。通过最小化G(λ),求得最优正则化系数点,从而确保使用尽可能少的奇异值来拟合数据项。其他优化方法见文献[80]。
L(x)为对流场三维分布先验编码的正则化项,又称惩罚项。正则化项的设计是为了鼓励流场的平滑性和连续性。当未知量数量小于方程个数时,正则化项充当“过滤器”的角色,通过最小化流场空间奇异值来抑制误差放大。当未知量数量大于方程个数时,正则化项通过提供额外的先验信息来“填补未知量空白”,进而获得唯一解。常见的正则化项有吉洪诺夫(Tikhonov)正则化项,即
L(x)=‖Lx‖2 (13)
式中L为m阶梯度的离散算子▽(m)·x[14],例如m=2时,L矩阵可表示为
式中ηi表示与i临近的体素数。对“临近”的常见定义为相邻体素面接触[14,54,82],因此对三维流场而言ηi≤6,对二维流场而言ηi≤4,依次类推。当然“临近”也可以为相邻体素点接触[83],此时ηi≤26。另外需要注意当i处在流场角点、边沿、表面时,临近的体素数会有所减少。
另一种常见的正则化项为全变分(total variation, TV)。其中,1阶TV正则化的表达式为[84]
吉洪诺夫正则化与TV正则化的侧重点有所不同。前者强调流场的平滑度并减少噪声影响,但可能导致流场的过度平滑进而丢失流动中潜在的不连续性。后者着重于保持流场锐利的边缘和具有恒定梯度的区域。这两种方法目前在BOS的应用都十分广泛。法国航空航天实验室Lanzillotta等[85]通过1阶吉洪诺夫正则化技术对135°范围内等间距布置的8台相机的投影数据进行分析,重构了超音速过膨胀射流密度场的平均和瞬时三维结构。南京理工大学Cai等[59]采用1阶TV正则化重构了160°范围内等间距布置的15个相机所拍摄的扩散火焰的投影数据。法国航空航天实验室Nicolas等[19]对比了两种正则化方法后,发现对于相对密度梯度较小的对流流动,不需要TV正则化来得到适定解。该结论在文献[1,86]的仿真试验中也得到了验证。相较吉洪诺夫正则化,TV正则化更适合捕捉流场的小尺度特征,但对流场整体结构的把握度差。因而在旋流预混火焰试验中(小密度梯度的对流流动),吉洪诺夫正则化能更好地反映流场结构。美国宾夕法尼亚州立大学Grauer等[21]通过贝叶斯求解框架将两种正则化整合,取吉洪诺夫正则化求得的流场作为初始值,代入TV正则化中对流场做进一步细化。法国航空航天实验室Todoroff等[87]将正则化目标函数的最小化问题通过变量分裂法分解为两个子正则化问题,以解决吉洪诺夫正则化容易过平滑的问题。通过比对6个、13个、23个投影数据的重构结果,发现正则化后流场重构误差可以减少50%。此外,吉洪诺夫正则化与TV正则化技术也广泛应用在其他光学层析重构技术中[82,84]。如青岛科技大学Zhang等[88]提出一种整合吉洪诺夫和TV的混合正则化技术。上海交通大学[83]在研究不完全角度化学发光层析技术时,发现当投影角度较少时(倾向于秩缺陷逆问题),由于需要额外的信息应优先选择能提供平滑性的吉洪诺夫正则化,而随着投影角度的增多(离散不适定问题占主导),截断奇异值分解法是首选,因为其能有效消除由小奇异值引起的误差放大。
解决欠定方程的第二种方法是尽可能减少未知量个数。BOS中常见的办法是添加掩膜(mask),即根据先验知识、获得的图像偏移测量值或已知流场的特性,来识别待测空间中折射率已知的体素(一般为折射率等于环境折射率的体素),进而约束实际未知量的体积,将解限制在期望区域内。文献[19]对比了添加和不添加掩膜前后流场的重构结果。图 12为有、无掩膜下的重构结果,其中左图为密度场三维等值图,中间的图为z向密度场切片,右图为x向密度场切片。从图 12可以看出,添加掩膜不仅降低了计算复杂性、减少了由稀疏相机布置引起的条状伪影,还提高了重建结果的数值准确性。类似的结论也在文献[65,89]得到证实。
图 12 有、无掩膜下的重构结果[19]
Fig.12 Reconstruction results with or without mask[19]
图 13(a)以小孔成像模型为例解释了掩膜添加原理。光线穿过待测流场的相关体素后投影在相机传感器上。若该光线上流场无梯度变化(绿线),那么相机传感器上该像素点的图案不发生变化,此时图像不发生偏移。因此可认为投影在不发生图像偏移的像素点上的相关流场体素,其折射率与环境折射率保持一致。对多相机测量系统,其掩膜添加的一般流程如下。
1)通过阈值法[65,89]或用户自定义[19,24],对每个投影面上的图像偏移结果划分二维掩膜区域,如图 13(b)黑色部分所示。对阈值法而言[65,89],二维掩膜的判定准则一般为图像偏移量小于0.1像素或者更小。澳大利亚蒙纳士大学Amjad等[89]在此基础上还增加了10%的非掩膜余量,以修正阈值判断误差、确保待测流场落在非掩膜区域。加拿大英属哥伦比亚大学Atcheson等[72]还指出不能简单地通过图像偏移阈值将流场投影与静态背景分割,因为即使图像偏移为0,三维流场梯度也可能不为0。例如当光线与折射率梯度平行时就会发生该情况。因此,文献[72]先对图像偏移场进行泊松积分,在求解折射率投影后,根据折射率投影的阈值来分割掩膜。
图 13 掩膜说明[19,90]
Fig.13 Mask illustration[19,90]
2)设定支持相应体素设为掩膜的最小投影数量。由于空间中流场任意体素都会投影到各相机传感器上,只有当该体素在超过给定数量的投影面上均判定为非掩膜,才能保留体素,以此构建三维掩膜。美国奥本大学Davis等[65]判定三维非掩膜成立的阈值条件为:穿过特定体素的所有光线中,需要有超过80%的光线图像偏移大于阈值,才能将该体素判定为包含在期望区域内。此外,该团队还考虑了单 视角下不同焦深和角度处理出的偏移图像, 现随着图像增多,其平均流场的三维掩膜结果更接近真实,其重构解也更优。
需要注意的是,掩膜创建完成后还需按照图4(b)准则判断二维非掩膜区域内的光线有效性,剔除无效光线及其相关方程组[19]。为了避免阈值法中掩膜结果受测量数据噪声的影响,德国杜伊斯堡-埃森大学Unterberger和Mohri[10]将进化式的遗传算法与BOS三维重构结合,提出在迭代过程中采用随机的方式逐步进化掩膜,进而使所建立的三维掩膜不依赖阈值。
减少未知量个数的另一种方式为利用一些先验信息或基函数对未知量做降维处理。如文献[20]利用傅里叶-贝塞尔基函数,在单视角下重构了层流非完全对称单火焰和多火焰的三维温度场,见图 14中模型UA-N和MA-N。不仅如此,在假定流场非对称性沿光线方向均匀的基础上,利用奇偶函数变换将流场分解为轴对称和非轴对称两部分,重构时仅需求解轴对称流场的未知量,见图 14中模型UA-Y和MA-Y。以上重构结果与多视角下的重构结果相比,还原度均达到90%以上。日本东北大学Lee等[79]利用方位角傅里叶基分解了待测流场,实现用7个相机的投影数据重构喷嘴压力比为2.42的欠膨胀射流三维密度场。更重要的是,文献[79]利用方位角傅里叶基还分解了具有高时间分辨率的麦克风数据,建立了BOS数据分解信息和麦克风数据分解信息之间的关系,仅利用非时间分辨的BOS层析和时间分辨的麦克风测量了具有时间分辨率的三维流场,实现了时空超分辨重构。其求解思路如图 15所示。
图 14 基于降维处理的单视角重构非对称流场案例[20]
Fig.14 Single-view reconstruction of asymmetric flow fields using dimensionality reduction-based techniques[20]
图 15 基于方位角傅里叶基的超时空分辨BOS层析测量[79]
Fig.15 Azimuthal fourier-based super spatiotemporal resolution BOS tomographic measurement[79]
正交分解除了可以对待测流场做降维处理外,也能对图像偏移场做降维处理,进而在流场重构时减少投影误差。法国航空航天实验室Lanzillotta等[85]利用图像偏移场的前两个正交分解模态信息,完成不同喷嘴压力比下超音速过膨胀射流脉动密度场的三维重构。
事实上,这种降维处理的方式类似于压缩感知(compressed sensing, CS):若信号在某个正交空间具有稀疏性,则可用少量测量值精确地恢复出原始信号[91]。笔者发现当未知量降维时,BOS层析求解所需的光线数量可以大大降低,即在减少计算成本的同时可以用少量的测量值完成流场重构[20]。用CS理论解决不完全角度引起的欠定流场重建问题时,其核心思想是在重建过程中增加关于稀疏变换后流场的稀疏性最小化约束,表达式为
式中:m为范数阶数或TV范数; x为要被重建的原信号; s为非零项远小于零项的稀疏信号; Ψ为投影矩阵。将非稀疏信号投影到稀疏性空间,找到合适的Ψ至关重要,常见的选择傅里叶变换、离散余弦变换、离散小波变换、字典学习等。目前CS理论主要应用在CT不完全投影数据重建中[91],特别是以图像全变分最小化为代表的一系列扩展算法。类似的算法也出现在其他光学层析重构技术中[92]。
无论是增加未知量约束还是减少未知量个数,都属于解决欠定病态方程的建模方法。在此基础上,需要合适的算法求解方程。目前公认的看法是代数重建算法(algebraic reconstruction technique, ART)比FBP更适合处理不完全投影重建问题,因为ART对噪声和不完整数据的鲁棒性更强[93],因而在BOS三维重构中得到了广泛应用[6]。ART的基本思想是基于投影数据和待测流场间的相互关系,利用迭代过程逐光线改进对流场的估计。其迭代表达式为
假定BOS病态矩阵Ax=B,,B∈,x∈中包含m个光线方程和n个未知量。式(17)中i为迭代次数k除以m的余数加1; ai为矩阵A的第i行; bi为观测向量B的第i个分量; λk为可选松弛系数,控制迭代的收敛速度和收敛性,且0<λk≤2。较大的松弛因子可以加速收敛但可能导致收敛发散; 较小的松弛因子可以保证收敛但收敛速度慢。
澳大利亚蒙纳士大学[32]、上海交通大学[37]、法国-德国圣路易斯研究所[94]采用ART算法完成了BOS三维流场重构。考虑到ART每次更新只用到一个投影方程,美国宾夕法尼亚州立大学[14,21,95]、东南大学[86]、美国NASA[24]在BOS三维重构时采用了联合迭代重建算法(simultaneous iterative reconstruction technique, SIRT),在每次迭代中考虑所有的投影信息为
xk+1=xk+λkTATM(B-Axk) (18)
式中T和M为对称正定矩阵。不同方法下T和M的选择不同,包括Landweber法、Cimmino法、组分平均法(component averaging, CAV)、对角放松正交投影法(diagonally relaxed orthogonal projections, DROP)、联合代数重建法(simultaneous algebraic reconstruction technique, SART)等,具体公式可见文献[96]。
澳大利亚蒙纳士大学[89]以360°范围内等间距布置15个相机的BOS测量系统为例,通过数值仿真测试了3种重建算法,即FBP、ART和顺序FBP-ART(将FBP作为初始值代入ART)。结果表明ART和FBP-ART求解精度相似,且误差都比FBP小50%。上海交通大学[16]用8个相机的投影数据研究火焰密度场时发现相较SIRT算法中的CAT和Landweber法,ART在重构中取得的结果更好。
除了代表性的ART和SIRT迭代方法外,可应用在BOS层析中的线性方程求解器还包括奇异值分解法(singular value decomposition)、最小二乘法(least squares method)、共轭梯度法(conjugate gradient, CG)、分裂布雷格曼法(split Bregman method)等[59,93]。CG是另一种常见的BOS三维重构算法[85,97-99],具体的算法流程见文献[100]。加拿大英属哥伦比亚大学[60,72,26]采用共轭梯度最小二乘法(conjugate gradient linear square, CGLS)完成了三维流场重构,并且指出当投影数量不多时可以用CGLS解决线性系统,但当投影数量超过8时,由于系数矩阵A所需的内存量较大,可以用SART算法,以矩阵-向量乘积的形式对矩阵做分块处理。法国航空航天实验室Nicolas等[19]同样利用矩阵-向量乘积的形式在GPU上实现了CG算法,可以在几小时内从超过1亿个投影数据中重建出超过1.3亿体素的三维流场。此外,该文献指出不完全数据的层析问题中迭代类算法需要在算法收敛前停止迭代,以避免噪声放大。德国杜伊斯堡-埃森大学[10]提出上述算法在解决非线性层析问题时(如大密度梯度流场重构)会遇到困难,因而提出了进化式的遗传算法来求解非线性重构问题。
需要注意的是,上述算法的求解对象均是式(11)所示无约束线性方程。即使求解对象是有约束的方程式(16),若想利用上述算法对其求解,则需要先将其按拉格朗日乘子法转为无约束方程。鲜有BOS将重构问题直接看成有约束方程求解。而在不完全投影CT测量时常将重构问题转换为有约束方程,并利用ART-TV、TVM-POCS(全变差最小化-凸集投影)等算法求解[91,101]。
近年来,基于数据驱动的深度学习模型凭借其强大的特征提取和非线性映射能力在不完全投影的层析重建问题中取得了许多突破性的成功应用。文献[102]将基于深度学习的不完全投影CT重建方法划分成如图 16所示的5类:①利用深度学习对不完全投影下的重构结果去伪影; ②利用深度学习补充缺失的投影数据,形成完全投影下的重构结果; ③利用深度学习先补充不完全投影数据,而后对其去伪影,实现投影域和空间域的双域联合处理; ④将深度学习嵌入迭代类算法中,代替原先的正则项与平衡参数的更新过程,或对目标优化方程中的子问题求解; ⑤利用深度学习直接从不完全投影数据映射得到重构结果。尽管深度学习模型在不完全投影CT重建中得到了长足的发展,但其在BOS三维重建中的应用仍处于起步阶段。现有基于深度学习的BOS层析测量方法主要集中在第5类端到端映射重建。
图 16 基于深度学习的不完全投影CT重建方法分类[102]
Fig.16 Classification of incomplete projection CT reconstruction methods based on deep learning[102]
第一篇利用深度学习解决BOS三维重构的研究为文献[68],该文献利用物理信息神经网络(physics-informed neural network, PINN)整合流场的基本物理规律与实验结果,实现仅根据时间和空间信息就能推断完整连续的三维温度、速度和压力场。该网络由全连接网络和残差网络组成:全连接网络用于联系时间、空间与温度、速度、压力场间的关系; 残差网络整合不可压缩Navier-Stokes和热传递方程完成温度、速度、压力场解的物理约束。整个网络的训练通过最小化物理约束残差项和温度预测不匹配项完成,如图 17所示。需要注意的是,温度预测不匹配是指由PINN预测的温度场与BOS三维重构测量结果间的差距,而BOS三维重构是通过经典的正则化迭代算法求得。因此,尽管该模型在预测时可以不受限于不完全投影数据,但在训练期仍依赖准确的BOS三维重构结果。换言之,由不完全投影数据引起的重构误差会被带进PINN中,造成网络预测也存在伪影、不适定等问题。类似的PINN网络也被德国柏林工业大学Rohlfs等采用,实现马赫数为2的湍流激波边界层相互作用的速度场重构[103]。
利用经典迭代算法求得的重构结果作为训练集,实现端到端式预测重构的深度学习模型还应用在化学发光层析测量中。如上海交通大学[104]结合卷积神经网络(convolution neural network, CNN)和长短期记忆网络(long short-term memory, LSTM)模型[见图 18(a)],通过前10个时刻的9个投影面数据来预测下一个时刻的火焰三维锋面。实验表明该模型能在毫秒级的时间尺度内准确预测层流扩散火焰和非预混湍流涡流稳定火焰的三维演变过程。在此基础上,该文献还提出了两种CNN框架(记为VT-Net1和VT-Net2),证明了当投影数量小于7时,CNN重构结果优于传统的ART [105][见图 18(b)]。南京理工大学分别采用CNN[106]和考虑正则化项的残差网络[107-108]从实时捕获的投影中可靠地重建出三维火焰结构,并根据循环神经网络框架(recurrent neural network),利用相邻投影数据间的空间相关性,提出一种基于门控循环单元(gated recurrent unit, GRU)神经网络的金字塔型BOS层析重建模型[109],如图 19所示。该模型利用残差网络提取投影数据的多层特征,从稀疏提取的特征开始,把二维特征图反投影到三维网格上,然后通过三维稀疏卷积输入到GRU网络。GRU整合来自不同角度的全局先验知识,所获得的稀疏流场作为下一个GRU的输入,最终得到稠密的三维预测结果。结果表明,该模型在重建速度和准确性方面优于传统的迭代重建法。
图 17 PINN-BOS模型[68]
Fig.17 PINN-BOS model[68]
图 18 CNN-LSTM模型与不同投影数下的算法重构性能对比[104-105]
Fig.18 Comparison of CNN-LSTM model and reconstruction performance under different numbers of projections[104-105]
图 19 基于GRU的BOS层析重构模型[109]
Fig.19 GRU based tomographic BOS model[109]
如前所述,上述模型的训练均依赖完全投影数据下的迭代重建结果作为真实值。因此不完全投影引起的以及迭代重建算法本身引入的噪声是不可避免的,这成为上述模型预测准确性和泛化能力的主要限制。为了不受该限制,美国宾夕法尼亚州立大学[63,69]将PINN模型中重构数据不匹配损失改为投影数据不匹配损失,并增加边界惩罚L=‖Ax^-B‖2),并增加边界惩罚b>bound,既减少了计算成本又规避了迭代重建的影响还能不依赖边界条件(见图 20),结合不同的物理约束方程(如Euler方程、无旋方程、可压RANS方程)实现了超音速流动、湍流膨胀射流的三维速度、密度、压力场的同时测量。针对PINN模型在存在高水平噪声时容易出现半收敛,提出了贝叶斯PINN来解决这个问题[110]。
图 20 PINN模型[63]
Fig.20 PINN model[63]
除了从算法上解决不完全角度BOS层析问题外,也可从硬件角度优化测量系统。光场相机利用微透镜阵列或其他特殊的光学配置来捕捉光线传播的更多信息,包括不同深度、角度和位置上的光线强度和颜色信息。图 21(a)和图 21(b)展示了光场相机的光学配置和成像原理[111-112]。可以看出,空间中任意一点在光场相机成像时,会被记录成一组具有不同角度和位置信息的微像素。每个微像素捕捉了特定方向和位置上的光线,组成光场。这种光场信息允许在后期处理时构建出不同焦深和角度的图像,从而实现对图像的动态调整、聚焦和视角改变,如图 21(c)所示。
图 21 光场相机成像原理[111-112]
Fig.21 Imaging principle of a plenoptic camera[111-112]
美国奥本大学[111]于2017年首次将光场相机应用在BOS图像偏移场的测量中,实现了单视角下同时捕捉不同深度位置的密度扰动。之后,文献[65,114]比较了光场相机和传统相机对图像偏移测量的影响,发现相较于传统BOS,光场BOS的优势在于扩展的景深和高效的光线收集。然而,这些优势是以信噪比降低了80%倍、像素尺寸比传统系统大了14倍(基于微透镜直径)以及空间分辨率降低引起的稀疏测量为代价,其结果如图 22(a)~图 22(b)所示。而传统BOS因较小的像素尺寸具有高空间分辨率,但由于其景深小,只对特定位置敏感。并且传统BOS还由于光圈较小而需要更长的曝光时间。长曝光时间可能导致传统BOS图像无法显示流动中的某些瞬态特征。采用4个光场相机重构了两个圆柱状的液体三维折射率,如图 22(c)~图 22(d)所示。实验表明:即使是光场相机,单相机中的有限角度范围也无法满足层析重构要在较大角度范围内进行投影的数据需求; 从各光场相机中渲染出不同焦深和角度下的多个透视视图有利于提升重构质量并减少伪影。因此当相机数量有限或光学访问受限时,光场相机的使用可能优于相同数量的传统相机。
图 22 光场BOS结果与系统[65,114]
Fig.22 Results of plenoptic BOS and systems[65,114]
显然提升光场相机的测量性能是扩大光场BOS应用范围的关键。美国桑迪亚国家实验室Guildenbecher等[113]开发了一款高倍率、远工作距离的光场BOS系统,可在超过0.25 m的工作距离下实现光场5倍的放大倍率,形成具有约25 μm的空间分辨率和三维成像能力。该技术被用于对马赫数为3.3的自由空气射流中的激波结构进行三维成像的展示。
在光学访问受限的情况下,相机硬件布置方式,包括相机数量、相机分布与拍摄角度,也会影响BOS层析系统测量的准确性。法国航空航天实验室Nicolas等[19]比较了3种相机布置方案,发现:①对于共面型布置[见图 23(a)],相机越多重建效果越好; ②当相机布置于垂直流动主轴的平面上时,相机共面的重建效果是最佳的,但不共面的分散型布置也不会显著改变重建质量[见图 23(b)]; ③当相机视角过于局限时,重建效果非常差[见图 23(c)]。
图 23 3种相机布置方案[19]
Fig.23 Three camera layout configurations[19]
德国柏林工业大学Lang等[55]进一步研究了共面型布置下的相机个数和视角,如图 24所示。发现:①投影数量需超过10个才能保证重构结果的准确性; ②投影为全圆型奇数布置时,重建误差小且收敛速度快; ③当投影数量小于8个,半圆型布局的误差与全圆型奇数布局的误差一致; ④半圆型布局会引起伪影的不对称性以及重建流场中心的轻微不对称性。
图 24 6种相机布置方案[55]
Fig.24 Six camera layout configurations[55]
在其他光学层析研究中,上海交通大学Wang等[115]提出了一种基于最小化各视角相关性的相机优化布置算法,并发现当光学访问存在障碍物遮挡导致共面型布局无法布置在理想平面时,以图 25所示的4π空间型布置能得到更优的火焰结构。类似的结论得到瑞典隆德大学[116]证实,发现当投影数量超过6个时,任意相机布置的重建质量差于共面半圆布置方案,反之则任意布置方案更优。西安工业大学Wang等[117]发现实现层析重建所需的最小投影数取决于摄像机像素的数量和大小,以及重建区域中离散体素的数量和大小,并通过Mojette变换的精确重建条件来确定最小投影数。北京理工大学Gao等[118]提出了一种同时考虑相机角度和流场分布影响的新布局优化方法,并通过层流和湍流火焰锋面的三维重建进行了验证。南京理工大学荣韶华[107-108]指出在有限光学窗口测量条件下,在窗口内尽可能布置多的阵列相机也能提升重构质量,如图 26所示。
图 25 两种相机布置方案[115]
Fig.25 Two camera layout configurations[115]
图 26 相机布置方案[107]
Fig.26 A camera layout configuration[107]
本文介绍了BOS成像及层析测量的基本原理和硬件设备,分析了当前BOS层析测量中面临的不完全角度问题。从数学角度看,投影数据的缺失、巨大的方程规模、不可避免的观测值噪声以及理想假设下系数矩阵与实际的差异性,使得BOS层析重构方程变得欠定和病态,这对重构算法的高效性、正确性和鲁棒性提出了挑战。本文围绕增加未知量约束、减少未知量个数、求解算法、深度学习、增加图像信息、优化有限角布置6个主题,从软件和硬件两部分回顾了近年来所发展的不完全角度BOS层析重建算法。可以看到,正则化约束的引入、三维掩膜的施加、未知量的降维处理、无约束型代数迭代算法的使用、端到端式神经网络的深度学习以及光场相机的多维信息采集和多相机布局优化等工作都有效地解决了不完全角度下BOS重建后流场质量降低的问题。但现阶段不完全投影层析技术在BOS重构中的应用深度和广度还远远不够。可预见,利用压缩感知和其他非端到端式深度学习等理论,结合丰富的先验信息和其他层析技术,可以从去伪影、增投影、减未知、优化迭代、优化布局全方位为不完全角度BOS重建问题的解决提供新思路。
除不完全角度问题外,BOS重构还面临视场不能覆盖流场全部空间的内问题和流场内部存在光线无法穿透物体的外问题。如何求解这两类问题引起的欠定病态方程,系统地提高不完全投影数据下BOS的重建质量,是亟待解决的应用难点。总之,经过20余年的发展,BOS已从可视化工具慢慢发展成三维、定量的折射率、密度、温度、速度测量技术,可应用于不同的流动环境。但在大视野、高频动态、低成本、好测量的优势背后,仍要着眼于其在高时空分辨率、高精度、高敏感度、高实时性等测量性能的挑战,以推动背景纹影技术的发展和产业化落地,更好地服务于基础和应用研究。