Visual monitoring method for atmospheric disturbance of moving objects

返回式航天器、民航客机、降落伞等是在大气层内飞行的快速运动物体，在其运动过程中都会与大气相互作用，形成大气扰流，其是在大气层内动力飞行的必然产物。对这些运动物体引发的大气扰流进行可视化光学监测可有效获取流动参数、监测运动状态等，可为返回式航天器平稳着陆提供数据支持，同时也可为民航客机超音速飞行、降落伞动力开伞的动力外形优化提供支撑。所以，大气层中运动物体大气扰流的可视化光学监测方法是动力学研究一项重要方向。

密度分布是压缩流体流动显示的主要因素，为了实现气体流动可视化，A. Toepler于1864年最早提出了纹影方法，沿用至今。但纹影成像方法需要复杂、精密的光学系统以及明亮的光源，利用折射光线揭示测试目标周围的空气密度梯度。1999年，Meier^[1]提出了背景纹影成像（Backgroud Oriented Schlieren，BOS）技术，具有专利权。2001年，Raffel^[2]对背景纹影成像技术进行精确描述，并在十几年的流动测量方面取得较大进展。背景纹影成像利用随机图案的背景以及成像相机来显示背景前面区域密度变化。相对于纹影成像技术，背景纹影成像技术取得了更广泛的应用，从风洞内拓展到实验室外。2010年，Michael^[3]等人利用高帧频相机以树林为背景开展了喷雾羽流、引擎盖热流的可视化成像，并以玉米地为背景开展了步枪射击、爆炸形成大气扰流的可视化成像，阐述了背景选择、成像设备、成像距离、背景纹影探测精度等与所拍摄运动过程之间的相互关系。但这些基本上应用在外场近距离、超高分辨率、水平观测方面。2013年，Raffel^[4]、Edward^[5]等人以树林为背景，对飞行中的全尺寸直升飞机旋翼尖端涡旋进行了可视化成像，并与大尺寸风洞内试验结果进行了对比分析，确认了背景纹影成像可对大气扰流进行清晰、精准成像。直升机飞行速度较慢、飞行速度较低，NASA将其应用到远距离、地物为背景的超音速冲击波显示领域。2011年，NASA在2014~2015、2019年AirBOS1试验基础上，开展了多次空对空飞行试验，以地物为背景分别获取了单架^[6-7]、双架^[8]超音速飞机大气扰流的图像，实现了冲击波相互作用过程的可视化。2017年，NASA开展了钙-K日食背景纹影^[9]（CaKEBOS）技术论证，在地面以日盘为背景，获取了飞机进出太阳边缘的大气扰流图像。这些重要技术进展都使得纹影图像质量取得重大改进，且每次技术进展都有独特的应用潜力。但是NASA并未公开相关技术细节，对相关技术研究有着广泛的应用价值。

文中以背景纹影成像原理为基础，建立了对运动物体在大气层内运动引发的大气扰流进行光学监测的方法。监测方法包括大气扰流光学监测机理、大气扰流光传输方法、扰动检出方法三个部分。适用于对在真实飞行条件下的返回式航天器、民航客机、降落伞等的大气扰流进行高精度大气扰流显示。

大气扰流光学监测主要基于平面波经过湍流流场产生波前畸变的原理，如图1所示。当光束通过随机流场时，光束受到某种方式的扰动，与无扰动时相比，光线会偏离原来的传播方向。

Figure 1.  Schematic diagram of wavefront distortion generated by a plane wave passing through turbulent flow field

入射光线平面波S在折射率为n₀(x,y,z,t)的无扰动流场传播，穿过流动速度为V_c、折射率为n₁(x,y,z,t)的湍流场后，由于流场中各种漩涡对光线的作用，出流场时，光线路径变化为S´，不仅位置发生移动，而且在畸变的波前位置光线方向改变了ε。

假设厚度为L的湍流层对通过它的波前S(x,y)=1（波长为λ）引入相位畸变，产生畸变的波前S´(x,y)。

上式为一阶几何光学近似（假定│n₁(x,y,z,t)│<<1），（x，y）为孔径坐标，z为光波传播方向坐标，k为波数（2π/λ），n为局地折射率，n´(x,y,z,t)为n₁(x,y,z,t)与n₀(x,y,z,t)的差值，OPL为光线通过湍流层路径的光学长度，ΔOPL为光线通过湍流层路径偏差。空气密度ρ(x,y,z,t)与ΔOPL满足下式：

式中：K_G-D为格拉斯通−戴尔常数。因此，波前畸变与光线传播方向的空气密度相关。

利用大气扰流引起大气密度发生变化，从而引起大气扰流区折射率变化的原理，通过探测地物光线通过气流密度变化区域造成的图像变化，检测波前变化，从而实现大气扰流的可视化光学监测。

大气扰流光学监测方法主要包括以下3步：

（1）利用目标大气扰流光传输方法，获取目标大气扰流光传输特性；

（2）利用大气扰流光偏折监测方法，获取包含大气扰流信息的背景图像；

（3）利用高精度扰流检出方法，获取可视化的目标大气扰流。

如图2所示，大气扰流光传输方法，利用结构设计软件，建立运动物体的精细三维外形模型，对运动物体及大气扰流区进行计算网格划分。利用流场分析软件，设定运动物体飞行高度、飞行速度、飞行姿态等参数，以及大气环境参数，选定流场计算模型，以获取运动物体引发的大气扰流区的密度、压力、温度等参数的三维分布。此参数计算方法较为成熟，各类相关书籍及文献中均有详细介绍，文中不再赘述。利用格拉斯通−戴尔公式，建立大气扰流场内各点折射率与该点密度、压力、温度等参数之间的关系。再利用下面所述的流场折射率、折射率梯度以及扰流场光线传输路径计算方法，获得大气扰流光传输特性。目前国内暂未开展相关研究。

Figure 2.  Composition of method for optical transmission in atmospheric disturbance

本节主要介绍大气扰流场折射率、折射率梯度以及扰流场光线偏折传输路径计算方法。设定为以P（x,y,z）为球心、r为半径的球体为流场区域，采用对P点附近、在球体范围的各点折射率进行距离加权平均插值的方法，求解空间任一点P（x,y,z）的折射率n_p（x,y,z）。

式中：n_i（i=1,2，…）为球体内空间点上的折射率。根据大气扰动流场中任一点P（x,y,z）折射率值n_p(x,y,z)计算得到空间内任一点折射率梯度值。折射率梯度的一般求解公式为：

采用Barron算子利用该点周围流场中其他多个点，求解流场内任一点折射率梯度值。

以x方向为例，取点沿x方向上邻近流场四点（x_i-2，y_j，z_k），（x_i-1，y_j，z_k），（x_i+1，y_j，z_k）和（x_i+2，y_j，z_k）的折射率用Barron算子法进行插值，得到流场中任一点（x_i，y_j，z_k），沿（x_i-2，y_j，z_k）方向的梯度值表达式：

同理，可得到沿y和z方向的折射率梯度：

将折射率梯度代入下式求解x和y方向的偏折角。

利用上述方法，对同一运动物体在相同的飞行高度、飞行速度，以及相同大气背景条件下，不同谱段地物光线以相同入射角进入大气扰流区的传输特性进行仿真分析，仿真结果如图3所示。

Figure 3.  Simulation result of deflection characteristics of ground light with different optical wavelength passing through the same atmospheric disturbance area

由仿真结果可知，不同谱段（0.532、1.06、4.0、10.0 μm）地物光线经过大气扰流区后，光线偏折特性分布情况相似、同一区域偏折大小相似。所以，光线穿过运动物体大气扰流区的偏折特性对波长选择性不大，但这一结论有待试验验证。

大气扰流光偏折监测方法，以上节得到的大气扰流光传输特性作为输入，利用背景纹影成像原理，建立监测系统数学模型，仿真分析各系统参数对监测性能的影响。

由格拉斯通-戴尔定律可知，气体折射率与密度的关系可用下式表示：

式中：n为气体折射率;ρ为气体密度;K_G−D为格拉斯通−戴尔常数，取决于气体的特性，具体形式可参考文献[10]。

背景纹影成像原理如图4所示，Z_i为像平面到相机透镜组的距离，Z_B为背景到相机透镜组的距离，Z_D为背景到流场中心的距离，ε_y为光线偏折角，Δy为参考图像和实验图像上对应点之间的y方向相对位移，Δy´为Δy在背景上的虚位移。同理，Δx为参考图像和实验图像上对应点之间的x方向相对位移，Δx´为Δx在背景上的虚位移。

Figure 4.  Schematic diagram of background oriented schlieren

由于光束的偏折包含了空间折射率梯度场沿光程的积分效应，因此，背景图像斑点发出光线的偏折角可表示为：

式中:ΔZ_D为测量流场的半宽度。由原理图的几何关系可知，图像平面位移量Δy与虚拟图像平面位移量Δy′之间关系可表示为：

式中：f为镜头焦距。因此，对于小偏折角而言，可得偏折角表示如下：

则

同理，在x方向的位移△x为:

当光线在非均匀介质中传输时，根据费马原理，如果光线偏移量远远小于流场宽度，则有：

式中：C为常数，与实验配置有关；Δx, Δy为测得的斑点在不同方向的位移量。对整个位移矢量场x向和y向求偏导，则可获得如下泊松方程：

对于给定的位移矢量场以及给定的边界条件，上式可通过有限差分或有限元法求解，进而获得测量区域的投影积分效果的定量折射率分布，并通过格拉斯通−戴尔公式得出定量密度场信息，从而实现大气扰流场的可视化反演。相关模型较为成熟，但高空背景纹影成像各因素对成像效能的影响分析，暂未见公开报道。

利用上述方法，对监测系统焦距、探测器像元分辨率、图像处理能力等各系统参数对系统监测性能的影响进行仿真分析。

由大气扰流光偏折监测方法数学模型可知，监测系统可监测的光线偏折角大小与目标飞行高度及监测高度相关。假设监测系统焦距为0.6 m、探测器像元大小为5 μm，图像处理精度为1/10 pixel，则运动物体飞行高度与监测高度对系统监测性能的影响见图5。

由图5可知，在相同监测高度下，运动物体飞行高度越高，系统可监测的偏折角越小，系统监测能力越强；在相同运动物体飞行高度下，监测高度越低，系统可监测的偏折角越小，系统监测能力越强。

由大气扰流光偏折监测方法数学模型可知，在运动物体飞行高度及监测高度一定的情况下，监测系统可监测的偏折角大小与探测器性能及图像处理能力相关。假设运动物体高度为6 km、监测高度为12 km、监测系统焦距为0.6 m，则探测器性能及图像处理能力对监测系统监测性能的影响如图6所示。

Figure 5.  Relationship between target flight altitude and performance of monitoring system in different monitoring heights

Figure 6.  Relationship between image processing capability and performance of monitoring system with different detector pixel sizes

由图6可知，在采用相同探测器件情况下，图像处理能力越强，系统可监测的偏折角越小，系统监测能力越强；在相同图像处理能力情况下，探测器像元越小，系统可监测的偏折角越小，系统监测能力越强。

由大气扰流光偏折监测方法数学模型可知，在运动物体飞行高度及监测高度一定的情况下，监测系统可监测的偏折角大小与探测器性能、图像处理能力及光学成像系统性能相关。其中探测器性能与图像处理能力可以用像面位移识别能力来综合表示。所以，假设运动物体高度为6 km、监测高度为12 km，则光学系统成像性能及像面位移识别能力对监测系统监测性能的影响如图7所示。

Figure 7.  Relationship between focal length and performance of monitoring system with different recognition ability in image surface displacement

由图7可知，在光学成像系统焦距相同时，像面位移识别能力越高，系统可监测的偏折角越小，系统监测能力越强；在相同像面位移识别能力情况下，光学成像系统焦距越大，系统可监测的偏折角越小，系统监测能力越强。

如图8所示，高精度扰流检出方法利用大气扰流光偏折监测方法获取的一幅包含运动物体大气扰流的地物图像和一幅不包含运动物体大气扰流的同一地物图像，进行整像素十字搜素^[11]和Newton-Raphson亚像素定位^[12]实现高精度配准，再通过边界条件设定、网格划分等，对大气扰流光偏折监测方法中得到泊松方程进行求解，以实现大气扰流反演，从而获得可视化的目标大气扰流信息。

Figure 8.  Composition of method for high-precision disturbance monitoring

本节主要介绍整像素搜索和亚像素定位方法的数学模型，其他模型已在上节详述。整像素搜索采用的是十字搜素方法，如图9所示。十字搜索算法根据相关系数分布曲面具有的单峰特性，将相关系数分布从H维的曲面变为一条只有单峰的曲线，搜索路径从二维降到一维，在计算精度不降低的情况下，大大减少了计算时间。假设C点为起始搜索点，u_v是起始点沿u方向的坐标值，假设沿此点对H维的相关系数曲面作一平行v方向的截面，则该截面与相关系数曲面相交的曲线就是一条单峰曲线，其顶点e₀所在位置，就是沿直线计算的相关系数极值点。然后经过e₀并垂直上一截面作一新的截面，交v轴于v₁，此截面与相关系数曲面相交的曲线顶点为e₁。如此重复搜索几次，最后在单峰曲线的顶点e_n为同一点时，认为该点就是平面内的相关系数极值点，结合C点位置即可得到整像素位移。

Figure 9.  Schematic diagram of integral pixel search

相关系数为：

式中：M、N分别为求解的x和y方向的像素数；（x，y）为初始图像中某一像素位置；（xʹ，yʹ）为求解图像中任意迭代窗口位置；f（x，y）为初始图像中准备求取位移量的诊断窗口；g（xʹ，yʹ）为求解图像中任意迭代窗口。

亚像素定位采用的是Newton-Raphson法（N-R法），其核心是通过牛顿法在相关系数极值附近迭代，最终收敛至相关函数极值。假设变形后图像存在刚体位移、伸缩、剪切等变形。因此，相关函数的自变量设定为含有6个变形参数（u、ν、∂u/∂x，∂u/∂y，∂ν/∂x，∂ν/∂y）的矢量。变形后图像子区各点可用如下形函数表示：

进行迭代的相关函数为：

式中：$\vec p = \left[ {u,v,\dfrac{{\partial u}}{{\partial x}},\dfrac{{\partial u}}{{\partial y}},\dfrac{{\partial v}}{{\partial x}},\dfrac{{\partial v}}{{\partial y}}} \right] $为待求矢量。变形前后图像子区的相似度最高的时候。相关系数${C_{f,g}}(\vec p) $应取极值，即相关系数${C_{f,g}}(\vec p) $梯度趋于0，则：

上式便可用N-R法进行迭代求解，N-R法求解极值是基于泰勒级数的前两阶展开式。

为了验证高精度扰流检出方法计算精度，采用若干随即分布高斯光斑强度迭加的方法模拟发生位移前后的散斑图像，建立了可以精确控制位移的数学仿真散斑图，如图10所示。

Figure 10.  Speckle pattern in Gaussian distribution got by numerical simulation

设定图10所示的散斑图X、Y方向各0.01个像元的随机位移，采用高精度扰流检出方法进行计算，可得如图11所示计算结果。图中u₀，v₀为给定的沿X和Y方向的位移，左列为计算所得的各像素沿X方向的位移值分布，右列为计算所得的各像素沿Y方向的位移值分布。将各分布图采用正态分布N(μ,σ²)进行拟合，整理计算如表1所示。

Figure 11.  Given displacement in X direction u₀=0.01 and given displacement in Y direction ν₀=0.01

由表1可知，对于标准高斯分布的散斑图，本方法采用的位移搜索算法精度，理论上可达1/50 pixel。

文中在大气扰流光学监测机理介绍基础上，分别从数学模型、仿真分析两个方面，详细介绍了由大气扰流光传输方法、大气扰流光线偏折监测方法以及高精度扰动检出方法组成的大气层中运动物体大气扰流的可视化光学监测方法。

大气扰流光传输方法中详细介绍了大气扰流场任意一点折射率、折射率梯度以及光线偏折追痕的方法，并通过仿真分析得出光线穿过运动物体大气扰流区的偏折特性对波长选择性不大。但这一结论有待试验验证。大气扰流光线偏折监测方法中详细介绍了监测系统各参数对监测效能的影响规律。高精度扰动检出方法中详细介绍了“整像素十字搜索”+“N-R亚像素定位”的方法，理论上可实现1/50 pixel的大气扰流检出精度。

文中的成果为返回式航天器、降落伞、民航客机等是在大气层内飞行的快速运动物体形成的大气扰流可视化监测提供了一种重要的方法与途径。