碎片快速撞击网格剖分重构

本文作者：葛双全 肖永浩 赵士操 单位：中国工程物理研究院 计算机应用研究所 仿真中心

0引言

地面超高速碰撞模拟试验是研究航天器抗空间碎片撞击性能的主要方法之一，撞击试验需要配备气体炮、电磁炮等超高速发射装置和摄像机系统、测量系统等辅助设备．为获取碰撞后产生的碎片云和防护材料背面喷射／表面溅射（又称为反溅）效果，国内外通常使用X光高速照相技术［1－2］或激光阴影照相技术［3］来记录碰撞过程．而脉冲激光全息技术用于研究强冲击载荷条件下靶板材料粒子的喷射／溅射现象，具有图像直观、信息量大、非接触和抗干扰等优点．尤其在同轴Fraunhofer激光全息技术［4－5］中，采用超短脉冲激光作为记录光源，可以在粒子场演化过程的某一时刻将其“固定”，从而达到测量的目的．通过激光全息照相获得碎片撞击防护结构表面材料喷射和（或）溅射粒子场全息图像，模拟其演化过程，进而获取粒子的各项参量，对于研究航天器防护结构及其材料在超高速撞击条件下的损伤破坏机理，实现航天器防护结构的优化设计［6－8］具有重要的工程价值．本文基于碎片超高速撞击防护结构所产生的粒子场激光全息图像研究撞击过程三维重构技术，对粒子堆叠区域网格剖分、粒子轮廓提取和粒子场演化过程模拟等关键步骤进行了详细的阐述，为研究空间碎片对航天器防护结构的损伤效应提供了一种新的思路．

1激光全息图像预处理

脉冲激光全息图像实际上是参考光与衍射物光形成的干涉条纹，需建立相应的同轴全息再现系统［9］，得到再现图像后才能进一步处理，本文所述激光全息图像均为全息再现图像．全息图像预处理包括边缘剪切、缩放和子图像分割三个步骤．边缘剪切可以去除图像周边的噪音区域；缩放是指根据工作站显示和运算性能调整图像显示比例，使之在可视化效果和系统运行效率之间取得平衡；子图像分割将缩放后的图像再划分为若干子图像．划分子图像的原因有两个：其一，全息图像既有粒子分散甚至空白区域，又有大量粒子堆叠的区域，考虑到在粒子场三维重构过程中不同区域的粒子在三维空间中形态各异，必须针对粒子的不同形态分别划分二维空间区域，对应在全息图像上就是大小不同的子图像；其二，如果单个图像尺寸过大，则粒子堆叠区域网格剖分和粒子轮廓提取所消耗的时间将非线性大幅度增加．经过多次比对，将单个子图像分辨率保持在300×300左右将是一个比较合适的尺寸．

2粒子场堆叠区域网格剖分

碎片超高速撞击防护结构所形成的表面材料粒子场密度和粒子与防护屏表面之间的距离成反比，越靠近防护屏表面，粒子密度越大，堆叠程度越高，当密度达到一定程度后，全息图像显示为大片的白色区域［10］（如图1左侧边缘）．显然，这些白色区域并非单一粒子投影而成，必须对其进行网格剖分，剖分得到的每一个多边形子区域表示为单个粒子．粒子堆叠区域网格剖分由粒子平均面积和粒子场密度决定，对于由稀疏区域和堆叠区域所组成的粒子混合区域，可参考稀疏区域的粒子平均面积和密度设定参量．对于子图像所有区域均为堆叠区域的参量数值，可参考临近稀疏区域或混合区域的面积与密度参量，也可根据撞击试验数值仿真分析结果设定．将Ω划分为凸子区域的方法是将其划分为三角形，可借助二维Delaunay三角形网格生成算法，二维Delaunay三角化的定义为任意点集PS中所有顶点对应的Delaunay三角形的并集是对点集的一个三角形剖分，该算法较为成熟，本系统采用Triangle［11］．在实现过程中需要判断空间点是否在区域内部，首先将Ω划分为若干个凸子区域C＝｛Ci｝，其中Ci由线段表示，以顺时针给出，则可以根据给定点及点到线段的两条直线的叉积来确定点是否在区域内部．假设构成Ci的某线段起始点a的坐标为（ax，ay），终点b的坐标为（bx，by），给定点h的坐标为（hx，hy），此时计算叉积。通过将ε中的正方形通过组合形成ω，根据值k可以确定由多少个εi组成一个ωi，然后这些εi中属于区域内部的那些正方形的外边界构成ωi的边界．由于构成ωi的εi集合不一定是连通的，因此，首先将要将连通的εi组成一个集合，然后再找这些集合的外包络边界．形成连通集合ξ时采用贪婪算法，其过程为：找到一个出发点，将其加入待搜索集σ；将σ加入到ξ；搜索所有属于σ的正方形相邻的正方形集合；将属于ε，但不属于ξ的正方形加入领域集η；若η＝φ，则搜索结束，否则将σ＝η，转至b；得到ξ，可获得其外包络边界，则这个外边界构成所需的ωi．图2给出一个堆叠区域划分示例．最终结果对原始区域边界的逼近度取决于图2（c）中规则化的细分程度，细分程度越高，对边界拟合度越好，但需要花费更长的时间．

3粒子轮廓提取

粒子图像轮廓实质上是粒子边缘所组成的封闭区域，图像信号的一阶微分有最大值，二阶微分过零点即为图像边缘．梯度对应于一阶导数，基于梯度运算的算子是一阶微分算子．对于一个连续函数f（x，y），它在位置（x，y）的梯度可表示为一个矢量（两个分量分别为沿x／y方向的一阶导数）由于数字图像是离散的，计算偏导数Gx和Gy时，可用差分代替微分．为计算简便，采用小区域模板或图像卷积来近似计算差分．采用不同的模板计算Gx和Gy可产生不同的边缘检测算子．由于Sobel结合了Gaussian平滑和微分，对图像噪音具有一定的鲁棒性，因此采用Sobel算子计算图像差分．对于每一个轮廓，可通过轮廓逼近的方法减少轮廓描述点，从而达到减少三维形体三角面片的目的。4粒子场三维空间投影经过网格化和粒子轮廓提取后，得到一系列单一粒子的二维轮廓点坐标，还需要将其投影到矩形、锥形或柱形等特定的三维投影空间中去．图

4为锥形空间投影示意．设锥形投影空间高度为H，宽度W，在某一XZ圆锥截面上有两个粒子P1、P2，其中心坐标分别为（x1，y0）和（x2，y0），通过式（9）可计算其圆锥截面半径如果粒子横坐标x小于等于r（P1），则表示该粒子位于圆锥体内，则其z坐标为对于每一个轮廓，可通过轮廓逼近的方法减少轮廓描述点，从而达到减少三维形体三角面片的目的化示意．

5粒子场演化过程模拟

经过三维空间投影后，得到每一个粒子的轮廓点数及其三维空间坐标．可通过3ds　MAXScript语言在3ds　max中读取、导入和显示粒子场，并构造每一个粒子的运动轨迹，模拟粒子场演化过程．获得每一个粒子的轮廓点数及其三维空间坐标，该坐标实际上仅为其在XY平面上的投影，如果要构成三维形体，还必须赋予粒子一定的厚度．首先构造由轮廓点连接而成的闭合区域，然后计算前两个轮廓点之间的距离作为粒子三维形体的厚度．如果粒子轮廓点数大于等于2，则类推，如果仅有一个轮廓点，则构造单位立方体．在3ds　max中，由二维平面轮廓构造三维形体通过挤出（Extrude）效果实现．粒子运动轨迹需要用户输入起始和最终运动位置，每一个粒子从起始位置按设定的运动轨迹移动到最终位置的过程即为粒子场演化过程．可设置粒子旋转等效果．

6验证

图6（a）为零时图像，图6（b）为摄像系统所拍摄的某时刻脉冲激光全息图像．图7为该粒子场演化过程模拟结果，当t＝100时，与拍摄时刻的全息图像吻合较好．

7结论

本文基于同轴激光全息照相系统拍摄的防护结构喷射和（或）溅射粒子场全息图像，对超高速撞击防护结构粒子场进行了三维重构技术研究，并对重构方案及若干关键步骤进行了详细阐述．重构结果表明，无论是粒子场静态三维重构结果还是其演化过程均与撞击试验所拍摄的全息图像吻合较好，验证了该技术的有效性．该技术存在的主要问题是高分辨率全息图像（20　000×10　000）在三维重构之后生成的粒子数量极大，从而导致其演化模拟和可视化显示对计算机提出了极高的要求．通过调整全息图像大小，减少生成的粒子数量可以缓解此问题，但同时势必降低可视化效果．因此，建议系统运行在拥有较高内存和主频的64位机器上，同时需要综合可视化效果和系统运行效能之间的平衡．