一种基于本征正交分解的速度场坏矢量识别和修正方法

摘要

本发明公开了一种基于本征正交分解的速度场坏矢量识别和修正方法，包括：A：对原始流场进行分区，得到大样本速度场；B、对大样本速度场进行本征正交分解（POD），判断不满足迭代收敛条件时，依据计算所得的最优重构阶数以及POD得到的模态和系数进行流场重构；C、计算原始流场与重构流场间的速度矢量绝对误差的模，判断所述速度矢量绝对误差的模大于预设的阈值时，确定对应的速度矢量为坏矢量，并剔除；D、在下次迭代之前，在已剔除坏矢量的位置进行基于POD的插值，修正被剔除矢量；之后返回B，继续对修正后的大样本速度场进行POD，直至满足迭代收敛条件，从而实现单一瞬时速度场的坏矢量识别和修正。

说明书

技术领域

本发明涉及流体力学速度测量技术领域，尤其涉及一种基于本征正交分解 的速度场坏矢量识别和修正方法。

背景技术

粒子图像测速（Particle Image Velocimetry，简称PIV）是一种现代激光测 速技术，主要运用于流场速度测量，通过追踪示踪粒子在流场中的运动来得到 流场速度场。PIV技术通过片光源或者体光源能实现二维或者三维速度场测量。 通常情况下，PIV实验测量得到的流场速度场会出现异常速度矢量，即常说的 坏矢量。流场速度场中出现坏矢量的原因很多，主要是由于粒子图像成像质量 不好，在PIV获取速度矢量的互相关计算时发生相关性分析失败所致。因此， 对实验测量得到的速度场进行坏矢量识别、并修正是保证实验数据正确性和精 度的一个非常重要的数据后处理环节。

现有速度场后处理技术中，对速度场坏矢量的识别和修正通常是分开实现 的，即先识别坏矢量，然后对坏矢量进行修正。在坏矢量识别方法中，一般通 过对比被检测点与其周围邻域内速度分布的统计关系来判定该点是否为坏矢 量。常见的识别方法包括：根据当地速度矢量均方差来进行偏差检验的识别方 法，还有利用当地速度矢量的中值或归一化中值来进行检验的方法等。其中， 所述归一化中值检测方法（normalized median test）是一种普遍被采用的计算效 率较高的方法。但是，由于归一化中值检测方法基于局部相邻速度矢量进行检 测，因此存在两个突出的问题：第一、检验方法是一种针对矢量场的纯数学处 理方法，并没有考虑到矢量场背后的流动机理，使得速度场坏矢量剔除和后续 的修正会引入非物理的偏差；第二、局部检验的方法往往在速度场出现较大坏 矢量区域时会失效，特别是当算法没有引入迭代识别机制时问题尤为突出。

另一方面，对于坏矢量识别后的修正，常见的修补方法有双线性插值、样 条插值、克利金插值等。这些方法虽然便于应用，但针对流场插值的正确性并 不理想，只是一种单纯的数学算法，不能反映真正流场的信息。当流场有强剪 切区域存在，即速度的梯度变化剧烈时，或者当坏矢量连通区域较大时，这些 插值方法无法有效地还原流场的真实信息。为了弥补这个缺陷，基于本征正交 分解（POD）的Gappy POD方法被人们应用于修正速度场。但是Gappy POD 方法在对流场进行修正时，需要事先知道坏矢量的位置，因此，它不具备坏矢 量识别功能，且目前常见的Gappy POD方法在迭代过程中对模态阶数的选择具 有一定难度，影响计算效率。

此外，也有基于流场大样本数据的POD速度场后处理技术。此类方法通过 全流场本征正交分解获得的模态信息，对坏矢量进行识别并完成修正，但它只 适用于处理由大量瞬时速度场组成的大样本实验数据，不适用于单一瞬时速度 场的后处理，因此该方法不能被应用于PIV数据处理中多层互相关分析中间步 骤的坏矢量修正，只能在得到全部PIV速度场后进行后处理。另外，现有此类 方法在低阶模态重构参考速度场时，采用了通过能谱曲线交点来动态选取模态 阶数的方法。但是这种方法很不稳定，会出现所选模态数不单调增加的现象， 尤其是在迭代后期、当能谱曲线逼近真实流场能谱曲线时，很难找到清晰的交 点，这为模态阶数的选取带来很大不确定度，往往造成算法失效。

发明内容

为解决现有存在的技术问题，本发明实施例提供一种基于本征正交分解的 速度场坏矢量识别和修正方法。

本发明提供的一种基于本征正交分解的速度场坏矢量识别和修正方法，该 方法包括：

A：对原始流场进行分区，得到大样本速度场；

B：对大样本速度场进行本征正交分解POD，判断不满足迭代收敛条件时， 依据计算所得的最优重构阶数以及POD得到的模态和系数进行流场重构；

C：计算原始流场与重构流场间的速度矢量绝对误差的模，判断所述速度 矢量绝对误差的模大于预设的阈值时，确定对应的速度矢量为坏矢量，并剔除；

D：在下次迭代之前，在已剔除坏矢量的位置进行基于POD的插值，修正 被剔除矢量；之后返回B，继续对修正后的大样本速度场进行POD，直至满足 迭代收敛条件。

上述方案中，所述原始流场为二维时，设流场在（x，y）方向上的网格节 点数分别为N,M，子域的网格节点数分别为n_B,m_B，并设置参数 R_B=(n_B×m_B)(N×M)，那么，

所述对原始流场分区满足的条件为：所述R_B的取值在0.2到0.5之间。

上述方案中，所述分区后得到的大样本速度场，即流场子域记为V_s，所述 s代表第s个子域，所述s作为索引的最大值为（N-n_B+1）*（M-m_B+1），最小 值为1。

上述方案中，所述对大样本速度场进行POD后，得到模态Φ_n和系数a_n， 则V_s可表示为其中，所述n_p为POD分解得到的总模态数，n表示 模态序列索引，其取值范围为（1，n_p），a_n,s表示系数a_n中针对V_s的系数。

上述方案中，所述迭代收敛条件为：所述e⁴代表前四阶的 相对能量，k代表循环次数。

上述方案中，所述计算最优重构阶数的方法，为：

以分区后得到的第一个子域V₁作为基准，采用高斯滤波或中值滤波进行平 滑处理得到平滑子域然后，用不同的前n阶模态重构所述平滑子域；计算 所述重构子域与平滑子域之间的相关系数R₁，选择使R₁最大的重构阶数作为 最优重构阶数n_opt；其中，所述R₁是重构阶数的函数。

上述方案中，所述流场重构的方法为：

用前n_opt阶模态重构所述流场V_s，重构后的流场记为s代表第s个子域， 重构公式为：${\tilde{V}}_s=\underset{n=1}{\overset{n_{\mathrm{opt}}}{\Sigma }}{a}_{n,s}{\Phi }_n.$

上述方案中，所述判断所述速度矢量绝对误差的模大于预设的阈值时，确 定对应的速度矢量为坏矢量，为：

设置针对第s个子域的残差阈值为残差阈值如果各网格 点上速度矢量绝对误差的模大于则该矢量被认为是坏矢量；

其中，所述c代表所取阈值的强度，取值为c=3；所述表示 每一子域误差的平均值；所述表示每一子域的标准差。

上述方案中，所述在已剔除坏矢量的位置进行基于POD的插值，修正被剔 除矢量的方法，包括：

对分区后的流场进行POD；

计算最优重构阶数；

依据所述最优重构阶数以及POD得到的模态和系数进行流场重构；

计算流场中各子域的截断残差；

利用重构流场和坏矢量周围相邻点的截断残差中值来修正被剔除的坏矢 量。

上述方案中，所述流场为二维时，所述利用重构流场和坏矢量周围相邻点 的截断残差中值来修正被剔除的坏矢量，为：

按公式计算坏矢量处的修正值；

其中，所述和分别代表坐标点（i，j）处坏矢量的修正后速度矢量和 重构速度矢量；函数median(r_s,ij)代表坐标点（i，j）周围3×3区域内的8个点相 邻点上截断残差的r_s,ij的中值。

本发明提供的基于本征正交分解的速度场坏矢量识别和修正方法，包括： A：对原始流场进行分区，得到大样本速度场；B、对大样本速度场进行本征正 交分解（POD），判断不满足迭代收敛条件时，依据计算所得的最优重构阶数以 及POD得到的模态和系数进行流场重构；C、计算原始流场与重构流场间的速 度矢量绝对误差的模，判断所述速度矢量绝对误差的模大于预设的阈值时，确 定对应的速度矢量为坏矢量，并剔除；D、在下次迭代之前，在已剔除坏矢量 的位置进行基于POD的插值，修正被剔除矢量；之后返回B，继续对修正后的 大样本速度场进行POD，直至满足迭代收敛条件，从而实现单一瞬时速度场的 坏矢量识别和修正。

可见，相比较传统的基于局部流场信息的速度场坏矢量识别和修正方法， 本发明将坏矢量的识别和修正很好的统一起来，通过对全流场流动结构的迭代 逼近来实现速度场优化，使得修正算法具有更高的精度和更强的物理含义。另 一方面，相对于现有采用POD的全流场修正方法，本发明能克服POD方法只 适用于大样本数据分析、而不适用于PIV互相关分析的中间步骤数据修正这一 缺陷，采用流场分区技术并结合更稳定的流场模态选取和坏矢量修正方法来实 现单一瞬时速度场的POD修正算法，使得速度场修正技术能应用于PIV互相 关分析的中间步骤。

附图说明

在附图（其不一定是按比例绘制的）中，相似的附图标记可在不同的视图 中描述相似的部件。具有不同字母后缀的相似附图标记可表示相似部件的不同 示例。附图以示例而非限制的方式大体示出了本文中所讨论的各个实施例。

图1为本发明实施例所述的基于本征正交分解的速度场坏矢量识别和修正 方法实现流程示意图；

图2为本发明一具体实施例所述的基于本征正交分解的速度场坏矢量识别 和修正方法实现流程示意图；

图3为本发明具体实施例所述瞬时速度场的划分示意图；

图4为本发明具体实施例所述最优重构阶数的选择示意图；

图5为本发明具体实施例所述速度场坏矢量修正方法流程示意图；

图6为本发明一具体应用场景中的初始速度场云图；

图7为本发明一具体应用场景中的坏矢量修正后的速度场云图。

具体实施方式

本发明的实施例中：对原始流场分区后得到的大样本速度场进行本征正交 分解（POD），判断不满足迭代收敛条件时，依据计算所得的最优重构阶数以及 POD得到的模态和系数进行流场重构；计算原始流场与重构流场间的速度矢量 绝对误差的模，判断速度矢量绝对误差的模大于预设的阈值时，确定对应的速 度矢量为坏矢量，并剔除；在下次迭代之前，在已剔除坏矢量的位置进行基于 POD的插值，修正被剔除矢量；之后继续对修正后的大样本速度场进行POD， 直至满足迭代收敛条件。

下面结合附图及具体实施例对本发明作进一步详细说明。

图1为本发明实施例所述的基于本征正交分解的速度场坏矢量识别和修正 方法实现流程示意图，如图1所示，包括：

步骤101：对原始流场进行分区，得到大样本速度场；

步骤102：对大样本速度场进行本征正交分解（POD），判断不满足迭代收 敛条件时，依据计算所得的最优重构阶数以及POD得到的模态和系数进行流场 重构；

步骤103：计算原始流场与重构流场间的速度矢量绝对误差的模，判断速 度矢量绝对误差的模大于预设的阈值时，确定对应的速度矢量为坏矢量，并剔 除；

步骤104：在下次迭代之前，在已剔除坏矢量的位置进行基于POD的插值， 修正被剔除矢量，之后返回步骤102，继续对修正后的大样本速度场进行POD， 直至满足迭代收敛条件。

图2为本发明一具体实施例所述的基于本征正交分解的速度场坏矢量识别 和修正方法实现流程示意图，如图2所示，包括：

步骤201：对原始流场进行分区，得到大样本速度场；

具体为：将待处理的已知瞬时速度场（即原始流程，二维、三维均可，本 实施例以二维流场为例）按图3所示的遍历扫描顺序划分成统一尺寸的流场子 域，各子域可以有重叠区域，其在x,y两个方向上的偏移为常数个网格节点，最 简单的可取1个网格节点。

这里，假设速度场在（x，y）方向上的网格节点数分别为N,M，子域的网 格节点数分别为n_B,m_B。设置参数R_B=(n_B×m_B)(N×M)，R_B越大说明子域的尺寸 越大，相应子域的个数越少，反之亦然。通过实验测试证明，所述R_B的取值应 在0.2到0.5之间。将分区得到的速度场子域依照分区扫描的顺序组合成大样本 速度场，记为V_s，其中s代表第s个子域，所述s作为索引的最大值为（N-n_B+1） *（M-m_B+1）。

步骤202：对分区后的流场进行POD，并判断前四阶的能量是否满足迭代 收敛条件，确定满足时，结束当前流程，停止迭代；否则执行步骤203；

具体的：由于坏矢量对流场的POD分解有影响，因此算法需要通过迭代来 逼近最优速度场。每个迭代步中将大样本速度场子域数据进行本征正交分解（第 一次迭代时，对步骤201中所得的大样本速度场子域数据进行POD，后续迭代 过程对前次迭代所得的修正后的大样本速度场子域数据进行POD），得到模态 Φ_n和系数a_n，则大样本速度场，即流场子域V_s可表示为其中n_p为 POD分解得到的总模态数，n表示模态序列索引，其取值范围为（1，n_p），a_n,s表示系数a_n中针对V_s的系数。

这里，迭代收敛条件可取小于某一阈值，通常取此处e⁴代表前四阶的相对能量，k代表循环次数，当k=1时可不用进行本步骤所 述收敛条件的判断，如果满足收敛条件，则跳出循环；否则进行步骤203。

步骤203：计算最优重构阶数；

具体为：最优重构阶数是一个关键参数，它决定了坏矢量剔除和填补的效 率和精度。可选择任一子域作为基准，为了方便我们选择分区后得到的第一个 子域V₁作为基准。将这个子域进行较强地平滑（高斯滤波或中值滤波等均可）， 以平滑坏矢量的影响，获得平滑子域然后，用不同的前n阶模态重构该子 域，重构的过程就是对平滑子域不断逼近的过程。计算重构子域与平滑子域 之间的相关系数R₁，该相关系数是重构阶数n的函数，如图4所示。一般随着n 的增加，R₁先变大后减小。当相关系数最大时，说明重构子域既足够光滑又能 最大限度地逼近原流场。因此，选择使R₁最大的重构阶数做为最优重构阶数n_opt。

步骤204：依据所述最优重构阶数以及POD得到的模态和系数进行流场重 构；

具体为：用前n_opt阶模态重构流场V_s，重构后的流场记为s代表第s个 子域。重构公式为：Φ_n,a_n,s为在步骤102中本征正交分解获得的模 态和系数。

步骤205：计算原始流场与重构流场间的速度矢量绝对误差的模；

具体为：计算原始流场V_s与重构流场之间各网格点上速度矢量绝对误差 的模，计为假设子域在两个坐标方向上的网格节点数分别为n_B,m_B，且网格 节点的索引分别为i，j，其取值范围分别为（1，n_B）和（1，m_B），则坐标（i， j）上的速度矢量绝对误差的模为：该点速度矢量V_s,ij和差的二范数，即 同时，计算出每一子域误差的平均值和标准差 ${\sigma }_s=\mathrm{std}\left(r_s^{*}\right).$

步骤206：判断所述速度矢量绝对误差的模大于预设的阈值时，确定对应 的速度矢量为坏矢量，并剔除；

具体为：设置针对第s个子域的残差阈值按照该阈值剔除坏矢量。残 差阈值其中c代表所取阈值的强度，即：几倍标准差，通常可取 c=3。如果各网格点上速度矢量绝对误差的模大于则该矢量被认为是坏矢 量，进行剔除处理。

在不同的迭代步中，参数c为常数，而阈值为动态改变的，能自适应流 场坏矢量的变化。通过大量实例证明，c=3能适应不同的流动结构，是一个较 稳定的通用参数。

步骤207：在下次迭代之前，在已剔除坏矢量的位置进行基于POD的插值， 修正被剔除矢量，之后返回步骤202，进行下一次迭代；

具体为：在已剔除坏点的位置进行基于POD的插值，修正被剔除的数据。 实验证明该方法的精度和效率都比较高。此步骤也可以通过迭代来提高精度， 通常只需要迭代两步即能实现较高精度的数据修正。该方法按图5流程的具体 实施步骤如下：

步骤2071：对分区后的流场进行POD；

具体的，该算法同样是建立在通过POD最优重构迭代逼近真实流场的基础 上的。在从步骤2071到步骤2075的第一次迭代循环中，被处理的速度场子域V_s对应的模态Φ_n和系数a_n为步骤202中的计算结果（此时的速度场为已剔除坏矢 量的速度场，该速度场还未进行过修正）；在迭代过程中间步对修正后的速度场 子域进行处理时，需要对修正后的速度场进行本征正交分解，得到相应模态Φ_n和系数a_n，两种情况下速度坏矢量的位置均为已知的。

步骤2072：计算最优重构阶数；

具体的，在步骤2071中获得模态Φ_n和系数a_n后计算最优重构阶数n_opt，其 计算方法与速度场坏矢量识别中步骤203一致。

步骤2073：依据所述最优重构阶数以及POD得到的模态和系数进行流场 重构；

具体的，采用前n_opt阶模态重构得到重构流场其计算方法与速度场坏矢 量识别中步骤204完全一致。

步骤2074：计算流场中各子域的截断残差；

按公式计算各子域初始流场和POD重构流场之间的截断残差r_s。

步骤2075：利用重构流场和坏矢量周围相邻点的截断残差中值来修正被 剔除的坏矢量；

具体的，按公式计算坏矢量处的修正值，和分 别代表坐标点（i，j）处坏矢量的修正后速度矢量和重构速度矢量；函数median(r_s,ij) 代表坐标点（i，j）周围3×3区域内的8个点相邻点上截断残差的r_s,ij的中值。

优选的，如果采用迭代算法，该方法还包括：返回步骤2071，对修正后的 流场再次进行POD，以修正被剔除的坏矢量。

这里，通常迭代循环的次数可以设置为2次，完成整个步骤207后进入下 一循环，重复步骤202。

下面结合一具体应用场景对本发明所述的方法进行描述。

速度场坏矢量的识别和修正的实施是针对如图6所示的具有120×120个网格 节点的已知瞬时速度场进行的。该速度场包含14400个速度矢量，记为速度场 A。

步骤一：将速度场A按步骤201划分子域，子域的网格点数为70×70，子 域与子域在x,y两个方向上的偏移均为1。根据步骤201中的设置 R_B=(n_B×m_B)(N×M)，可知R_B=0.34，落在[0.2,0.5]的范围之内。将划分得到的 速度子域记为V_s，s代表第s个子域，所述s的最大值是（N-n_B+1）*（M-m_B+1） =51*51=2601，因此，s的取值范围为（1,2601）。

步骤二：通过迭代实现流场的坏矢量识别和修正，每个迭代步中将初始子 域V_s或上一个迭代得到的修正后速度子域进行本征正交分解，得到相应的模态 Φ_n和系数a_n。迭代的收敛条件为e⁴代表前四阶的相对能量，k 代表循环次数，当k=1时不用判断。如果满足收敛条件则跳出循环，不满足则 继续步骤三。

步骤三：计算最优重构阶数。选择第一个速度子域V₁做为基准，对其进行 过滤窗口为5×5个节点的中值滤波，得到光滑的速度场用不同的阶数n重构 该速度场并计算重构速度场与平滑速度场之间的相关系数R₁，该相关系数是重 构阶数n的函数。图4给出了本实施例计算的结果。随着n的增加，R₁先变大后 变小。这说明随着模态的增加，重构的流场开始逼近真实解，后来逐渐受坏点 的影响偏离真实解。选择使R₁最大的重构阶数作为最优重构阶数n_opt。本实例中 第一次迭代选择n_opt=16，第二次迭代是26，第三次迭代是54，图4中的a、b 和c代表这三次迭代。

步骤四：对子域V_s处理时，用前n_opt阶模态重构流场，重构后的流场记为重构公式为：所述a_n,s,Φ_n已在步骤二的本征正交分解中得到。

步骤五：根据步骤205计算速度子域V_s与重构流场之间的绝对误差的模 同时计算出对应的平均值和标准差

步骤六：设置针对第s个子域的残差阈值残差阈值按照 该阈值剔除坏矢量，其中参数c=3。

步骤七：根据步骤207进行两次迭代算法，在剔除坏点的位置进行基于本 征正交分解的插值，修正缺失的数据。

步骤八：迭代次数k加1，重复步骤二，直到满足迭代收敛条件结束迭代过 程。

本实施例方法在本征正交分解的基础上识别和修正坏矢量，能够在保证原 始正确数据不被影响的情况下，通过POD低阶重构流场对比的方法识别流场中 的坏矢量，同时通过对重构流场叠加坏矢量周围网格点上速度截断残差中值的 方法来实现坏矢量的修正。本方法实施后，如图6所示的原始带有坏点的流场 被修正成如图7所示的较光滑流场。

本发明通过本征正交分解将原始待处理的流场分解成各阶模态的线性组 合，低阶模态包含流动的大尺度结构，高阶模态包含流场噪音的影响。通过本 发明中的模态重构方法可以最大限度地迭代逼近原始流场。计算原始待处理的 流场和重构流场之间的绝对误差，当绝对误差大于事先设置的阈值时，该速度 矢量被判断为坏点。在进入下一次迭代之前，已被判断出的坏矢量需要进行插 值运算。该插值方法同样基于本征正交分解，被填充的值是最优重构值和误差 修正值的和，实验证明该插值方法具有很高的精度。经过有限次迭代之后，流 场的模态收敛且不再变化，流场趋于稳定，几乎不再含有坏点。

以上所述，仅为本发明的较佳实施例而已，并非用于限定本发明的保护范 围。