基于本征正交分解(POD)的PIV数据坏点剔除方法

基于本征正交分解（POD）的PIV数据坏点辨识方法

摘要：

1.引言：

粒子图像测速技术是一种基于图像互相关的激光测速技术。在进行PIV实验的过程的由于示踪粒子浓度不均匀、激光强度分布不均、粒子成像质量差等原因容易造成互相关峰值不确定1，出现所谓的“坏点”。最近几年发展起来的3D3CPIV2, 3，由于在重构空间粒子场的时候会不可避免的引入虚假粒子，进一步增加了流场中了坏矢量。因此，坏矢量的剔除和重建是PIV数据后处理的一个重要内容。一般通过对比该点误差与周围3*3或5*5邻域内平均误差来确定该点是不是坏点。目前比较通用并且效果很好的是Westerweel和Scarano 提出的归一化中值检测方法4（normalized median test）。考虑一个位移矢量U0和其周围3⨯3的相邻矢量{U1,U2,U3,U4,U5,U6,U7,U8}。用Um表示{U1,U2,U3,U4,U5,U6,U7,U8}的中值，用ri表示残差，定义为ri=Ui-Um

（i=1,2,„8）。U0的残差用ri的中值rm归一化后得到下面的公式：

r0*=0-Um

rm+ε

其中ε是与流场的平均噪音有关的量，一般设为0.1-0.2 pixel。对于平面PIV，r0*的阈值设为24，大于该阈值的点被认为是坏点。其他坏点识别的方法与该方法原理相同，基本上都是基于相邻点的统计特性判断该点的性质。这种方法有缺陷，一是只能用到少量相邻点的值，忽略了流场的全场特性；二是局部相邻点的速度矢量很有可能也是坏点，导致判断的结果不可靠；三是r0*的阈值需要根据流中坏点的多少确定，在实际运用中并不能精确的知道坏点所占的比例。因此，本文作者提出了基于POD的坏点检测方法。该方法对周期性或类周期性流场PIV数据的后处理有很好的效果。

本文的结构安排如下：首先在第二节简单的介绍了本征正交分解的原理，推导了POD与流场湍动能的关系；在第三节对流场的误差进行了数值模拟，讨论了误差对POD分解的影响以及详细的介绍了基于POD的坏点剔除方法；一个真实PIV计算出的流场运用该方法进行了坏点剔除并和归一化中值检测方法进行了比较，这在第四节给出；第五节是结论。

2.POD原理

POD是一种从统计意义上提取流场中主要流动结构的方法，可以实现对复杂非线性系统的线性降维处理。这种方法最早由Lumley（1967年）引入湍流的研究，用于辨识大尺度拟序结构，后来Berkooz5等人对POD方法给出了系统的介绍。POD方法就是要找出和原流场最相似的流动模态5，对于复杂流场这样的流动模态不止一个，于是可以假设流场是一系列POD基模态的线性组合，即：

其中U(X,t)是不同时刻的速度场，Φk(X)代表POD的第k阶模态，ak(t)代表与第k

阶模态相关的时间系数。POD的模态满足正交关系：

⎧0, i≠jΦi,Φj=δi,j=⎨⎩1, i=j

从数学上可以从下面的公式解出Φk(X)： k=1U(X,t)≈∑ak(t)Φk(X)K

CΦk=λkΦk （1）

其中C是U(X,ti)和U(X,tj)的两点互相关系数，即C=(X,ti),U(X,tj)。这里我们着重讨论一下特征值λk的取值。对于两点互相关矩阵C中的元素Cij有计算公式：

Cij=(U(X,t),U(X,t))=K

一般都是用速度的脉动场进行POD分解的，式中U(X,t)代表的是速度脉动。当拥有足够多的不同时刻流场的数据时，不论是对于classical POD还是snapshot POD，Kii都代表流场的湍动能，而且不会有太大的变化，即在Cij的公式中分母基本都相同，这也正是可以直接用Kij替换Cij带入公式8计算的原因。因此，当对Φk归一化后，特征值λk就是一个与该阶模态的湍动能有关的量。λk越大，该阶模态所含的湍动能越多，也就是说POD可以按能量提取流场的主要特征。详细的介绍可以参考文献5。

3.误差模拟分析

3.1 模拟方法

评价PIV计算结果好坏的一个重要指标就是坏点所占的百分比Q。一般来说，在PIV中存在两类误差。一类是以单个点形式出现的错误矢量，这类误差大小和分布都带有明显的随机性，主要是由于查询窗内互相关信噪比太低峰值不明确造成的；另一类误差成片形式出现，通常有好几个误差向量集中在同一个区域，这种误差很有可能是由于图像在这一区域粒子太少或质量太差造成的。设误差向量的连通区域的大小为nerr，当nerr=1时表示单个误差矢量，nerr=2表示两个误差矢量连在一起，以此类推。为了得到nerr的分布，对一个PIV计算出来的真实的流场用归一化中值检测方法检识别坏点，然后统计坏点连通区域的分布。得到图 1中的统计结果。从直方图中可以看出单个误差矢量出现的概率最高，连通区域越大出现的概率越小。对数据进行高斯拟合得到光滑的拟合曲线：

f(nerr)=0.3439exp(-2nerr3.8152)(nerr=1,2,3...)

f(nerr)代表该类型误差占所有向量的比例。误差的大小在查询窗内随机分布。当流场误差总的个数确定后，按该概率密度函计算不同类型的误差并添加到基本流场之上，完成误差场的模拟。

图 1 误差向量连通区域nerr的概率密度分布

图 2 Q=3%时模拟结果

基本流场取自一套平滑过后的圆柱扰流的平面PIV数据。该数据在北京航空航天大学低速水槽完成，自由来流速度为35mm/s。圆柱水平放置在水槽中间直径为10mm。雷诺数为250。在流场中撒播直径为5μm ，密度为1.05gmm-3的空心玻璃微珠作为示踪粒子，用一台2w的激光器照亮测量区域。相机的分辨率为640*480，采样频率为

80Hz。在进行互相关计算时采用了窗口变形算法，查询窗最终大小为16*16,50%重叠区。 本文模拟了8中不同的工况，分别是Q=0.5%,1%,2%,3%,4%,5%,7%,10%。对其中Q=3%的情况进行了详细的分析。图 2给出了Q=3%时模拟的结果。

3.2 误差向量对POD分解的影响

POD分解按流场的湍动能提取拟序结构，而误差向量直接影响流场湍动能。直观上判断流场中坏点越多，流场的湍动能就会越大，而且这种湍动能的增加并不是流动结构的变化引起的。设Q代表误差向量的个数占流场全部向量个数的百分比，rQ代表引入误差后流场湍动能的值与真实流场湍动能值的比值，当Q=0时rQ=1。图 3给出了rQ随Q的变化规律。从图中可以看出，rQ明显随Q呈线性变化，增加1%的误差向量流场的湍动能增加23.27%。但是值得注意的是这种湍动能是由于随机误差造成的，并不包含流动结构。

图 3 流场湍动能随Q的变化规律

rQ=23.27Q+1

接下来我们考虑误差向量对POD分解的影响。在前面一节我们已经知道基于两点互相关的POD按湍动能提取流场的流动特征。那么，由于引入了误差导致湍动能增加会影响POD分解的模态吗？怎样影响？为此，我们选择包含0%，0.5%，1%，3%，5%，10%的误差的流场进行本征正交分解，根据拟合的公式可知流场的湍动能分别为：100%，110.75%，121.5%，164.5%，207.5%，315%。图 4给出了分解后的能量谱。由于流场整体能量的增加，各阶模态所占的相对能量百分比下降。设各阶模态的绝对能量为e

i，

E0为原始没有误差向量的流场所含总湍动能的绝对值。那么，原始场各阶模态的相对能量为：

λi=eie=i E0ei

当引入误差矢量后，流场的湍动能会按照线性规律增加。增加的比例用rQ表示。用λi'表示该流场模态的相对能量，于是得到：

ei'ei'ei' λ='=='EeirQE0'i

为了和原始模态进行比较，在上面的公式两端都乘以rQ，进行归一化处理。得到的结果如图所示。从图中可以看出，归一化后的各阶模态的相对能量曲线在低阶时很好的重合在一起。虽然流场的湍动能在不断增加，但是主要流动结构（低阶模态）的能量并没有发生变化。这说明误差向量并没有影响低阶模态的提取，用较少的低阶模态同样能重构出和原流场非常近似的流场。但是，从图中我们发现，随着误差向量的越来越多，对流场的影响逐渐深入到低阶高能模态，而且模态的阶数越来越多。这与流场的相关性降低有关。高阶模态反映了流场的随机误差，虽然随机误差场占有很大一部分湍动能，但这些湍动能并不包含流动结构并且每阶模态杂乱无章。低阶模态包含流场的大尺度流动结构，这些模态所含的能量越低，越容易受随机误差的影响。但是值得注意的是，随机误差在这些流动结构中所含的能量是十分有限的，大部分误差的能量仍然集中在高阶模态，并且各阶模态之间能量差别很小。正是这样的原因导致了POD模态的急剧增加，达到上千阶。对于含误差3%的PIV数据，如果选择前15阶和100阶模态进行重构得到的流场如图 5所示。从图中可以看出，前15阶得到的流场受随机误差的影响要小很多，重构阶数越多流场受到的扰动越大。而且随机误差会传播到流场中的每一个点，这种误差均匀分布于全流场，通过对比很容易判断出误差的分布。 图 4 不同Q下流场POD分解的能量谱

图 5 Q=3%流场与POD重构对比

通过上面的讨论得知，随机误差会快速的增加流场的湍动能，但是在一般PIV误差量级下，其能量在POD分解时还不足以影响流场的大尺度结构，即低阶高能模态。虽然各个点的误差会传播到整个流场，导致POD的模态不光滑，但并不影响对流场的近似重构。

3.3 基于POD的坏点剔除

通过上面的分析我们知道，运用较少的低阶高能模态重构可以引入较少的随机误差。对比重构的流场和原始流场，在存在坏点的地方速度矢量的差别很大，而不存在坏点的地方绝对误差与周围的绝对误差相近。由于重构采用了固定的阶数，不可避免的产生了阶段误差，但是从局部3*3 或 5*5 的区域来看，这种重构的阶段误差应该是非常相近的，这为我们判断坏点提供了方便。与Westerweel和Scarano4相同，考虑速度矢量V0和其周围3*3区域的相邻矢量{V1,V2,V3,V4,V5,V6,V7,V8}，与之相对应的前m阶重构的速度矢量用V0m和{V1m,V2m,V3m,V4m,V5m,V6m,V7m,V8m}。在平面PIV中每个速度矢量有两个分量u,v

。重构的绝对误差

riPOD=i-Vim=可以看出rPOD为两矢量差的长度。该误差与重构阶数、流场特性有关。一般而言在强剪切区或湍流度比较大的区域该误差较大，这与POD分解的特性有关，这些区域需要更多的重构阶数才能达到一定的精度。另外，该误差还与该点速度的大小有关。由此可见用绝对误差rPOD作为误差判断的标准并不能适应各种流动工况，一种简单的解决方法就是用周围点的绝对误差对其归一化。具体公式为：

r*POD

0r0POD =PODmin(ri)+ε

式中r0*POD代表V0点误差归一化后的残差，分母上为周围8个点的绝点误差的最小值。和Westerweel

和

Scarano文献不同，之所以选择最小值，是加速收敛和对成片坏点有更好的剔除效果。ε为一常数，反映流场中整个误差的平均水平。可以预见在光滑的区域r0*POD应该是趋近于1的，而在存在坏点的地方残差应该很大。需要说明的是，之所以分母上选择中值而非平均，是因为中值滤波可以减少坏点的干扰，算法上更稳定。这样可以通过迭代剔除坏点不断的逼近真实流场。算法的实现过程如下：

a) 准备第k次迭代

b) 对前一次迭代得到的流场进行POD分解，若k等于1，则对原始PIV数据进行POD

分解。这里要求原流场拥有足够多的帧数和周期，在时间上能够很好的解析整个流动过程。本文选择snapshot POD方法。

2c) 选择前4阶模态的相对能量作为收敛条件的判断。设ek代表第k次迭代后前4阶

模态的能量。随着坏点的不断剔除，低阶模态的能量会越来越大，最终趋于稳定。因此可以认为当22(ek-ek-1)2ek

g的取值范围为1~2%。

d) 计算重构阶数mk并用前mk阶模态对流场进行重构。当k=1时，mk=4。当k>=2

时，取原始数据能量谱和第k-1次迭代完成并插值后数据能量谱的交点对应的阶

数。如图XX所示，该点位置在迭代过程基本不变。

e) 按照公式计算残差r0*POD。式中εk=mean(riPOD)，即该次迭代所有误差的平均值。

并计算每一个时刻的残差的均值mean(r0*POD)和脉动std(r0*POD)。

f) 设置残差的阈值rth=mean(r0*POD)+3*std(r0*POD)，若残差大于该阈值则认为是坏

点被剔除。这里3是较为合理的，若太小则会剔除正确的矢量，若太大则要迭代

更多的次数。

g) 对已经剔除的点进行插值。插值点附近依然可能存在坏点，故插值的方法尽量对

坏点不敏感。线性插值、样条插值等一般插值方法都是基于周围数据的统计信息

进行插值计算的，很容易受到噪音的影响。Gappy POD是一种基于全流场信息的插值方法，在这里非常适用。但由于插值并不是本文的重点，这里并不引入该方法。本文选择了一种基于POD重构流场的插值方法。考虑第k次迭代并剔除坏点后的

速度矢量V0k和其周围3*3区域的速度{V1k,V2k,V3k,V4k,V5k,V6k,V7k,V8k}，与之相

mkmkmkmkmkmkmkkk对应的前mk阶重构的速度矢量用V0mV1mk,V2k,V3k,V4k,V5k,V6k,V7k,V8k}表k和{

mkmk示。用rkmk表示误差，即rik。 =Vik-Vik

若V0k被当做坏点剔除，则可用下面的公式计算其插值：

mkkV0'k=V0m+median(rkik)

h) 插值后的流场作为新的原始流场，进入第a步。

3.4 数值评估

按照前面介绍的误差模拟方法，对基本流场添加5%的误差向量并运用POD方法识别坏点。按照上面所描述的方法，设置收敛条件g=2%，6次迭代后满足收敛条件。由于基本流场是一个周期性很好的流动，前两阶模态的能量占流场湍动能的90%左右。因此计算出的重构阶数mk较小，大概为16，并且随迭代次数的变化很小，如图XX中国黑色圆圈所示。红色线代表真实流场的能谱，从图XX可以看出，随着迭代次数的增加，流场中的坏点被剔除，流场的能谱逐渐逼近真实解。同时由随机误差造成的高阶低能模态越来越少，最终POD模态的阶数也会逼近真实解。从中也可以看出，前几阶模态在真实流场中的能量应该是最大的。图XX给出了前四阶能量的收敛曲线和对应坏点比例曲线。在迭代剔除坏点的过程中，坏点会逐渐被判断出来，对应的前四阶能量会逐渐变大，直到达到收敛条件。最终整个流场识别出5.2%的坏点，非常接近预设的5%。

通过上面的分析，我们从整体上判断出流场的收敛性，但是并不清楚流场坏点识别和

填充的细节。从流场取出某一帧的数据，用红色标示出识别出的坏点，得到图XX。基本上所有坏点都被识别。按照POD插值后得到图XX的结果。图中颜色标示速度矢量的大小。通过计算插值结果与原数据的相对误差，给出了相对误差的概率密度分布（PDF），如图XX。相对误差落在[0,2%]范围内的概率为85%，落在[0,3%]的概率为91%。和一般的基于周围点的数学插值方法不同，本文所用的基于POD的插值方法利用了流场所有的时空信息，能够有效的还原流动结构，对成片的坏点和强剪切区域的插值有很好的效果。为了说明这一点，在真实流场中随机剔除速度矢量然后用两种方法插值，比较相对误差的大小。得到的结果如图XX所示。横坐标表示插值区域的大小，例如5，则代表5*5的插值区域。纵坐标是统计了上万个点后相对误差的平均值。可以清楚的看到，随着区域的增大，基于POD的插值方法的相对误差明显低于线性插值。图中仍然取前16阶进行插值。需要说明的是，之所以进行了统计分析，是因为插值的精度与插值的位置有关，速度梯度变化大的地方精度会变低。同时本文所用的插值方法与重构阶数的选择有关，尽量选择没有被误差影响的模态，如步骤d所示。只要阶数选择合理，基于POD的插值方法在精度和稳定性上更具有优势，能够加速迭代过程的收敛。

相对误差的概率密度分布 两种插值方法对比

从上面的结论中可以看出，在流场中坏点较少的情况下，该方法能够很好的识别坏矢量并且能同时插值出很好的结果。为了进一步发掘该方法的优势以及测试该方法的有效性，需要建立更加苛刻的模拟条件。对基本流场添加15%的随机误差和一个5*5大小的朝一个方向的坏点区。依然按照给定的参数进行迭代。迭代7次后流场收敛，收敛过程如图XX和XX所示。能谱依然能够收敛到真实解。最终识别的坏点有15.2%，稍微高于设定的15%。

这一次，我们与归一化中值检测方法做了对比。图XX是POD识别坏点的结果，图XX

是归一化中值检测的结果。总体上看归一化检测方法只识别出13%的坏点，仍有2%坏点被漏掉。注意图中绿色圆圈内的区域，POD识别出全部的坏点，但中值检测方法在这里失效，这是因为基于相邻区域统计的方法无法计算出正确的参考向量和误差，导致判断失效。总之，对于单个坏点两种方法识别的结果相同，但对于成片的坏点，基于POD的坏点剔除明显优于中值检测。对于插值的结果，图给出了相对误差分布的云图。图XX是原始流场，图XX是POD插值后的结果。在成片误差区域相对误差较大，单点误差较小，与前面的分析一致。

原始速度场（矢量长度） 插值后速度场 相对误差

通过上面的分析，基于POD的坏点识别和其他方法4, 6相比有如下优势：

a) 该方法基于全流场的POD分解，主要流动结构的重构能够很好的揭示流动的物理

含义，不再单纯的依靠数学统计上的意义判断坏点；

b) 将[u,v]两个分量耦合求解，改进了归一化中值检测方法，更加容易识别出坏点，

对成片的坏点也能很好的识别；

c) POD插值方法使得在剔除坏点的同时能够高精度的填充坏点，得到的流场既可以直

接用作数据分析，也可以作为流场平滑的参考。

不过该方法不能用在多层互相关的程序中，比较适合对数据进行后处理。

4.实例

本节将给出用该方法识别真实PIV原始数据中坏点的例子。本实验在北京航空航天大学流体力学研究所的回流式水槽中完成，研究低速条带在层流边界层中的演化过程。条带由固定位置的圆柱干扰丝产生，在下游垂直于干扰丝放置光滑的有机玻璃板。条带进入边界层后，由于行波扰动失稳，将产生Sinuous模式和Varicose模式两种模态7。粒子图像的大小为640*480，如图XX所示。该粒子图像的粒子浓度偏低，而且出现了较多的单个像素，容易导致互相关计算时峰值锁定。从粒子图像上就可以判断互相关计算的结果应该含有较多的坏点。在实际计算的过程中，不做图像的前处理，采用了图像变形算法计算速度场，最终查询窗口的大小为16*16,50%的重叠区。每帧图像有60*80个速度矢量，总共有8497帧。为了便于显示，只取其中40*40个向量。

图XX给出本实例最重要的收敛曲线。由于层流边界层中的条带的演化是一个类周

期的过程，故前四阶模态所占的能量比较低，迭代10次后收敛到38%左右。而坏点的比例达到了13%。这么高的坏点比例，正如前面所讲，是因为粒子图像的质量不高同时又没进行图像前处理造成的。坏点识别以及插值填充的结果在图XX中给出。

针对以上问题，可以发展出一种迭代POD坏点识别程序。即在前一步识别坏点的

基础上对流场进行高精度数学插值，也可以运用 gappy POD方法对缺失的点进行高精

8, 9度重构，得到新的含坏点较少的流场。在该流场的基础上继续前一步剔除坏点的工作，循环直到流场中无法检测到坏点。

1. Hart DP. PIV error correction. Exp Fluids. 2000; 29(1): 13-22. 2. Elsinga GE, van Oudheusden BW, Scarano F. Experimental assessment of Tomographic-PIV accuracy. 2006. 3. Gao Q, Wang H, Wang J. A single camera volumetric particle image velocimetry and its application. Science China Technological Sciences. 2012. 4. Westerweel J, Scarano F. Universal outlier detection for PIV data. Exp Fluids. 2005; 39(6): 1096-100. 5. Berkooz G, Holmes P, Lumley JL. The Proper Orthogonal Decomposition in the Analysis of Turbulent Flows. Annual Review of Fluid Mechanics. 1993; 25(1): 539-75. 6. Garcia D. A fast all-in-one method for automated post-processing of PIV data. Exp Fluids. 2011; 50(5): 1247-59. 7. MATSUBARA M, ALFREDSSON PH. Disturbance growth in boundary layers subjected to free-stream turbulence. J Fluid Mech. 2001; 430: 149-68. 8. Raben SG, Charonko JJ, Vlachos PP. Adaptive gappy proper orthogonal decomposition for particle image velocimetry data reconstruction. Meas Sci Technol. 2012; 23(2). 9. Everson R, Sirovich L. Karhunen-Loève procedure for gappy data. J Opt Soc Am A. 1995; 12(8): 1657-64.