0   引　言

粒子图像测速技术（Particle Image Velocimetry,PIV）具有光路系统搭建简单、数据后处理技术相对完善、能够获取流场中的丰富结构等优势，受到了研究者关注^[1-2]并广泛应用于流动的二维^[3]、三维^[4-5]和时空^[6]结构研究。PIV的基本工作原理是在2个相邻时刻之间，通过粒子图像中有限尺度的问询域窗口的互相关求得瞬时速度场。因此，在PIV测量中，小于窗口尺度的流场信息会被过滤^[7]，导致测量结果偏离真实值。研究者对PIV有限空间分辨率导致的误差进行了研究。Saikrishnan等^[8]研究发现对原始槽道流DNS数据进行二维和三维PIV处理后，速度瞬时场云图轮廓更为平滑。Atkinson等^[9]也发现，平面PIV对槽道流中单点雷诺应力的空间滤波程度高达50%。

特别地，这种误差还体现在与速度相关的其他空间多点统计量上，尤其是需要计算速度梯度信息的统计量，例如湍流耗散率^[10-13]、应变率张量和旋转张量等。Lavoie等^[14]使用时空积分模型观察到窗口过滤作用使得PIV解析小速度和速度梯度的能力显著衰减。de Jong等^[15]对比了多种根据PIV测量结果计算湍流耗散率的方法，总结了不同方法受PIV有限空间分辨率的影响程度。Xu和Chen^[16]通过四参数的经验拟合函数定量表征了PIV问询域窗口大小对湍流耗散率测量值的影响。Buxton等^[17]通过对照直接平均的二维混合层DNS数据和三维立体PIV射流实验数据^[18]发现，较低的PIV空间分辨率会降低速度梯度张量分量的大小，并且低估剪切流中“极端事件”发生的频率，影响对湍流中精细尺度特征结构的更深入理解。

PIV有限空间分辨率导致的误差十分普遍，对其进行修正尤为必要。一个自然的思路是直接修正PIV滤波的影响，恢复瞬时速度场^[19]。Scarano^[20]提出了一种窗口变形的自适应图像处理方法，通过局部细化问询域窗口，提高PIV测量不均匀流场中大速度梯度的分辨率，以逼近真实速度场。Gao等^[21]使用神经网络优化和修正互相关计算的初始速度场，并通过反卷积层实现了空间超分辨率计算，获得了单像素速度场。但直接修正方法通常是不稳定的，会放大其他误差的影响（比如随机误差或有限像素大小带来的截断误差等）。因此，更实用的方法是对基于PIV测量结果得到的统计量进行修正。例如，Scharnowski等^[22]提出了一种由PIV互相关函数确定窗口内速度概率密度分布函数以计算雷诺应力的理论方法，并通过数值合成和实验数据验证了该方法具有较高的空间分辨率和精度。Tanaka和Eaton^[23]提出组合不同高分辨率方法来修正PIV数据中随机误差对由速度梯度计算的湍流耗散率造成的影响。

当前湍流研究越来越关注高雷诺数流动^[24-25]。随着雷诺数增大，湍流运动向更细小的尺度发展^[6]；而PIV分辨率受限于系统硬件（相机像素、激光片厚度、激光功率等）未能快速提高。为解析Kolmogorov尺度，Tanaka和Eaton^[26]使用延长筒和高分辨率镜头等设计了一种精细的高分辨率粒子图像测速系统，但却受限于较小的测量范围（仅约4 mm²）。为了获得较大尺度的流场信息，大多数情况下实验采用的PIV窗口尺寸远大于湍流耗散尺度（Kolmogorov尺度）甚至进入惯性区尺度^[27]。基于惯性子区内问询域窗口尺度的传输湍动能和亚窗口尺度的耗散湍动能之间的动态平衡假设，Sheng等^[28]提出了大涡PIV方法来估计高雷诺数下湍流耗散率的全局分布。而位于惯性区尺度的大问询域窗口内会存在明显的速度梯度，导致窗口内粒子非均匀运动十分显著。Westerweel^[29]通过理论推导和数据验证表明，当问询域窗口内示踪粒子位移的差异小于粒子图像平均直径时，流场速度梯度对PIV速度测量结果的影响才可以忽略。因此，探究大问询域窗口导致的低分辨率对湍流统计量的影响十分必要，尤其是对速度结构函数、湍流耗散率等表征湍流场特性的重要统计量的影响。

本文通过滤波模型探究PIV有限空间分辨率对均匀各向同性湍流中多点统计量的影响；理论推导得到二阶速度结构函数在耗散区、惯性区的误差一般形式，并提出修正方法；利用DNS数据生成合成PIV图像，对其进行PIV测速，检验本文误差理论及修正方法；开展近似均匀各向同性湍流中的PIV实验，并基于半经验公式对二阶速度结构函数和湍流耗散率的实验测量结果进行拟合修正。

1   PIV测量误差分析

1.1   PIV测量的误差模型

假设PIV测量中的示踪粒子完美跟随当地流体运动，忽略其他因素（例如有限像素数目、有限像素大小、粒子亮度的有限量化位数及随机误差等）的影响，通过PIV测量得到的1个问询域窗口中心点的速度$ {{\bar {\boldsymbol{u}}}}({{{\boldsymbol{x}}}}) $可视为窗口内各处流体速度的加权平均结果，即：

式中：$ {{\bar {\boldsymbol{u}}}}({{{\boldsymbol{x}}}}) $为滤波后的速度，在此代表仅考虑PIV有限问询域窗口影响的速度（后文带“–”的物理量都表示包含滤波误差的影响），${{ \boldsymbol{u}}}({\boldsymbol{x}}) $为空间点${\boldsymbol{x}} $的真实速度；${ρ} $为空间任意一点距离测量点（问询域窗口中心）的位移矢量；${G}({ρ}) $为某种形式的滤波函数。在平均意义上，${G}({ρ}) $满足3个特性：非负性，${G}({ρ}) $ ≥ 0；归一性，$ \int {G({{{ρ} }}) = 1} $；单调性，${G}({ρ}) $的值随$\left\| {\left. {ρ} \right\|} \right.$的增大而递减，即离中心越近的示踪粒子对速度测量的贡献越大。为分析更为简便，式（1）中的积分域为整个三维空间，有限窗口的大小通过${G}({ρ}) $的取值进行确定。${G}({ρ}) $的3个特性是普适的，任何实际可行的加权互相关计算均可以归结为某个合适的滤波函数${G}({ρ}) $。需要指出的是，这种以低通滤波作用来模拟PIV测量过程的思路由来已久，类似式（1）的描述也已经见诸文献^[9]。

以上是对单点速度测量滤波误差的数学描述，但多数情况下我们更关心多点速度统计量产生的误差，例如速度相关函数或与之对应的二阶速度结构函数。以二阶纵向速度结构函数为例：

式中：r为任意2个问询域中心点之间的空间间距，${\boldsymbol{e}}_1 $为单位矢量。将式（1）中滤波后的速度形式代入式（2），可得到包含PIV滤波误差的积分形式：

式中：${D_{11}}({{\boldsymbol{e}}_1}r + {{{ρ}}_2} - {{{ρ}}_1}) $表示相距$({{\boldsymbol{e}}_1}r + {{{ρ}}_2} - {{{ρ}}_1}) $、沿${\boldsymbol{e}}_1 $方向的空间任意两点的二阶速度结构函数，$ {D_{11}}({{{ρ}}_2} - {{{ρ}}_1}) $的含义类似。以$ {ρ}={{{ρ}}_2} - {{{ρ}}_1} $做变量替换，并整理式（3）为：

式中，${G^2}({{ρ}}) $为滤波函数自身卷积：

图1展示了对式（4）的直观理解：基于有限窗口大小，PIV得到的平行于${{\boldsymbol{e}}_1}r $方向的二阶纵向速度结构函数$ {\bar D_{{\mathrm{LL}}}}(r) $实际包含了在${\boldsymbol{r}}'={{\boldsymbol{e}}_1}r+{ρ} $方向的贡献（从技术细节上讲，计算$ {\bar D_{{\mathrm{LL}}}}(r) $时采用的是${{\boldsymbol{e}}_1}r $方向的速度差，因此对于沿任意${\boldsymbol{r}}' $方向的粒子对，其对$ {\bar D_{{\mathrm{LL}}}}(r) $的贡献包含平行于${\boldsymbol{r}}' $方向的速度差和垂直于${\boldsymbol{r}}' $方向的速度差，即展开式（4）后，$ {\bar D_{{\mathrm{LL}}}}(r) $的表达式中不仅会出现$ {D_{{\mathrm{LL}}}}({\boldsymbol{r'}}) $，也会出现${D_{{\mathrm{NN}}}}({\boldsymbol{r}}') $。但是，对于均匀各向同性不可压湍流，$ {D_{{\mathrm{LL}}}}({\boldsymbol{r'}}) $与${D_{{\mathrm{NN}}}}({\boldsymbol{r}}') $之间满足不可压条件约束，因此，仍然可以将式（4）写为仅含${D_{{\mathrm{LL}}}}({\boldsymbol{r'}}) $的形式）。式（4）等号右边的第2项表示图1（b）中同一问询域窗口内沿着${{\boldsymbol{e}}_1}r $方向、相距为任意${ρ} $的两点速度差的贡献，该项实际上是补偿等号右边第1项对2个示踪粒子都位于同一窗口内时重复计算的部分。

图  1  空间任意两点速度差对问询域中心点速度差的贡献

Fig.  1  The effects of velocity differences of two arbitrary points on the measurement of velocity differences of two interrogation windows

下载: 全尺寸图片 幻灯片

理论上，若已知真实的速度结构函数D_LL(r)的表达式，代入式（4）即可得到滤波后的速度结构函数$ {\bar D_{{\mathrm{LL}}}}(r) $及相应误差。然而，目前并无跨越耗散区、惯性区及大尺度区的速度结构函数解析形式，即使采用经验公式^[30]，其形式也过于复杂，不利于分析。因此，下文将分别讨论耗散区和惯性区尺度内二阶纵向速度结构函数受到的影响。

1.2   耗散区（r $\ll $η）速度结构函数误差

根据Kolmogorov推导的第一相似性假设^[31]（简称“K41理论”）可知：若两点间距在Kolmogorov特征长度η ≡ (ν³/ε)^1/4以内（其中，ν为流体的运动黏性系数，ε为湍流耗散率），以Kolmogorov特征速度u_η ≡ (εν)^1/4进行无量纲化的二阶纵向速度结构函数可写为:

式中系数a₁、a₂、a₃可表示为沿${{\boldsymbol{e}}_1}r $方向速度各阶导数的二阶矩：

系数a₁、a₂、a₃都为常数，至多与雷诺数Re_λ有关，如a₁ = 1/15^[32]、a₂ = S/(360·15^1/2)^[33]，其中，速度一阶导数的偏斜度因子$S\;{\equiv }\; \langle {( \partial {u_1}/ \partial {x_1})^3}\rangle / {\langle {( \partial {u_1}/ \partial {x_1})^2}\rangle ^{3/2}} $。

将式（6）代入式（4），即可得到耗散区受滤波影响的二阶纵向速度结构函数形式：

式中，系数A_i为关于a_i和滤波函数的表达式（i = 1, 2, 3），例如A₁可以表示为∫G²(${ρ} $)(a₁ + 6a₂${ρ} $² + 15a₃${ρ} $⁴)d${ρ} $。相应地，速度结构函数滤波前后的绝对误差δD_LL(r) = D_LL(r) − $ {\bar D_{{\mathrm{LL}}}}(r) $，则：

在式（8）和（9）中未设定滤波函数${G}({ρ}) $的形式，因此可以描述任意实际可行的PIV算法操作。式（8）和（9）的一个直接结论是：PIV系统的空间滤波并不影响二阶纵向速度结构函数在小尺度内随r²变化的标度律，仅r²的系数发生了改变。此外还可以证明，对于任一${G}({ρ}) $，滤波效应的主导因素使得系数A₁ ＜ a₁，即基于PIV结果得到的式（9）等号右边第1项r²的系数总是小于真实值。后文可以看到，这会使得利用PIV数据测量湍流耗散率产生了特定的误差。

1.3   惯性区（r$ \,\gg \,$η）速度结构函数误差

K41理论第二相似性假设^[31]为：若两点间距r远大于Kolmogorov尺度η，处于惯性区，则二阶纵向速度结构函数仅与两点间距r、湍流耗散率ε有关^[33]，可表示为无量纲形式：

式中：C₂为常数，通常取2.1^[34]。将式（10）带入式（4），即可得到惯性区内滤波后的二阶纵向速度结构函数表达式：

式（11）的第2项常数即为惯性区内滤波导致的二阶纵向速度结构函数的绝对误差：

式（12）也适用于任何形式的滤波函数。

以上为通过理论推导得出的两点速度统计量由于滤波误差在不同尺度内的表现形式。值得注意的是，二阶纵向速度结构函数在耗散区和惯性区尺度的表达式都可用于计算描述湍流的另一个重要统计量：湍流耗散率。因此，以上结果也可用于评估基于PIV得到的湍流耗散率所受滤波误差的影响。

1.4   对湍流耗散率误差的影响

湍流耗散率ε定义为：

湍流耗散率与耗散区二阶纵向速度结构函数之间的关系为：

在PIV实验中，通常先测量速度导数，再使用式（13）计算湍流耗散率，其本质与采用式（14）一致，因此，PIV测量结果中也会包含滤波误差。对比式（6）和（8）中二阶纵向速度结构函数r²前的系数，可得包含滤波误差的湍流耗散率$ \bar \varepsilon $与真实值ε的关系为：

如前所述，A₁ ＜ a₁，若不考虑其他误差，PIV测得的湍流耗散率总是小于真实值。

在PIV实验中，也可利用湍流耗散率与惯性区二阶纵向速度结构函数之间的关系进行计算：

式（16）中的速度结构函数测量值包含式（12）的常数滤波误差。常数滤波误差与基于真实速度结构函数的式（16）的计算结果之间的关系为：

由于滤波函数和二阶速度结构函数具有非负性，式（11）第2项的积分结果总是大于0，而式（17）的计算结果总是小于1，因此，滤波误差总会使PIV系统在惯性区测得的式（16）的湍流耗散率小于真实值。

1.5   特定滤波函数的影响

1.1节～1.4节的分析对任何形式的滤波函数都成立。为进一步定量显式地分析PIV有限窗口大小的影响，可以选定滤波函数，代入式（4）、（9）和（12）进行计算，得到滤波误差的具体形式。有研究者^[9]的研究思路与本文类似，采用帽型滤波函数模拟PIV的有限窗口效应，但仅考察了帽型滤波对一维能谱的影响。本文采用三维高斯滤波函数G(ρ)模拟PIV的有限窗口效应：

式中的σ可以理解为与PIV系统问询域窗口大小Δx相对应的等效滤波尺寸。对于采用同一互相关算法的PIV系统，ζ ≡ σ/Δx为描述PIV系统本身滤波特性的参数。利用PIV系统对湍流场进行测量时，其对流场参数的影响通过Δx/η反映出来。因此，将式（18）代入式（4）、（9）和（12）并积分，可得到高斯滤波模型下各统计量误差的定量表示，其中显含ζ和Δx/η参数，具体表现形式如表1所示（忽略小尺度二阶（r/η）²以上的高阶项和惯性区四阶（Δx/η）⁴以上的高阶项）。

表  1  高斯滤波函数对多点统计量的显式误差影响

Table  1  Explicit filtering errors of Gaussian filter function on measured turbulence multi-point statistics

多点 统计量	误差形式
二阶纵向结构函数绝对误差	耗散区	$ \dfrac{{{\text{δ}} {D_{{\mathrm{LL}}}}}}{{u_\eta ^2}} = \left[ { - 28{a_2}{\zeta ^2}{{\left( {\dfrac{{{\text{Δ}} x}}{\eta }} \right)}^2} - 756{a_3}{\zeta ^4}{{\left( {\dfrac{{{\text{Δ}} x}}{\eta }} \right)}^4}} \right]{\left( {\dfrac{r}{\eta }} \right)^2} $
	惯性区	$ \dfrac{{{\text{δ}} {D_{{\mathrm{LL}}}}}}{{u_\eta ^2}} = 10{a_1}{\zeta ^2}{\left( {\dfrac{{{\text{Δ}} x}}{\eta }} \right)^2} + 140{a_2}{\zeta ^4}{\left( {\dfrac{{{\text{Δ}} x}}{\eta }} \right)^4} $
二阶横向结构函数绝对误差	耗散区	$ \dfrac{{{\text{δ}} {D_{{\mathrm{NN}}}}}}{{u_\eta ^2}} = \left[ { - 56{a_2}{\zeta ^2}{{\left( {\dfrac{{{\text{Δ}} x}}{\eta }} \right)}^2} - 1512{a_3}{\zeta ^4}{{\left( {\dfrac{{{\text{Δ}} x}}{\eta }} \right)}^4}} \right]{\left( {\dfrac{r}{\eta }} \right)^2} $
	惯性区	$ \dfrac{{{\text{δ}} {D_{{\mathrm{NN}}}}}}{{u_\eta ^2}} = 8{a_1}{\zeta ^2}{\left( {\dfrac{{{\text{Δ}} x}}{\eta }} \right)^2} + 100{a_2}{\zeta ^4}{\left( {\dfrac{{{\text{Δ}} x}}{\eta }} \right)^4} $
湍流耗散率比值	$ \dfrac{{\bar \varepsilon }}{\varepsilon } = 1 + 420{a_2}{\zeta ^2}{\left( {\dfrac{{{\text{Δ}} x}}{\eta }} \right)^2} + 11340{a_3}{\zeta ^4}{\left( {\dfrac{{{\text{Δ}} x}}{\eta }} \right)^4} $

下载: 导出CSV 

| 显示表格

使用高斯函数模拟PIV计算中的滤波作用的做法值得深入讨论。首先，对于间隔时间很短的图片对，真实PIV速度测量所给出的瞬时速度取决于摄像时刻示踪粒子在图像上的位置分布。因此，对于每一图像对，若以式（1）表示测量结果，其中的瞬时等效滤波函数${G}({ρ}) $都是不同的，而正如前文引入式（1）时所言，${G}({ρ}) $应被视为平均意义上的滤波函数，这也正是本文着重分析速度统计量而非瞬时值的原因。在此情形下，很多不同的瞬时等效滤波函数平均后更趋向于高斯函数。其次，正如式（4）所揭示的，对于速度结构函数的误差而言，其中出现的并非${G}({ρ}) $，而是其自身卷积${G^2}({ρ}) $。由卷积性质可知，任一积分收敛的函数与自身不断卷积的结果会越来越接近高斯函数。例如一维平顶帽型滤波函数，其与自身卷积一次的结果为三角尖顶帽型函数，形状上已接近高斯函数。因此，在讨论速度结构函数、湍流耗散率等高阶统计量时，采用高斯函数比处理速度本身误差更为有利。最后，高斯函数具有非常优异的数学性质，采用高斯函数描述时，式（4）的积分可以简便地获得解析结果，有利于对误差的理论分析。从这个角度出发，ζ ≡ σ/Δx仅为模拟PIV系统的一个模型参数而已，若事先确定该参数的取值，后续则可在应用时对PIV系统测量误差进行先验分析或校正。

至此，本文已经通过滤波模型理论分析了PIV测量对二阶速度结构函数和湍流耗散率产生的误差影响。在将此理论结果应用于分析真实PIV系统之前，先基于数值模拟数据对其进行检验。

2   误差理论分析验证

2.1   基于DNS数据的PIV测量过程模拟

本文DNS数据来自约翰霍普金斯大学湍流数据库（Johns Hopkins Turbulence Databases, JHTDB, http://turbulence.pha.jhu.edu/）中的“Isotropic1024”算例^[35]。该算例计算网格数为1024³，计算域为(2π)³，采用周期性边界条件，通过大尺度的随机外力，维持计算域内统计稳态均匀各向同性湍流，基于泰勒微尺度定义的雷诺数Re_λ = 418，Kolmogorov特征长度为2.8 × 10⁻³或2π/(2.24 × 10³)。相关统计量如表2所示。

表  2  DNS数据库流场数据（DNS单位）

Table  2  Parameters of the DNS database (DNS unit)

泰勒 雷诺数Reλ	运动黏性 系数ν	Kolmogorov 特征长度η	Kolmogorov 特征速度uη	湍流 耗散率ε
418	1.85 × 10−4	2.8 × 10−3	6.60 × 10−2	1.03 × 10−1

下载: 导出CSV 

| 显示表格

为模拟PIV测量过程，在数据库提供的总时长约为5个大涡翻转周期的时间区间内（0～10.056）、边长为2π（DNS单位）的三维立方体空间中，等时间间隔下载了100个时刻、类似于正方形薄片的1024 × 1024 × n_z个网格点上的瞬时三维速度场，其中n_z对应的间距模拟激光片厚度Δh，该厚度也与PIV问询域窗口大小对应。该DNS湍流场积分长度尺度为1.1（即约为2π/6），因此本文将每个瞬时速度场划分为36个（6 × 6）独立的、网格点为170 × 170 × n_z的正方形子区域作为统计独立的流场区域。进行上述处理后，100个时刻的瞬时三维速度场共产生3600个进行PIV模拟的、相互独立的流场，类似于实验中在一个空间均匀的流动装置内部多个小区域进行PIV测量。为生成模拟PIV图像，在每个算例对应的瞬时时刻和170 × 170 × n_z个网格点对应的区域内随机选取示踪粒子位置。在每个示踪粒子所在位置，通过三维线性插值得到粒子速度，并根据此速度计算得到前后各dt/2时刻的粒子位置（时间间隔dt由平均粒子移动距离8～12像素确定）。在这前后2个时刻，分别生成单个粒子平均大小为3像素 × 3像素^[36]、粒子数密度为每个问询域窗口10～15个粒子、光强分布为高斯分布且同一粒子图像光强在前后2个图片上保持不变^[36]的模拟粒子图像。为考察PIV问询域窗口相对大小（或称之为“PIV分辨率”）的影响，本文选取了5种不同的PIV分辨率Δx/η（对应像素量为32像素 × 32像素或48像素 × 48像素），如表3所示。计算速度场时，每2个相邻问询域窗口之间的重合度为50%，以减小所得速度结构函数数据点之间的间距。最后，针对每一种PIV分辨率Δx/η下的3600组粒子图像对，利用开源软件PIVlab^[37]对每一对合成图像进行PIV运算，获得问询域窗口中心点的平均速度。图2展示了1对合成图像中的部分区域及计算得到的PIV速度矢量。

表  3  基于DNS数据模拟PIV图像的参数设置

Table  3  Parameter settings for DNS data to simulate PIV images

模拟激光片 相对厚度Δh/η	问询域窗口 相对大小Δx/η	问询域 窗口大小Δx	问询域窗口 之间的重合度
2	3.5	48像素 × 48像素	50%
4	4.7	32像素 × 32像素	50%
6	6.6	32像素 × 32像素	50%
8	9.1	32像素 × 32像素	50%
10	11.8	32像素 × 32像素	50%

下载: 导出CSV 

| 显示表格

图  2  基于DNS数据模拟得到的PIV合成图像及速度矢量

Fig.  2  A pair of PIV synthetic images using DNS data and diagram of velocity vectors

下载: 全尺寸图片 幻灯片

2.2   滤波误差理论分析和验证

2.2.1   二阶纵向速度结构函数

图3（a）中的离散数据点表示基于上述模拟PIV测量的速度场计算得到的二阶纵向速度结构函数${\bar D_{{\mathrm{LL}}}}(r)$，黑色实线表示基于DNS原始速度场得到的二阶纵向速度结构函数D_LL(r)。从图中可以看到，在耗散区和惯性区内，基于PIV得到的速度结构函数均比真实值偏小，在耗散区内的相对误差尤其显著，且这种误差随着问询域窗口的增大而增大。值得指出的是，包含误差的速度结构函数$ {\bar D_{{\mathrm{LL}}}}(r) $在小尺度上仍然满足r²标度律，在惯性区也近似符合r^2/3标度律，这些都与前文的理论分析结果一致。

图  3  利用DNS数据模拟PIV测量所得二阶纵向速度结构函数及其绝对误差

Fig.  3  The longitudinal structure function and its absolute errors obtained from the DNS data used to simulate PIV measure-ments

下载: 全尺寸图片 幻灯片

在图3（b）中，离散数据点表示不同PIV分辨率下二阶纵向速度结构函数的绝对误差δD_LL(r)。从图中可以看到，在小尺度耗散区内，绝对误差表现为理论公式（9）预测的主导项二次项（r/η）²，而惯性区内的绝对误差为常数，与式（12）的预测一致。

为进一步对误差进行定量分析，需确定理论模型中的滤波函数。如前所述，本文选用了高斯函数，因此需确定模型参数ζ ≡ σ/Δx。由图3（b）可知，惯性区内的绝对误差为常数，其形式较为简单，可与表1中的理论结果对比确定ζ。如图4所示，图中蓝色正方形表示不同PIV分辨率下$ {\bar D_{{\mathrm{LL}}}}(r) $在惯性区的绝对误差δD_LL(r)平均值随PIV分辨率Δx/η的变化，红色虚线表示利用表1第2行给出的δD_LL(r)在惯性区的理论形式（式（19））进行最小二乘法拟合得到的结果。

图  4  绝对误差δD_LL(r)在惯性区的平均值随分辨率的变化及拟合结果

Fig.  4  The trend of the mean values of absolute errors δD_LL(r) of second-order structure function in the inertial region with the resolution of PIV and its fitting

下载: 全尺寸图片 幻灯片

从图4可以看到，基于理论公式得到的绝对误差与模拟PIV测量得到的绝对误差吻合得相当好。采用最小二乘法拟合得到的模型参数为ζ = 0.21。实际上，当PIV分辨率较高时（即Δx/η较小时），仅保留式（19）中的（Δx/η）²项并取ζ = 0.21，即可获得令人满意的结果，如图4中的黑色曲线所示。这一特点有利于后续采用较简单的函数形式对数据进行拟合及校正。

如上所述，确定模型参数ζ后，即可利用理论结果分别考察惯性区和耗散区的误差。为定量描述δD_LL(r)在整个r/η区间的变化，受Batchelor^[30]二阶结构函数插值公式的启发，本文采用以下插值形式：

式中${{{d}}_1}$和${{{d}}_2}$取值见表4，以满足${r \mathord{\left/ {\vphantom {r \eta }} \right. } \eta } \ll 1$或${r \mathord{\left/ {\vphantom {r \eta }} \right. } \eta } \gg 1$时式（20）给出的耗散区和惯性区的绝对误差理论结果。

表  4  高斯滤波作用下的二阶速度结构函数绝对误差函数形式

Table  4  The specific forms of absolute errors of second-order structure function under the role of Gaussian filtering

多点统计量	绝对误差函数形式（忽略包含系数a3的高阶项）
二阶纵向 速度结构函数	d1	$ - 28{a_2}{\zeta ^2}{\left( {\frac{{{\text{Δ}} x}}{\eta }} \right)^2} $
	d2	$ {{{d_1}} \mathord{\left/ {\vphantom {{{d_1}} {\left[ {10{a_1}{\zeta ^2}{{\left( {\frac{{{\text{Δ}} x}}{\eta }} \right)}^2} + 140{a_2}{\zeta ^4}{{\left( {\frac{{{\text{Δ}} x}}{\eta }} \right)}^4}} \right]}}} \right. } {\left[ {10{a_1}{\zeta ^2}{{\left( {\frac{{{\text{Δ}} x}}{\eta }} \right)}^2} + 140{a_2}{\zeta ^4}{{\left( {\frac{{{\text{Δ}} x}}{\eta }} \right)}^4}} \right]}} $
二阶横向 速度结构函数	D1	$ - 56{a_2}{\zeta ^2}{\left( {\frac{{{\text{Δ}} x}}{\eta }} \right)^2} $
	D2	$ {{{D_1}} \mathord{\left/ {\vphantom {{{D_1}} {\left[ {8{a_1}{\zeta ^2}{{\left( {\frac{{{\text{Δ}} x}}{\eta }} \right)}^2} + 100{a_2}{\zeta ^4}{{\left( {\frac{{{\text{Δ}} x}}{\eta }} \right)}^4}} \right]}}} \right. } {\left[ {8{a_1}{\zeta ^2}{{\left( {\frac{{{\text{Δ}} x}}{\eta }} \right)}^2} + 100{a_2}{\zeta ^4}{{\left( {\frac{{{\text{Δ}} x}}{\eta }} \right)}^4}} \right]}} $

下载: 导出CSV 

| 显示表格

图3（b）中的实线即为式（20）给出的δD_LL(r)在不同PIV分辨率下随r/η的变化。从图中可以看到，采用插值公式（20）能够较好地描述绝对误差δD_LL(r)随r/η变化的整体趋势。

更进一步，还可以利用式（20）对PIV测量得到的${\bar D_{{\mathrm{LL}}}}(r)$进行修正并以$ {\hat D_{{\mathrm{LL}}}}(r) $近似真实值D_LL(r)：

修正结果和修正前后误差对比如图5所示。图中虚线表示修正前的直接测量值，实线表示修正后的值。从图中可以发现：PIV分辨率较高（Δx/η ＜ 7）时，整体绝对误差明显减小，由式（21）可以很好地修正滤波误差；当PIV分辨率较低（Δx/η ＞ 7）时，采用式（20）修正的效果不是很理想，其原因在于式（20）推导过程中忽略了小尺度中的四阶（Δx/η）⁴及以上的高阶项。

图  5  二阶纵向速度结构函数及其绝对误差

Fig.  5  The trend of the longitudinal second-order structure function and its absolute errors

下载: 全尺寸图片 幻灯片

需要说明的是，在以上理论推导及后续计算过程中，本文都是以二阶纵向速度结构函数为例，但这些分析方法也完全适用于二阶横向速度结构函数D_NN(r)，对于${\bar D_{{\mathrm{NN}}}}(r)$和δD_NN(r)也会得到类似的预测公式（如表1和4所示），与${\bar D_{{\mathrm{LL}}}}(r)$和δD_LL(r)的区别仅在于系数差异。

2.2.2   湍流耗散率误差分析及修正

如上节所述，理论分析结果能够较好地描述乃至修正模拟PIV测量得到的二阶速度结构函数中的误差，可以预计，理论分析结果也能较好地预测与结构函数密切相关的湍流耗散率的误差。

图6的空心方形和空心圆形数据点表示分别由式（14）和（16）得到的湍流耗散率测量值随PIV分辨率Δx/η的变化，水平实线表示基于DNS数据直接得到的湍流耗散率真实值。将经过式（21）修正后的速度结构函数惯性区数据，再利用式（16）计算得到图中实心圆形数据点，与真实值十分接近，最大相对误差仅为2%，如表5所示。

图  6  模拟PIV得到的湍流耗散率DNS数据及二阶和四阶拟合结果

Fig.  6  The variations of DNS turbulent dissipation rate obtained from the PIV calculations and the second- and fourth-order fitting curves

下载: 全尺寸图片 幻灯片

表  5  模拟PIV得到的湍流耗散率计算值与真实值之间的相对误差

Table  5  The relative errors between the calculated values and the true value of the turbulent dissipation rate obtained by simulating PIV

Δx/η	相对误差
	耗散区				惯性区
	式（14）	式（22）	式（23）		式（16）	修正后再以 式（16）计算
3.5	9.8%	3.7%	2.0%		2.6%	1.8%
4.7	19.6%	3.7%	2.0%		3.4%	2.0%
6.6	29.4%	3.7%	2.0%		4.3%	1.8%
9.1	41.5%	3.7%	2.0%		5.7%	1.3%
11.8	52.7%	3.7%	2.0%		7.9%	1.8%

下载: 导出CSV 

| 显示表格

对比可以发现，随着Δx/η增大，基于PIV测得的速度梯度（或耗散区D_LL(r)斜率）计算得到的湍流耗散率比惯性区数据更加急剧地偏离真实值。这种偏离也可以由理论预测的函数形式（表1最后一行结果）来描述：

式中的ε₀、F₂、F₃都是拟合参数。利用式（22）拟合得到的结果如图6中的红色虚线所示，可以看出该式很好地描述了PIV分辨率对湍流耗散率测量值的影响，当Δx/η = 0时，拟合预测得到的湍流耗散率ε₀与真实值相当接近，仅相差3.7%。对于Δx/η ≤ 7的情况，采用仅含二次项（Δx/η）²误差的函数，即：

代入参数a₂和ζ = 0.21进行拟合，可以得到图6中绿色点划线所示的拟合结果。从图中可以看出，对于PIV分辨率稍高的情况，式（23）对湍流耗散率测量结果的预测精度也很高。需要指出的是，利用式（23）进行拟合时，仅有1个待定参数ε₀，即湍流耗散率估计值，对于本文算例，其与真实值仅相差2.0%。以上观察提供了一种对湍流耗散率的估计思路：在PIV分辨率不是很差的情况下，可以在不同大小的问询域窗口下，采用式（22）或式（23）对湍流耗散率进行拟合，得到较好的估计值。

通过以上实例，对前文推导出的描述PIV系统测量误差的滤波模型进行了验证，并与DNS数据进行了比较，确定了高斯滤波函数下模型参数ζ的取值。将以上误差分析和修正方法应用于真实PIV实验测量数据，对其进行深入评估。

3   实验数据分析及修正

实验测量对象为冯·卡门旋流系统产生的近似均匀各向同性的湍流^[38]。实验装置如图7（a）所示。该装置主体为内径480 mm、外径636 mm、厚度10 mm的圆柱形有机玻璃筒，筒体上下两端用铝合金盖板密封。为避免图像畸变，在筒侧面安装1个平行于相机拍摄平面，直径为100 mm的圆形平面窗口；在该窗口90°方向上，安装1个相同大小的圆形平面窗口，以便进行激光照明。装置内部的工作流体（实验中为去离子水）在2个直径220 mm、高度50 mm且反向旋转的圆柱形叶轮作用下形成具有丰富结构的强湍流。以往测试表明，中心区域的流动为近似均匀各向同性湍流^[39-42]。冯·卡门旋流系统实验装置详情请参阅文献[38]第二章。

图  7  冯·卡门旋流系统流场示意图及实验布置图

Fig.  7  Schematic diagram of flow field and experimental arrangement of von Kármán swirling system

下载: 全尺寸图片 幻灯片

在叶轮转速f = 0.85 Hz、对应泰勒雷诺数Re_λ = 353时，对叶轮旋转轴对称平面中心区域的脉动速度场进行二维PIV测量。光路布置如图7（b）所示，实验参数见表6，流场拍摄范围为4.14 cm × 3.31 cm。通过PIV测量获得了400对幅面为3388像素 × 2712像素、时间间隔为320 μs的图片。采用5组大小为32像素 × 32像素、48像素 × 48像素、64像素 × 64像素、96像素 × 96像素、128像素 × 128像素的问询域窗口进行PIV互相关计算，得到5组不同PIV分辨率下的速度矢量场，对应的Δx/η分别为4.3、6.4、8.6、12.8和17.1。另外，还进行了相关统计量的计算和误差分析。

表  6  冯·卡门旋流PIV实验参数设置

Table  6  Parameter settings of von Kármán swirling PIV experiment

实验装置	描述
PIV系统 运行平台	DynamicStudio
光源设备	Nd:Yag双脉冲固体激光器 Vlite‒200，波长532 nm，脉冲能量2 × 200 mJ。示踪粒子为聚酰胺颗粒，平均直径5 μm
数据采集 设备	IMPERX B3440相机，3388像素 × 2712像素，像素尺寸3.69 μm，有效像素物理尺寸12.2 μm/像素

下载: 导出CSV 

| 显示表格

3.1   随机误差分析与修正

图8（a）展示了在不同PIV解析度下得到的二阶纵向速度结构函数${\tilde D_{{\mathrm{LL}}}}(r)$。一个明显的特点是${\tilde D_{{\mathrm{LL}}}}(r)$在小尺度未呈现出r²标度律，这是因为PIV实验数据不仅包含滤波误差，还包含不可避免的随机测量误差，而随机测量误差对二阶统计量的贡献可能在小尺度占主导^[27]。因此，为更准确评估滤波误差的作用，首先需要消除实验中随机误差的影响。

图  8  实验中随机误差修正前后对结构函数的影响

Fig.  8  The influence of random errors on longitudinal second-order structure function before and after correction in the experiments

下载: 全尺寸图片 幻灯片

记PIV系统测量得到的速度场为${{{\tilde {\boldsymbol{u}}}}_{{\mathrm{m}}}}({\boldsymbol{x}})$，与前述包含滤波误差的速度场$ {{\bar {\boldsymbol{u}}}}({\boldsymbol{x}}) $之间的关系为：

式中${\mathbf{\kappa }}({\boldsymbol{x}})$表示随机误差。测量得到的二阶纵向速度结构函数可写为：

假设每一点速度的随机误差与其他点之间统计无关，即为所谓“白噪声误差”，式（25）可化简为：

假设随机误差在各个方向的强度相等，则类似推导可得二阶横向速度结构函数为：

式（26）和（27）表明，白噪声误差对二阶速度结构函数的贡献是在包含滤波误差的基础上叠加了1个正比于误差强度的常数，这种贡献导致的随机误差在小尺度耗散区测量中尤为显著，使得当r→0时，二阶速度结构函数测量值完全由随机误差主导，逐渐偏离理论预测的r²而趋向于常数，与图8（a）中结构函数测量值的表现一致。

为消除随机误差的影响，需要估计其强度κ²。对此，可利用耗散区内不可压流场中纵向和横向二阶速度结构函数D_LL(r)与D_NN(r)之间的关系式^[32]：

易证，在滤波条件下$ {\bar D_{{\mathrm{LL}}}}(r) $和$ {\bar D_{{\mathrm{NN}}}}(r) $之间的上述关系仍然成立，因而可得：

参考表1中第1行和第3行的形式，由式（26）～（28），可以将随机误差强度表示为：

得到随机误差强度后，可对PIV原始测量值进行修正，得到去除随机误差后的PIV实验数据：

以上修正在整个r范围内进行。图8（b）中为修正后的不同分辨率下的结构函数$ {\bar D_{{\mathrm{LL}}}}(r) $，其在小尺度符合预期的r²标度律，并表现出与上节基于DNS数据模拟PIV测量结果相同的规律，即其数值随着Δx/η增大越来越小，反映了PIV计算中滤波作用的影响。下面，利用本文发展的误差理论，对通过实验得到的$ {\bar D_{{\mathrm{LL}}}}(r) $进行分析和校正（实际上，基于DNS数据模拟PIV测量的结果也包含了上述随机误差的影响，虽然强度很弱，但也可以使用上述类似方法进行修正）。

3.2   二阶速度结构函数实验误差分析及修正

为获得准确的D_LL(r)估计值，采用式（21）对实验得到的$ {\bar D_{{\mathrm{LL}}}}(r) $进行修正，消除PIV计算中滤波作用的影响，但使用该公式需要确定模型参数ζ = σ/Δx。也表明该参数具有较大的普适性。在其他研究者的PIV实验中，可能采用了不同的激光片厚度、互相关算法等，也许更希望确定对其系统更合适的参数ζ。在此，本文给出一种通过实验数据确定最优模型参数ζ_m的方法。

由理论分析式（12）可知，二阶速度结构函数在惯性区的绝对误差为常数，在分析DNS模拟数据时即可依此简单形式，利用式（19）确定参数ζ。对于实验数据，困难之处在于惯性区内的D_LL(r)真实值无法预知，不能直接采用式（19）。对此，本文给出的方案是：利用PIV分辨率最高的惯性区数据$ {\bar D_{{\mathrm{LL}}}}(r; {{{\text{Δ}} {x_1}} / \eta }= 4.3) $减去其他PIV分辨率较低的数据$ {\bar D_{{\mathrm{LL}}}}(r; {{{\text{Δ}} {x_i}} /\eta }) $，上述差值δ$ D_{{\mathrm{LL}}}^{1i}(r; {{{\text{Δ}} {x_i}} / \eta }) $在惯性区仍为常数，在选定的惯性区r/η范围内求得其平均值$ {{\text{δ}} D{_{{\mathrm{LL}}}^{1i}}_{\text{mean}}}/{u_\eta ^2}$后，由式（19）可知其与PIV分辨率Δx/η满足以下关系：

使用此式对${{\text{δ}} D{_{{\mathrm{LL}}}^{1i}}_{\text{mean}}}/{u_\eta ^2} $数据进行最小二乘拟合，可以得到针对特定实验数据的最优参数ζ_m = σ/Δx_m。对于本文实验数据，采用此方法处理后得到的拟合结果为ζ_m = 0.23。将其代入绝对误差插值函数式（20），并对实验数据$ {\bar D_{{\mathrm{LL}}}}(r) $和$ {\bar D_{{\mathrm{NN}}}}(r) $进行修正，得到如图9所示的结果。

图  9  实验测量的二阶速度结构函数修正前后对比

Fig.  9  Comparison of the second-order structure function measured by experiment before and after correction

下载: 全尺寸图片 幻灯片

从图9可以看到：当Δx/η ≤ 12时，小尺度耗散区纵向和横向结构函数由于滤波误差导致的偏离都得到了明显改善；但是，当Δx/η = 12.8和17.1时（即问询域窗口较大时），式（19）中未包含小尺度（Δx/η）⁴高阶项的影响就变得极为显著，修正效果并不理想。为修正这种情况下的误差，可以在后续研究中考虑更多的高阶项，进行更精细的理论计算。更值得指出的是：上述结果表明，这种PIV分辨率很低的情况不适宜用于研究小尺度湍流统计特性；为了对小尺度湍流进行实验探究，必须在实验时仔细考量测量系统的分辨率。本文工作为实验分辨率的选择提供了一个定量指导。

3.3   修正湍流耗散率测量值

图10中空心方形与空心圆形分别表示利用式（14）和（16）基于耗散区和惯性区结构函数$ {\bar D_{{\mathrm{LL}}}}(r) $获得的湍流耗散率随Δx/η的变化。可以发现：随着PIV问询域窗口增大，无论是耗散区还是惯性区估计值都减小，但耗散区估计值减小更为严重，这表明在实验测量之后对湍流耗散率实验数据进行修正十分重要。

图  10  冯·卡门旋流实验中湍流耗散率测量值及拟合结果

Fig.  10  The variations of measured turbulent dissipation rate with respect to PIV resolution in von Kármán swirling system and its fittings

下载: 全尺寸图片 幻灯片

采用上节方法确定参数ζ_m后，可以基于式（21）对通过实验得到的结构函数进行修正得到${\hat D_{{\mathrm{LL}}}}(r)$，再基于式（16）利用惯性区数据计算得到湍流耗散率，如图10中实心圆形所示，可以发现式（21）对惯性区的湍流耗散率有很好的修正效果。

本文将不同PIV分辨率下利用惯性区数据修正得到的湍流耗散率平均值作为真实湍流耗散率的参考值。之所以如此处理，而不是试图进一步提高PIV分辨率以通过速度梯度直接测得湍流耗散率，主要基于以下考虑：

1）通过2.2节的分析尤其是图6可知，为使得PIV空间分辨率对湍流耗散率的影响可以忽略，空间分辨率需接近Kolmogorov尺度，即Δx/η ≈ 1。如此高的空间分辨率在实验操作上有很大困难，例如需要采用类似文献[26]中的高倍率长焦镜头，还需要提高示踪粒子密度和激光功、减小激光片厚度，这些都已超出现有设备的性能范围。

2）2.2.2节利用DNS数据模拟PIV测量的结果表明，湍流耗散率在惯性区受滤波误差影响较小，且修正之后的数值与DNS数据真实值之间最大相对误差仅为1%。

基于以上考虑，本文以修正之后惯性区湍流耗散率的平均值作为真实值的参考，对小尺度下湍流耗散率测量修正后的结果进行评估。

对于受滤波影响更显著的小尺度，可以利用已知模型参数ζ_m = 0.23及式（23）对分辨率较高的数据点进行拟合，其中仅有1个拟合参数，即图10中绿色点划线与Δx/η = 0交点所示的湍流耗散率估计值ε₀（图中绿色实心方形）。ε₀的拟合结果与惯性区数据平均参考值之间的相对误差仅为1.5%。

使用式（23）进行估计时，需要确定模型参数ζ_m。本文基于惯性区数据得到该参数，但在某些情况下（如雷诺数较小、测量区域过小或流场空间不均匀性较大等），可能无法获得惯性区数据。此时若仅关心湍流耗散率，可以采用一种更为直接的方法对式（14）得到的湍流耗散率测量值进行处理。当参数ζ_m未知时，式（23）可写为：

式中有2个拟合参数ε₀、F₂。前文2.2.2节的分析表明，该式在PIV分辨率较高时，能很好地描述PIV滤波误差的影响，因此采用PIV分辨率较高（Δx/η ＜ 7）的数据，基于式（33）进行最小二乘拟合，拟合结果如图10中的红色虚线所示，其中红色实心方形为采用上述方法得到的湍流耗散率估计值ε₀。该估计值与湍流耗散率参考值之间的相对误差小于5%。因此，在真实PIV实验测量中，可以通过对同一流场的PIV图像采集结果在多组不同问询域窗口大小条件下进行互相关运算，再使用以上任一方法得到湍流耗散率小尺度的估计值。

4   结论与展望

本文采用滤波模型分析了PIV测量对单点速度矢量及相关统计量的影响，揭示了PIV测量过程对二阶速度结构函数和湍流耗散率的影响。针对高斯滤波模型获得了仅含1个模型参数ζ（ζ ≡ σ/Δx）的误差显式表达式。误差理论分析得到了基于DNS数据模拟的PIV测量结果的支持，并由此确定了模型参数ζ的具体数值。对冯·卡门旋流装置中的湍流流动开展了PIV实验测量，并对实验数据误差进行了分析和校正。研究结果表明：

1）当PIV空间分辨率Δx/η ＜ 12时，本文提出的基于高斯滤波的理论模型能够定量预测和矫正PIV测量对二阶速度结构函数及湍流耗散率的影响，相关研究结果为对在PIV实验中获得更准确的多点统计量测量值提供了基于理论指导且含有极少参数、易于使用的修正方案。

2）本文的理论分析结果为利用PIV实验研究小尺度湍流时选择PIV测量系统的分辨率提供了一个定量指导。

将本文分析方法推广至更一般流场的情形，有以下展望：就方法而言，本文的误差分析同样可用于评估PIV测量技术对其他类型湍流流场统计量的影响，其分析结果在数学形式上类似。然而，一般类型湍流流场并不具备均匀各向同性流场那样优良的数学性质，例如，对壁湍流而言（即便是充分发展的槽道壁湍流），其二阶速度结构函数不仅与两点间距有关，也与两点与壁面的距离以及两点连线与壁面的法向夹角有关，式（4）的数学形式将极为复杂。更重要的是，现有文献中关于壁湍流中二阶速度结构函数张量一般形式下的定量结果甚少，无法据之开展进一步的误差分析工作。因此，将本文工作完全照搬于壁湍流测量误差分析是不可取的，尤其是在黏性底层及缓冲区。但是，在对数区，湍流小尺度的各向同性逐渐恢复，本文中关于小尺度统计量（如湍流耗散率）的分析应该可以近似成立。进入中心区，本文对惯性区内速度结构函数测量误差的分析结果也应成立。对于平板边界层等，也可以有类似推理。至于更一般的情况，也可用本文方法分析PIV测量技术对某些特定统计量的影响，例如壁湍流对数区中沿流向方向的纵向或横向结构函数，但过程较为繁琐。

致谢：感谢西北工业大学航空学院陈鑫博士、张亿宝博士对在冯·卡门旋流流场PIV实验测量方面提供的支持和帮助。感谢国家自然科学基金委员会基础科学中心项目“非线性力学的多尺度问题研究”（项目号11988102）的支持。