基于旋转因子的曲面流场涡特征提取方法

摘要

本发明公开了一种基于旋转因子的曲面流场涡特征提取方法，目的是解决曲面流场涡特征快速、准确提取的问题。技术方案是首先通过顶点矢量数据的大小检测出矢量场中所有的临界点；然后根据临界点周围圆环的旋转系数确定临界点中哪些为涡核点；最后通过涡核点周围非临界顶点的旋转系数来确定涡区的范围。确定了流场中涡核点和涡区的范围也就完成了涡特征的提取。采用本发明可准确检测曲面涡特征的位置与范围，有效提高涡特征检测的准确率，可以有效去除鞍点区域。

说明书

技术领域

本发明涉及曲面流场涡特征提取方法，尤其是基于旋转因子的曲面流场涡特征提取方法。

背景技术

随着高性能计算机性能的飞速提升，科学计算产生的数据规模越来越大（TB-PB级）。面对如此海量数据，特征提取技术发挥了重要作用，该技术能够帮助研究人员迅速定位与提取数据场中的重要特征，剔除大量的冗余信息，准确保留数据的重要信息，使得用较小的存贮空间保存大规模数据的重要特征。

涡是流场中的重要物理特征，目前涡可视特征的主要提取方法有λ₂方法、Q-标准、M_z-标准、基于信息熵的特征提取等。这些方法不能够准确提取流场中的全部特征，容易出现错误提取的现象；同时，λ₂方法、Q-标准、M_z-标准等方法只适用于2维平面数据和3维体数据，对曲面数据无法进行特征提取；基于信息熵的特征提取适用于曲面特征提取但是存在严重的误检，这主要是因为该方法不能有效区分涡核区域和鞍点区域。

涡特征的定义最早于1972年提出，涡是指大量物质点围绕一个公共中心点的旋转运动状态。因此，涡特征包含涡核点和涡区两个方面，其中这一公共的中心点为涡核点，大量粒子所占的空间为涡区。在通常情况下，涡核点为流场中的临界点，其中临界点是指流场中速度为零的数据点，但并不是所有的临界点都是涡核点。

发明内容

本发明要解决的技术问题是：针对曲面流场提出一种基于旋转因子的曲面流场涡特征提取方法，解决三角网格曲面流场涡特征的准确提取问题。

本发明具体的技术方案是：首先，通过顶点矢量数据的大小检测出矢量场中所有的临界点；然后，根据临界点周围圆环的旋转系数确定临界点中哪些为涡核点；最后，通过涡核点周围非临界顶点的旋转系数来确定涡区的范围。确定了流场中涡核点和涡区的范围也就完成了涡特征的提取。

具体的技术方案是：

第一步：由曲面顶点数据vertexList[M][6]、网格数据faceList[N][3]，得到各个顶点的邻接面片索引vertexNeighberList[M]。其中vertexList[M][6]和faceList[N][3]为代表曲面流场的两个数组。vertexList[M][6]为曲面流场的顶点数据，M为曲面流场中的顶点数目，每一个顶点对应的6个浮点数分别为该顶点的三维坐标值和流场矢量的三个分量；faceList[N][3]为面片数据，N为曲面流场中的面片数目，每个面片对应3个整数，这三个整数为三角形面片三个顶点在vertexList[M][6]数组中的顶点编号；邻接面片索引数组vertexNeighberList[M]存储每个顶点的邻接面片链表（顶点的邻接面片是指以该点作为面片顶点的所有面片）的首指针，每个链表存储顶点的相邻面片，邻接面片链表的元素为描述面片信息的结构，其中包括面片号、链表中下一元素指针。获得vertexNeighberList[M]主要是为了方便以后的计算，如法向量计算等方面。得到各个顶点的邻接面片索引vertexNeighberList[M]的具体方法是：

1.1：初始化邻接面片索引数组vertexNeighberList[M]，将其元素均赋值为NULL。

1.2：初始化面片序号循环变量i₁等于0。

1.3：序号为i₁的面片是其顶点faceList[i₁][0],faceList[i₁][1],faceList[i₁][2]的邻接面片，因此，需要将序号为i₁的面片作为上述三个顶点的邻接面片分别插入到各自的邻接面片链表中。具体步骤为：对上述三个顶点的邻接面片链表分别进行遍历，如果三个顶点中任意一个顶点的邻接面片链表中有元素的面片号为i₁，则转到1.4步；如果不存在面片号为i₁的元素，则将序号为i₁的面片作为链表的最后一项插入到链表vertexNeighberList[faceList[i₁][x]]中，其中x为0、1或2，然后执行1.4步。

1.4：i₁增加1。

1.5：若i₁小于N则转到1.3步，否则执行第二步。

第二步：根据vertexList[M][6]、faceList[N][3]得到各个顶点的邻接顶点索引数组VNList[M]。VNList[M]中的元素为各个顶点的邻接顶点链表首指针，其中顶点的邻接顶点是指与该顶点有一条边（三角形面片中的边）相连的所有顶点，各个顶点对应的邻接顶点链表中存储该顶点的邻接顶点信息，邻接顶点链表的元素为描述顶点信息的结构，包括顶点号、下一邻接顶点指针。其中顶点邻接顶点是指与该顶点位于同一面片内的顶点，获得该数组主要是为了方便以后的计算。具体方法是：

2.1：初始化邻接顶点索引数组VNList[M]，将其元素均赋值为NULL。

2.2：初始化顶点序号循环变量i₂=0。

2.3：将序号为i₂的顶点邻接面片链表中所有面片的顶点均插入到VNList[i₂]中，如果顶点重复则放弃插入顶点。具体方法为：根据第一步的输出结果vertexNeighberList[M]可知序号为i₂的顶点的所有邻接面片链表为vertexNeighberList[i₂]，易知链表中所有面片的顶点都是序号为i₂的顶点的邻接顶点，因此，将这些面片的所有顶点都插入到VNList[i₂]中，即可得到序号为i₂的顶点的所有邻接顶点。每个顶点的插入方法为：遍历VNList[i₂]链表，如果该顶点已经存在，则转到2.4步；如果不存在则将该顶点作为链表最后一项插入到链表中。

2.4：循环变量i₂增加1。

2.5：若i₂小于M则转到2.3步，否则执行第三步。

第三步：根据vertexList[M][6]、faceList[N][3]计算各个顶点的单位法向量数，并将单位法向量数存到数组vertexNormalList[M][3]中，vertexNormalList[M][3]元素为浮点数，存储每个顶点法向量的3个分量。各个顶点的法向量可以作为矢量投影表面的法向量。各个顶点法向量数计算方法是：

3.1：初始化顶点序号循环变量i₃等于0。

3.2：根据序号为i₃的顶点的邻接面片链表vertexNeighberList[i₃]计算序号为i₃的顶点的各个邻接面片的单位法向量。并将所有面片单位法向量的平均向量作为序号为i₃的顶点的法向量，存入vertexNormalList[i₃]之中。面片的单位法向量计算方法为：设面片的三个顶点分别为v₀、v₁、v₂，则该面片的法向量为向量(vertexList[v₁][0]-vertexList[v₀][0],vertexList[v₁][1]-vertexList[v₀][1],vertexList[v₁][2]-vertexList[v₀][2])与向量(vertexList[v₂][0]-vertexList[v₀][0],vertexList[v₂][1]-vertexList[v₀][1],vertexList[v₂][2]-vertexList[v₀][2])的叉积，计算出法向量后，将该向量的各个分量分别除以该向量的长度即可得该面片的单位法向量。

3.3：循环变量i₃增加1。

3.4：若i₃小于M则转到3.2步，否则转到第四步。

第四步：根据vertexList[M][6]、faceList[N][3]，采用速度与方向的判断方法，检测流场中的临界点，并将临界点存入临界点链表，criticalPL为临界点链表的首指针，链表中元素为描述临界点信息的结构，该结构包括顶点号、下一个元素指针。

4.1：初始化临界点链表首指针criticalPL为NULL。

4.2：初始化顶点序号循环变量i₄=0。

4.3：根据序号为i₄的顶点速度大小判断该顶点是否为临界点，如果序号为i₄的顶点速度值等于0，该顶点即为临界点，转4.4；如果该顶点速度大于0，则直接执行4.5步。

4.4：将序号为i₄的顶点的顶点号插入到临界点链表最后一项的顶点号域，criticalPL向后移一个元素，转到4.5步；

4.5：循环变量i₄增加1。

4.6：若i₄小于M则转到4.3步，否则转到4.7步。

4.7：初始化面片序号循环变量i₅=0。

4.8：根据4.1步至4.5步的结果判断序号为i₅的面片的三个顶点中是否存在临界点，如果在4.1步至4.5步中已经检测到该三个顶点中存在临界点，则无需检测该面片，直接跳转到4.11步；否则执行4.9步。

4.9：采用基于速度方向判断的方法判断序号为i₅的面片是中否为临界点，具体判断方法为：如果该三角面片的三个顶点上的速度值在直角坐标系下按顺时针排列顺序与三角面片顶点的排列顺序相同且相邻向量之间对应方向的夹角小于180°，则判定三角面片内存在临界点。如果该面片内存在临界点，则跳转到4.10步；如果该面片内不存在临界点，则直接执行4.11步。

4.10：将该面片中速度最小顶点的顶点号插入到临界点链表的最后一项的顶点号域，criticalPL向后移一个元素。

4.11：循环变量i₅增加1。

4.12：若i₅小于N则转到4.8步，否则转到第五步。

第五步：根据第四步得到的临界点链表，计算每个临界点旋转因子，根据旋转因子判断临界点是否为涡核点，删除非涡核点的临界点，得到涡核点链表。

5.1：初始化临界点遍历指针Pr=criticalPL。

5.2：令采样半径r为曲面流场中面片平均边长的5倍。

5.2.1：建立采样点链表，并将其首指针SVList赋值为NULL。采样点链表存储了Pr所指向临界点的所有采样点，采样点链表元素由顶点号、下一个元素的指针组成。

5.2.2：采用采样点面片扩展搜索方法根据邻接顶点索引数组VNList[M]，搜索Pr所指向临界点在采样半径为r时的所有采样顶点，并将这些采样顶点存入采样点链表。采样点面片扩展搜索方法具体为：

5.2.2.1：设Pr所指向临界点为C，将临界点C作为采样点链表的第一项插入到采样点链表中。面片扩展搜索方法中一个重要的循环控制变量为当前面片扩展顶点，当前面片扩展顶点是指当前需要向外扩展的顶点，且需遍历采样点链表的所有顶点。设置当前面片扩展顶点为采样点链表的第一项，令Pr₂为指向采样点链表的第一项的指针。

5.2.2.2：设Pr₂指向的采样点链表项的顶点序号为i₆，根据邻接顶点索引数组VNList[i₆]，找到序号为i₆的顶点的所有邻接顶点。在这些顶点中，与临界点C的距离小于采样半径r的顶点即为临界点C采样顶点。因此，需将序号为i₆的顶点的所有邻接顶点中的临界点C采样顶点依次插入到采样点链表中。每个临界点C采样顶点的插入过程为：首先遍历采样点链表，如果采样点链表中无此临界点C采样顶点，插入到链表的最后一项之后，并转至5.2.2.3步；如果采样点链表中有此采样顶点，则放弃插入，直接执行5.2.2.3步。

5.2.2.3：Pr₂后移一个元素。

5.2.2.4：如果Pr₂为空，则采样点搜索完毕，完成面片扩展搜索，因为采样点链表中的第一项为临界点，不属于采样点，所以需将采样点链表第一项素删除，转至5.3步；如果Pr₂不为空，则跳转至5.2.2.2步。

5.3：根据采样点链表采用旋转因子计算方法计算临界点C的旋转因子，具体方法为：

5.3.1：根据5.2步计算的临界点C的采样点链表统计以下数据：采样点总数N_t，顺时针采样点数目N_s，逆时针采样点数目N_n。其中顺时针采样点和逆时针采样点分别是指：如果采样点的位置矢量（采样点位置与临界点C的位置的差矢量）到采样点矢量（采样点流场矢量在临界点C处的切平面上的投影，该切平面与临界点C的法向量垂直）的顺时针方向角度位于[π/2-ε,π/2+ε]范围内，则称该采样点为顺时针采样点，其中ε为角度判定的阈值，用以容忍矢量夹角的偏差，取值为π/8；反之，如果采样点位置矢量到采样点矢量的逆时针方向的夹角位于[π/2-ε,π/2+ε]范围内，则称该采样点为逆时针采样点。

5.3.2：计算顺时针采样点百分比P_s＝N_s/N_t，逆时针采样点百分比为P_n＝N_n/N_t。则旋转评定值为P_o＝max(P_s,P_n)，P_o为旋转因子的第一项。

5.3.3：根据临界点C的采样点链表，计算所有采样点的位置矢量和顶点矢量的夹角A_i，其中i取值为0,1,…N_t，分别表示N_t个采样点。则顶点矢量排列顺序和位置向量排列顺序的匹配程度值其中R_o为旋转因子的第二项。

5.4：判断临界点C是否符合涡核点的标准，也就是如果P_o＞1-ε_p且R_o＜ε_r，则临界点C为涡核点，其中ε_p和ε_r为旋转因子判定阈值，分别取值为0.25；否则临界点C为非涡核点，从临界点链表中删除该临界点。

5.5：Pr后移一个元素。

5.6：如果Pr为空，则临界点链表遍历完毕，非涡核点已经删除，此时的临界点链表即为涡核点链表，转至第六步；如果Pr不为空则跳转至5.2.2步。

第六步：检测所有涡核点的涡区范围：

6.1：初始化涡核点遍历指针Pr₃的初始值为涡核点链表首指针criticalPL，此时的临界点链表已将非涡核点删除，链表中的元素均为涡核点，链表首指针也就为criticalPL。

6.2：根据采样半径R（R为5.2步中r的2倍），利用5.2.2中的采样点面片扩展搜索方法，搜索Pr₃所指向涡核点的周围所有采样点，并将这些采样点作为涡区候选顶点。

6.3：采用5.3步的旋转因子计算方法计算所有涡区候选顶点的旋转因子第一项P_o。

6.4：如果涡区候选顶点的P_o＞E_p，其中E_p为涡区判断阈值，取值为0.5，则将该涡区候选顶点标记为涡区顶点。

6.5：Pr₃后移一个元素。

6.6：如果Pr₃为空，则涡核点链表遍历完毕，涡区检测完毕，转第七步；否则跳转至6.2步。

第七步结束。

与现有技术相比，采用本发明可以达到以下技术效果：

准确检测曲面涡特征的位置与范围，有效提高了涡特征检测的准确率，相比基于信息熵的检测方法采用旋转因子这一评价标准可以有效去除鞍点区域，准确率有了非常大的提高。

附图说明：

图1是本发明总体流程图。

图2是本发明第四步中临界点检测流程图。

图3是邻接面片索引数组vertexNeighberList[M]结构示意图。

图4是根据速度方向判断面片中是否存在临界点的判断原理示意图。

具体实施方式：

图1是本发明总体流程图。

1：由曲面顶点数据vertexList[M][6]、网格数据faceList[N][3]，得到各个顶点的邻接面片索引vertexNeighberList[M]。

2：根据vertexList[M][6]、faceList[N][3]得到各个顶点的邻接顶点索引数组VNList[M]。

3：根据vertexList[M][6]、faceList[N][3]计算各个顶点的单位法向量数，并将单位法向量数存到数组vertexNormalList[M][3]中。

4：根据vertexList[M][6]、faceList[N][3]，采用速度与方向的判断方法，检测流场中的临界点，并将临界点存入临界点链表。

5：根据临界点链表计算每个临界点旋转因子，根据旋转因子判断临界点是否为涡核点，删除非涡核点的临界点，得到涡核点链表。

6：检测所有涡核点的涡区范围。

7：结束。

图2是本发明第四步：临界点检测流程图。

4.1：初始化临界点索引链表首指针criticalPL为NULL。

4.2：初始化顶点索引变量i₄=0。

4.3：判断顶点i₄速度是否为0，如果是则转到4.4步；否则转到4.5步。

4.4：将顶点i₄插入到临界点链表中。

4.5：令i₄=i₄+1。

4.6：判断i₄是否小于总顶点数M，是则转到4.3步；否则转至4.7步。

4.7：初始化面片索引变量i₅=0。

4.8：判断面片i₅的三个顶点中是否存在4.2步至4.6步中检测到的临界点，如存在，则转到4.11步；否则转到4.9步。

4.9：根据速度方向判断面片i₅中是否存在临界点，如存在，则转到4.10步；否则转到4.11步。

4.10：将面片i₅的三个顶点中速度最小的一个顶点插入到临界点链表中。

4.11：令i₅=i₅+1。

4.12：判断i₅是否小于总面片数N，是则转到4.8步；否则转至4.13步。

4.13：结束。

图3是邻接面片索引数组vertexNeighberList[M]结构示意图。

该数组的元素为链表的首指针，每个指针指向相应顶点的邻接面片链表。包括面片号、链表中下一元素指针。

图4是根据速度方向判断面片中是否存在临界点的判断原理示意图。

设ΔABC三个顶点A、B、C上投影到三角面片的向量分别为如果在直角坐标系中相邻两个向量的夹角都小于180°的话，根据插值方法可知，三角面片内必然存在一个速度为0的数据点。同时，由涡旋的旋转特征可知三个向量在直角坐标系中存在特定的排列关系，也就是逆时针排列的顶点A、B、C上对应的向量在直角坐标系下必须为逆时针排列；顺时针排列的顶点A、B、C上对应的向量在直角坐标系下必须为顺时针排列。因此，我们得到判断插值离散数据场中三角面片内含有临界点的充分必要条件：顶点A、B、C上对应的向量在直角坐标系下按时针方向的排列与三角面片顶点的排列相同且相邻向量之间对应时针方向的夹角小于180°。