CN105095555A - 一种速度场的无散平滑处理方法及装置 - Google Patents

Abstract

本发明公开了一种速度场的无散平滑处理方法，包括：将三维速度场三个速度分量依据坐标索引展开成一个列向量速度场；构造无散光滑的标准正交基向量组；计算速度场平滑参数；根据所述无散光滑的标准正交基向量组以及计算得到的速度场平滑参数对所述列向量速度场进行无散平滑；将平滑后的无散的列向量速度场还原为与原始速度场结构相同的光滑三维速度场。本发明还公开了一种速度场的无散平滑处理装置，通过本发明，既能消除高频噪音误差，又能使平滑后速度场满足不可压缩流体的连续方程。

Description

一种速度场的无散平滑处理方法及装置

技术领域

本发明涉及流体力学速度测量技术领域，尤其涉及一种速度场的无散平滑处理方法及装置。

背景技术

粒子图像测速(PIV，ParticleImageVelocimetry)是一种现代图像测速技术，通过追踪流场中的示踪粒子来实现流场速度的测量，近年来发展起来的层析PIV(TomographicParticleImageVelocimetry)可以实现体空间内速度场三维三分量(3D3C)的测量。

在PIV实验测量中，实验数据经常因为受到噪音误差的干扰从而影响速度场的测量精度。并且，速度场中的噪音误差往往会传播到其他基于速度场计算得到的物理量中，例如：速度梯度、加速度场和压力场等。因此，在对PIV实验速度场进行分析之前，有必要对实验测得的速度场进行数据处理，以削弱噪声的影响。常见的速度场处理方法主要包括对速度矢量场进行坏点剔除、矢量场平滑、以及基于流动控制方程的速度场修正等。在这些处理中，矢量场平滑处理是消除高频噪音误差最有效的方法。但是，采用简单的数学上的平滑处理来优化速度场，往往又会引入非物理速度分量，使得平滑后的速度场不能严格满足流动的控制方程。

对于不可压缩流动的层析PIV的实验结果，即3D3C的速度场，不仅需要具有较高的精度，同时速度场光滑需要满足不可压缩流体的连续方程，即速度场无散。不可压缩流体的连续方程如式1所示：

$\▿\·U=0---\left(1\right)$

其中，U为三维速度矢量，为微分算子。

传统的平滑方法不受矢量场无散的约束，因此，无法保证平滑后速度场满足不可压缩流体的连续方程。

发明内容

为解决现有存在的技术问题，本发明实施例提供了一种速度场的无散平滑处理方法及装置，既能消除高频噪音误差，又能使平滑后速度场满足不可压缩流体的连续方程。

为达到上述目的，本发明的技术方案是这样实现的：

本发明实施例提供了一种速度场的无散平滑处理方法，所述方法包括：

将三维速度场三个速度分量依据坐标索引展开成一个列向量速度场；

构造无散光滑的标准正交基向量组；

计算速度场平滑参数；

根据所述无散光滑的标准正交基向量组以及计算得到的速度场平滑参数，对所述列向量速度场进行无散平滑；

将平滑后的无散的列向量速度场还原为与原始速度场结构相同的光滑三维速度场。

上述方案中，所述构造无散光滑的标准正交基向量组包括：

计算散度算子矩阵D；

求解不定方程组Du＝0，得到Du＝0解空间中的一组标准正交基向量其中，u为三维速度场三个速度分量依据坐标索引展开成的列向量速度场；

计算二阶偏导数系数矩阵K；

根据得到的Φ和K构造无散光滑的标准正交基向量组Ψ。

上述方案中，所述求解不定方程组Du＝0，得到Du＝0解空间里的一组标准正交基向量包括：

通过散度算子矩阵D的奇异值分解Du＝0，计算得到解空间里的一组标准正交基向量

上述方案中，所述根据得到的Φ和K构造无散光滑的标准正交基向量组Ψ包括：

求解矩阵Φ^TKΦ的特征值Λ＝diag(λ₁,λ₂,...λ_N)以及对应的特征向量矩阵P＝[p₁,p₂,...p_N]：(Φ^TKΦ)P＝PΛ；

根据特征向量矩阵P＝[p₁,p₂,...p_N]和构造无散光滑的标准正交基向量组：Ψ＝ΦP；

其中，Ψ＝[ψ₁,ψ₂,...,ψ_N]，每一列ψ_i(i＝1,2,...,N)为一个无散光滑的标准正交基向量。

上述方案中，所述计算速度场平滑参数包括：通过求解GCV函数最小值的方法来确定权重系数中待定的参数s的最佳值。

上述方案中，所述根据所述无散光滑的标准正交基向量组以及速度场平滑参数对所述列向量速度场进行无散平滑包括：

将原始速度场u_exp在无散光滑的标准正交基向量组Ψ上投影获得初始投影系数矩阵a＝Ψ^Tu_exp；

对光滑性不同的基向量对应的初始投影系数乘以不同的权重，得到修订后的投影系数矩阵其中，(I+sΛ)^-1为初始投影系数对应的权重，I为单位阵，s为速度场平滑参数，Λ为矩阵Φ^TKΦ的特征值；

根据修订后的投影系数矩阵和无散光滑的标准正交基向量组Ψ，得到平滑后的速度场u_s：u_s＝Ψ(I+sΛ)^-1Ψ^Tu_exp。

本发明实施例还提供了一种速度场的无散平滑处理装置，所述装置包括：速度场处理模块、基向量构造模块、平滑参数计算模块、速度场平滑模块，其中，

所述速度场处理模块，用于将三维速度场三个速度分量依据坐标索引展开成一个列向量速度场；

所述基向量构造模块，用于构造无散光滑的标准正交基向量组；

所述平滑参数计算模块，用于计算速度场平滑参数；

所述速度场平滑模块，用于根据所述无散光滑的标准正交基向量组以及速度场平滑参数对所述列向量速度场进行无散平滑；

所述速度场处理模块，还用于根据三维速度场展开成列向量速度场的反向操作，将平滑后的无散的列向量速度场还原为和原始速度场结构相同的光滑三维速度场。

上述方案中，所述基向量构造模块构造无散光滑的标准正交基向量组包括：

计算散度算子矩阵D；

解不定方程组Du＝0，得到Du＝0解空间里的一组标准正交基向量其中，u为三维速度场三个速度分量依据坐标索引展开成的列向量速度场；

计算二阶偏导数系数矩阵K；

根据Φ和K构造无散光滑的标准正交基向量组Ψ。

上述方案中，所述基向量构造模块解不定方程组Du＝0，得到Du＝0解空间里的一组标准正交基向量包括：

所述基向量构造模块通过散度算子矩阵D的奇异值分解Du＝0计算得到解空间里的一组标准正交基向量

上述方案中，所述基向量构造模块根据Φ和K构造无散光滑的标准正交基向量组Ψ包括：

所述基向量构造模块求解矩阵Φ^TKΦ的特征值Λ＝diag(λ₁,λ₂,...λ_N)以及对应的特征向量矩阵P＝[p₁,p₂,...p_N]：(Φ^TKΦ)P＝PΛ；

根据特征向量矩阵P＝[p₁,p₂,...p_N]和构造无散光滑的标准正交基向量组：Ψ＝ΦP；

其中，Ψ＝[ψ₁，ψ₂，...，ψ_N]，每一列ψ_i(i＝1,2,...,N)为一个无散光滑的标准正交基向量。

上述方案中，所述平滑参数计算模块计算速度场平滑参数包括：所述平滑参数计算模块通过求解GCV函数最小值的方法来确定权重系数中待定的参数s的最佳值。

上述方案中，所述速度场平滑模块根据所述无散光滑的标准正交基向量组以及速度场平滑参数对所述列向量速度场进行无散平滑包括：

所述速度场平滑模块将原始速度场u_exp在无散光滑的标准正交基向量组Ψ上投影获得初始投影系数矩阵a＝Ψ^Tu_exp；

对光滑性不同的基向量对应的初始投影系数乘以不同的权重，得到修订后的投影系数矩阵其中，(I+sΛ)^-1为初始投影系数对应的权重，I为单位阵，s为速度场平滑参数，Λ为矩阵Φ^TKΦ的特征值；

根据修订后的投影系数矩阵和无散光滑的标准正交基向量组Ψ，得到平滑后的速度场u_s：u_s＝Ψ(I+sΛ)^-1Ψ^Tu_exp。

本发明实施例提供的速度场的无散平滑处理方法及装置，将三维速度场三个速度分量依据坐标索引展开成一个列向量速度场；构造无散光滑的标准正交基向量组；计算速度场平滑参数；根据所述无散光滑的标准正交基向量组以及计算得到的速度场平滑参数对所述列向量速度场进行无散平滑；将平滑后的无散的列向量速度场还原为与原始速度场结构相同的光滑三维速度场。本发明实施例的速度场的无散平滑处理方法及装置，同时对速度场的三个分量进行平滑修正，既能消除实验数据中的高频噪音误差，又能使得平滑后的速度场严格满足不可压缩流体的连续方程。

附图说明

图1为本发明实施例速度场的无散平滑处理方法实现流程示意图；

图2为本发明实施例构造无散光滑的标准正交基向量组方法流程示意图；

图3为本发明实施例按光滑性排列的基向量的矢量场形式示意图；

图4为本发明实施例测试体的截面内速度大小在平滑去噪前后分布云图；

图5为本发明实施例平滑后速度场散度的概率密度分布图；

图6为本发明实施例速度场的无散平滑处理装置结构示意图。

具体实施方式

本发明实施例中，将三维速度场三个速度分量依据坐标索引展开成一个列向量速度场；构造无散光滑的标准正交基向量组；计算速度场平滑参数；根据所述无散光滑的标准正交基向量组以及计算得到的速度场平滑参数对所述列向量速度场进行无散平滑；将平滑后的无散的列向量速度场还原为与原始速度场结构相同的光滑三维速度场。

由于实验测得的速度场是分布在规则结构网格结点上的离散数据，这时，应该考虑公式(1)所示的连续方程的离散形式，通常连续方程的离散形式可以表示为：

Du＝0(2)

这里u是三维速度场三个速度分量依据坐标索引展开成一个列向量速度场，D是仅与网格有关的散度算子矩阵。由于方程(2)的个数小于未知数的个数，因此满足(2)的解向量有无数多个，这些解向量的全体构成(2)的解空间。本发明旨在从方程(2)解空间中找到一个比较光滑的向量来逼近原始的速度场，使用方程(2)解空间中的光滑向量来对速度场进行平滑。

具体地，首先在方程(2)的解空间内构造一组无散并且按照光滑性最优排列的标准正交基向量组；这些基向量是方程(2)的解向量；方程(2)的任意一个解向量都可以通过这些基向量表示。

然后将原始速度场向这组基向量上投影，得到各组正交基对应的初始投影系数；对不同投影系数按照基向量的光滑性施以不同的权重，进行修改以削弱噪声的影响；最后，利用修订后的投影系数与基向量重构出平滑的速度场。

投影系数权重设定时包含一个待定的参数s，控制着平滑的强弱。s越大对原始速度场的平滑越强，s越小对原始速度场的平滑越弱。为了避免对s主观的选择，本发明采用了GCV(generalizedcross-validation)方法，即通过求解GCV函数最小值的方法来确定权重系数中待定的参数s的最佳值。这使得算法对噪声强度具有较好的自适应能力，能够实现速度场有效的去噪以及连续性的修正。此外，本发明中的无散光滑基向量组仅仅与数据的网格几何特性有关，与不同瞬时速度场无关。

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

图1为本发明实施例所述的速度场的无散平滑处理方法实现流程示意图，包括：

步骤101：将三维速度场三个速度分量依据坐标索引展开成一个列向量速度场；

不可压缩流动的层析PIV测量结果是分布在规则结构网格结点上的3D3C瞬时速度场其中，为三维速度矢量,下标(i,j,k)表示网格坐标，l(l＝1,2,3)表示速度的某个分量。在每个网格节点上，速度满足方程(2)的无散条件。假设网格节点数为n＝nx×ny×nz，网格的间距为Δx,Δy,Δz。将三维速度场三个速度分量依据坐标索引展开成一个列向量速度场u后，位于向量u的第i+(j-1)×nx+(k-1)×nx×ny+(l-1)×nx×ny×nz个索引位置上，u列向量的维度为3n。原始的待平滑列向量速度场用u_exp表示，平滑后的列向量速度场用u_s表示。

步骤102：构造无散光滑的标准正交基向量组；

图2为本发明实施例构造无散光滑的标准正交基向量组方法流程示意图，如图2所示，所述构造无散光滑的标准正交基向量组的步骤包括：

步骤1021：计算散度算子矩阵D；

采用不同精度的差分格式求解的散度算子矩阵D会有不同的形式。本发明实施例以三点中心差分格式为例说明D的计算过程。

首先生成一维m节点求导矩阵：

设n阶的单位阵记为I_n。则散度算子矩阵D为：

$D=[\frac{1}{\mathrm{\Δx}}{I}_{\mathrm{nz}}\⊗I_{\mathrm{ny}}\⊗M_{\mathrm{nx}},\frac{1}{\mathrm{\Δy}}{I}_{\mathrm{nz}}\⊗M_{\mathrm{ny}}\⊗I_{\mathrm{nx}},\frac{1}{\mathrm{\Δz}}{M}_{\mathrm{nz}}\⊗I_{\mathrm{ny}}\⊗I_{\mathrm{nx}}]---\left(4\right)$

这里表示矩阵的直积(Kronecker积)，I_nz为nz阶的单位阵，I_ny为ny阶的单位阵，I_nx为nx阶的单位阵，M_nx为一维nx个节点的求导矩阵，M_ny为一维ny个节点的求导矩阵，M_nz为一维nz个节点的求导矩阵。

这里，散度算子D的维度是n×3n。

步骤1022：解不定方程组Du＝0，得到Du＝0解空间里的一组标准正交基向量其中，u为三维速度场三个速度分量依据坐标索引展开成列向量速度场；

这里，通过散度算子矩阵D的奇异值分解Du＝0计算得到解空间里的一组标准正交基向量

具体地，作D的奇异值分解：

D＝UΓV^H(5)

假设Σ＝diag(σ₁,σ₂,...σ_r),σ_i≠0。则取出V的第r+1至3n列就得到一组基向量N＝3n-r。的排列顺序可以不同。数值计算的经验表明，Φ一般包含2n+1列基向量，即N＝2n+1。

步骤1023：计算二阶偏导数系数矩阵K；

这里，首先介绍二阶偏导数系数矩阵K的意义：为了表征速度场的光滑程度，引入光滑系数的概念。光滑系数是速度场的函数，它的大小反应速度场的光滑性。将速度场的光滑系数P(u)定义为各网格点速度三个分量所有二阶偏导数之平方和，即：

$\underset{\mathrm{ijklmn}}{\mathrm{\Σ}}{\left(\frac{{\∂U}_{\mathrm{ijk}}^l}{{\∂x}_m{\∂x}_n}\right)}^2(l,m,n=\mathrm{1,2,3})$

这里记x₁,x₂,x₃分别为x,y,z。在实际计算中，速度场光滑系数同样通过线性矩阵乘积的形式来表示：

P(u)＝u^TKu

其中，K为表征速度场二阶偏导数的系数矩阵。

下面介绍二阶偏导数系数矩阵K的计算过程：

首先计算一维m个节点的二阶求导矩阵：

二阶偏导数系数矩阵K为：

$\begin{array}{c}K=\frac{1}{{\mathrm{\Δx}}^2}{I}_3\⊗I_{\mathrm{nz}}\⊗I_{\mathrm{ny}}\⊗{N_{\mathrm{nx}}}^2+\frac{1}{{\mathrm{\Δy}}^2}{I}_3\⊗I_{\mathrm{nz}}\⊗{N_{\mathrm{ny}}}^2\⊗I_{\mathrm{nx}}\\ +\frac{1}{{\mathrm{\Δz}}^2}{I}_3\⊗{N_{\mathrm{nz}}}^2\⊗I_{\mathrm{ny}}\⊗I_{\mathrm{nx}}+\frac{2}{\mathrm{\Δx\Δy}}{I}_3\⊗I_{\mathrm{nz}}\⊗N_{\mathrm{ny}}\⊗N_{\mathrm{nx}}\\ +\frac{2}{\mathrm{\Δy\Δz}}{I}_3\⊗N_{\mathrm{nz}}\⊗N_{\mathrm{ny}}\⊗I_{\mathrm{nx}}+\frac{2}{\mathrm{\Δx\Δz}}{I}_3\⊗N_{\mathrm{nz}}\⊗I_{\mathrm{ny}}\⊗N_{\mathrm{nx}}\end{array}---\left(7\right)$

其中，I_nz为nz阶的单位阵，I_ny为ny阶的单位阵，I_nx为nx阶的单位阵，M_nx为一维nx个节点的求导矩阵，M_ny为一维ny个节点的求导矩阵，M_nz为一维nz个节点的求导矩阵。

步骤1024：根据Φ和K构造无散光滑的标准正交基向量组Ψ。

求解矩阵Φ^TKΦ的特征值Λ＝diag(λ₁,λ₂,...λ_N)以及对应的特征向量矩阵P＝[p₁,p₂,...p_N]：

(Φ^TKΦ)P＝PΛ(8)

根据特征向量矩阵P＝[p₁,p₂,...p_N]和构造无散光滑的标准正交基向量组Ψ；

为了保证最后得到的基向量的光滑性有序地排列，上述特征值应该从小到大排列，即λ₁≤λ₂,...≤λ_N。最后，待求的无散光滑基向量组构成的矩阵：

Ψ＝ΦP(9)

这里Ψ＝[ψ₁，ψ₂，...，ψ_N]，中的每一列ψ_i(i＝1,2,...,N)为一个无散光滑的标准正交基向量。λ_i的大小，表征了ψ_i的光滑性。λ_i越小，ψ_i越光滑。

步骤103：计算速度场平滑参数；

具体的，通过求解GCV函数最小值的方法来确定权重系数中待定的参数s的最佳值，其中，GCV函数为：

$\mathrm{GCV}\left(s\right)=\frac{{(u_s-u_{\exp })}^T(u_s-u_{\exp })/3n}{{(1-\mathrm{tr}\left({(I+\mathrm{s\Λ})}^{-1}\right)/3n)}^2};---\left(10\right)$

其中，u_s＝Ψ(I+sΛ)^-1Ψ^Tu_exp。

可以采用数学优化的方法比如牛顿法、共轭梯度法求s。

步骤104：根据所述无散光滑的标准正交基向量组以及速度场平滑参数对所述列向量速度场进行无散平滑；

具体的，将原始速度场u_exp在无散光滑的标准正交基向量组Ψ上投影获得初始投影系数矩阵a＝Ψ^Tu_exp；

这里，向量之间的投影，通过向量的点乘实现。

初始投影系数矩阵是一个列向量，它的第i个元素a_i代表原始速度场u_exp在第i个基向量ψ_i上的投影：a_i＝ψ_i ^Tu_exp。

将光滑性不同的基向量对应的初始投影系数乘以不同的权重，得到修订后的投影系数矩阵其中，(I+sΛ)^-1为初始投影系数矩阵对应的权重系数矩阵，I为单位阵，s为速度场平滑参数，Λ为矩阵Φ^TKΦ的特征值；基向量越光滑，相应的投影系数中噪声的影响越小，施加的权重系数越大。

根据修订后的投影系数矩阵和无散光滑的标准正交基向量组Ψ，得到平滑后的速度场其中Ψ＝[ψ₁，ψ₂，...，ψ_N]，即：

u_s＝Ψ(I+sΛ)^-1Ψ^Tu_exp(11)

步骤105：将平滑后的无散的列向量速度场u_s还原为和原始速度场结构相同的光滑三维速度场。

平滑后的无散的速度场u_s为一列向量，以网格节点数为n＝nx×ny×nz为例，平滑后的无散的速度场维度为3n；根据步骤101所述的三维速度场场展开成列向量速度场的反向操作，可以获得和原始速度场相同结构的被修正后的光滑三维速度场。

下面结合具体实例验证本发明实施例所述速度场的无散平滑处理方法的可靠性。这里，用各向同性湍流直接数值模拟(DNS，DirectNumericalSimulation)结果在30×30×6等间距网格上的数据来测试。首先，在DNS的结果上加强度为0.1的高斯分布的噪声，然后用本发明的方法进行平滑。

步骤一：将三维速度场三个速度分量依据坐标索引展开成一个列向量速度场后根据102所述的步骤构造无散光滑的基向量组矩阵Ψ；

这里，取nx＝30,ny＝30,nz＝6；Δx＝Δy＝Δz＝1；用公式(3)-(9)依次计算D,Φ,K,P,Λ,Ψ。本实例中计算得到的基向量组Ψ包含10801个基向量。

图3为本发明实施例按光滑性排列的基向量的矢量场形式，其中图3-a为本发明实施例按光滑性排列的第1阶基向量的矢量场形式，图3-b为本发明实施例按光滑性排列的第5阶基向量的矢量场形式，图3-c为本发明实施例按光滑性排列的第10阶基向量的矢量场形式，图3-d为本发明实施例按光滑性排列的第20阶基向量的矢量场形式。

步骤二：通过求GCV(s)最小来计算最佳平滑参数s；

本发明实施例中，采用matlab中的fminbnd函数计算GCV(s)的最小值，得到s＝0.1568。

步骤三：利用公式(11)来对速度场进行平滑，并将平滑后的无散的列向量速度场还原为和原始速度场结构相同的光滑三维速度场。

图4为本发明实施例在测试体的截面(z＝3)内速度大小在平滑去噪前后分布云图，其中，图4-a为平滑去噪前速度分布云图，图4-b为加强度为0.1的高斯分布的噪声后的速度分布云图，图4-c为平滑去噪后速度分布云图。图5为本发明实施例平滑后速度场散度的概率密度分布，其修正后的散度在10^-5量级。散度的残差是由数值离散格式造成的，提高离散格式的精度，残差将进一步减小。

本发明实施例还提供了一种速度场的无散平滑处理装置，图6为本发明实施例速度场的无散平滑处理装置结构示意图，如图6所示，所述装置包括：速度场处理模块61、基向量构造模块62、平滑参数计算模块63、速度场平滑模块64，其中，

所述速度场处理模块61，用于将三维速度场三个速度分量依据坐标索引展开成一个列向量速度场；

这里，不可压缩流动的层析PIV测量结果是分布在规则结构网格结点上的3D3C瞬时速度场其中，为三维速度矢量,下标(i,j,k)表示网格坐标，l(l＝1,2,3)表示速度的某个分量；

将三维速度场三个速度分量依据坐标索引展开成一个列向量速度场u后，位于向量u的第i+(j-1)×nx+(k-1)×nx×ny+(l-1)×nx×ny×nz个索引位置上，u列向量的维度为3n。原始的待平滑列向量速度场用u_exp表示，平滑后的列向量速度场用u_s表示。

所述基向量构造模块62，用于构造无散光滑的标准正交基向量组；

具体的，所述基向量构造模块62构造无散光滑的标准正交基向量组的步骤包括：

a：计算散度算子矩阵D；

所述基向量构造模块62采用不同精度的差分格式求解的散度算子矩阵D会有不同的形式。本发明实施例以三点中心差分格式为例说明D的计算过程。

首先生成一维m节点求导矩阵：

设n阶的单位阵记为I_n。则散度算子矩阵D为：

$D=[\frac{1}{\mathrm{\Δx}}{I}_{\mathrm{nz}}\⊗I_{\mathrm{ny}}\⊗M_{\mathrm{nx}},\frac{1}{\mathrm{\Δy}}{I}_{\mathrm{nz}}\⊗M_{\mathrm{ny}}\⊗I_{\mathrm{nx}},\frac{1}{\mathrm{\Δz}}{M}_{\mathrm{nz}}\⊗I_{\mathrm{ny}}\⊗I_{\mathrm{nx}}]$

这里，表示矩阵的直积(Kronecker积)，I_nz为nz阶的单位阵，I_ny为ny阶的单位阵，I_nx为nx阶的单位阵，M_nx为一维nx个节点的求导矩阵，M_ny为一维ny个节点的求导矩阵，M_nz为一维nz个节点的求导矩阵。

这里，散度算子D的维度是n×3n。

b：解不定方程组Du＝0，得到Du＝0解空间里的一组标准正交基向量其中，u为三维速度场三个速度分量依据坐标索引展开成列向量速度场；

这里，通过散度算子矩阵D的奇异值分解Du＝0计算得到解空间里的一组标准正交基向量

具体地，作D的奇异值分解：

D＝UΓV^H

假设 $\mathrm{\Γ}={\left[\begin{array}{cc}{\mathrm{\Σ}}_{r\×r}& 0\\ 0& 0\end{array}\right]}_{n\×3n},$ Σ＝diag(σ₁,σ₂,...σ_r),σ_i≠0。则取出V的第r+1至3n列就得到一组基向量N＝3n-r。的排列顺序可以不同。数值计算的经验表明，Φ一般包含2n+1列基向量，即N＝2n+1。

c：计算二阶偏导数系数矩阵K；

首先计算一维m个节点的二阶求导矩阵：

二阶偏导数系数矩阵K为：

$\begin{array}{c}K=\frac{1}{{\mathrm{\Δx}}^2}{I}_3\⊗I_{\mathrm{nz}}\⊗I_{\mathrm{ny}}\⊗{N_{\mathrm{nx}}}^2+\frac{1}{{\mathrm{\Δy}}^2}{I}_3\⊗I_{\mathrm{nz}}\⊗{N_{\mathrm{ny}}}^2\⊗I_{\mathrm{nx}}\\ +\frac{1}{{\mathrm{\Δz}}^2}{I}_3\⊗{N_{\mathrm{nz}}}^2\⊗I_{\mathrm{ny}}\⊗I_{\mathrm{nx}}+\frac{2}{\mathrm{\Δx\Δy}}{I}_3\⊗I_{\mathrm{nz}}\⊗N_{\mathrm{ny}}\⊗N_{\mathrm{nx}}\\ +\frac{2}{\mathrm{\Δy\Δz}}{I}_3\⊗N_{\mathrm{nz}}\⊗N_{\mathrm{ny}}\⊗I_{\mathrm{nx}}+\frac{2}{\mathrm{\Δx\Δz}}{I}_3\⊗N_{\mathrm{nz}}\⊗I_{\mathrm{ny}}\⊗N_{\mathrm{nx}}\end{array}$

其中，I_nz为nz阶的单位阵，I_ny为ny阶的单位阵，I_nx为nx阶的单位阵，M_nx为一维nx个节点的求导矩阵，M_ny为一维ny个节点的求导矩阵，M_nz为一维nz个节点的求导矩阵。

d：根据Φ和K构造无散光滑的标准正交基向量组Ψ。

求解矩阵Φ^TKΦ的特征值Λ＝diag(λ₁,λ₂,...λ_N)以及对应的特征向量矩阵P＝[p₁,p₂,...p_N]：(Φ^TKΦ)P＝PΛ；

根据特征向量矩阵P＝[p₁,p₂,...p_N]和构造无散光滑的标准正交基向量组Ψ；

为了保证最后得到的基向量的光滑性有序地排列，上述特征值应该从小到大排列，即λ₁≤λ₂,...≤λ_N。最后，待求的无散光滑基向量组构成的矩阵即Ψ＝ΦP；

这里，Ψ＝[ψ₁，ψ₂，...，ψ_N]，中的每一列ψ_i(i＝1,2,...,N)为一个无散光滑的标准正交基向量；λ_i的大小，表征了ψ_i的光滑性。λ_i越小，ψ_i越光滑。

所述平滑参数计算模块63，用于计算速度场平滑参数；

具体的，所述平滑参数计算模块63计算速度场平滑参数包括：所述平滑参数计算模块63通过求解GCV函数最小值的方法来确定权重系数中待定的参数s的最佳值，其中，GCV函数为：

$\mathrm{GCV}\left(s\right)=\frac{{(u_s-u_{\exp })}^T(u_s-u_{\exp })/3n}{{(1-\mathrm{tr}\left({(I+\mathrm{s\Λ})}^{-1}\right)/3n)}^2};$

其中，u_s＝Ψ(I+sΛ)^-1Ψ^Tu_exp。

所述速度场平滑模块64，用于根据所述无散光滑的标准正交基向量组以及速度场平滑参数对所述列向量速度场进行无散平滑；

具体的，所述速度场平滑模块64根据所述无散光滑的标准正交基向量组以及速度场平滑参数对所述列向量速度场进行无散平滑包括：

所述速度场平滑模块64将原始速度场u_exp在无散光滑的标准正交基向量组Ψ上投影获得初始投影系数矩阵a＝Ψ^Tu_exp；

这里，向量之间的投影，通过向量的点乘实现。

初始投影系数矩阵是一个列向量，它的第i个元素a_i代表原始速度场u_exp在第i个基向量上的投影：a_i＝ψ_i ^Tu_exp。。

对光滑性不同的基向量对应的初始投影系数乘以不同的权重，得到修订后的投影系数矩阵其中，(I+sΛ)^-1为初始投影系数对应的权重，I为单位阵，s为速度场平滑参数，Λ为矩阵Φ^TKΦ的特征值；基向量越光滑，相应的投影系数中噪声的影响越小，施加的权重系数越大。

根据修订后的投影系数矩阵和无散光滑的标准正交基向量组Ψ，得到平滑后的速度场其中Ψ＝[ψ₁，ψ₂，...，ψ_N]，即：

u_s＝Ψ(I+sΛ)^-1Ψ^Tu_exp。

所述速度场处理模块61，还用于将平滑后的无散的列向量速度场还原为和原始速度场结构相同的光滑三维速度场。

平滑后的无散的速度场u_s为一列向量，以网格节点数为n＝nx×ny×nz为例，平滑后的无散的速度场维度为3n；根据三维速度场场展开成列向量速度场的反向操作，可以获得和原始速度场相同结构的被修正后的光滑三维速度场。

本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用硬件实施例、软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器和光学存储器等)上实施的计算机程序产品的形式。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

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

Claims (12)

1.一种速度场的无散平滑处理方法，其特征在于，所述方法包括：

将三维速度场三个速度分量依据坐标索引展开成一个列向量速度场；

构造无散光滑的标准正交基向量组；

计算速度场平滑参数；

根据所述无散光滑的标准正交基向量组以及计算得到的速度场平滑参数，对所述列向量速度场进行无散平滑；

将平滑后的无散的列向量速度场还原为与原始速度场结构相同的光滑三维速度场。

2.根据权利要求1所述的方法，其特征在于，所述构造无散光滑的标准正交基向量组包括：

计算散度算子矩阵D；

求解不定方程组Du＝0，得到Du＝0解空间中的一组标准正交基向量其中，u为三维速度场三个速度分量依据坐标索引展开成的列向量速度场；

计算二阶偏导数系数矩阵K；

根据得到的Φ和K构造无散光滑的标准正交基向量组Ψ。

3.根据权利要求2所述的方法，其特征在于，所述求解不定方程组Du＝0，得到Du＝0解空间里的一组标准正交基向量包括：

通过散度算子矩阵D的奇异值分解Du＝0，计算得到解空间里的一组标准正交基向量。

4.根据权利要求2所述的方法，其特征在于，所述根据得到的Φ和K构造无散光滑的标准正交基向量组Ψ包括：

求解矩阵Φ^TKΦ的特征值Λ＝diag(λ₁,λ₂,...λ_N)以及对应的特征向量矩阵P＝[p₁,p₂,...p_N]：(Φ^TKΦ)P＝PΛ；

根据特征向量矩阵P＝[p₁,p₂,...p_N]和构造无散光滑的标准正交基向量组：Ψ＝ΦP；

其中，Ψ＝[ψ₁,ψ₂,...,ψ_N]，每一列ψ_i(i＝1,2,...,N)为一个无散光滑的标准正交基向量。

5.根据权利要求1所述的方法，其特征在于，所述计算速度场平滑参数包括：通过求解GCV函数最小值的方法来确定权重系数中待定的参数s的最佳值。

6.根据权利要求1所述的方法，其特征在于，所述根据所述无散光滑的标准正交基向量组以及速度场平滑参数对所述列向量速度场进行无散平滑包括：

将原始速度场u_exp在无散光滑的标准正交基向量组Ψ上投影获得初始投影系数矩阵a＝Ψ^Tu_exp；

对光滑性不同的基向量对应的初始投影系数乘以不同的权重，得到修订后的投影系数矩阵其中，(I+sΛ)^-1为初始投影系数对应的权重，I为单位阵，s为速度场平滑参数，Λ为矩阵Φ^TKΦ的特征值；

根据修订后的投影系数矩阵和无散光滑的标准正交基向量组Ψ，得到平滑后的速度场u_s：u_s＝Ψ(I+sΛ)^-1Ψ^Tu_exp。

7.一种速度场的无散平滑处理装置，其特征在于，所述装置包括：速度场处理模块、基向量构造模块、平滑参数计算模块、速度场平滑模块，其中，

所述速度场处理模块，用于将三维速度场三个速度分量依据坐标索引展开成一个列向量速度场；

所述基向量构造模块，用于构造无散光滑的标准正交基向量组；

所述平滑参数计算模块，用于计算速度场平滑参数；

所述速度场平滑模块，用于根据所述无散光滑的标准正交基向量组以及速度场平滑参数对所述列向量速度场进行无散平滑；

所述速度场处理模块，还用于根据三维速度场展开成列向量速度场的反向操作，将平滑后的无散的列向量速度场还原为和原始速度场结构相同的光滑三维速度场。

8.根据权利要求7所述装置，其特征在于，所述基向量构造模块构造无散光滑的标准正交基向量组包括：

计算散度算子矩阵D；

解不定方程组Du＝0，得到Du＝0解空间里的一组标准正交基向量其中，u为三维速度场三个速度分量依据坐标索引展开成的列向量速度场；

计算二阶偏导数系数矩阵K；

根据Φ和K构造无散光滑的标准正交基向量组Ψ。

9.根据权利要求8所述装置，其特征在于，所述基向量构造模块解不定方程组Du＝0，得到Du＝0解空间里的一组标准正交基向量包括：

所述基向量构造模块通过散度算子矩阵D的奇异值分解Du＝0计算得到解空间里的一组标准正交基向量。

10.根据权利要求8所述装置，其特征在于，所述基向量构造模块根据Φ和K构造无散光滑的标准正交基向量组Ψ包括：

所述基向量构造模块求解矩阵Φ^TKΦ的特征值Λ＝diag(λ₁,λ₂,...λ_N)以及对应的特征向量矩阵P＝[p₁,p₂,...p_N]：(Φ^TKΦ)P＝PΛ；

根据特征向量矩阵P＝[p₁,p₂,...p_N]和构造无散光滑的标准正交基向量组：Ψ＝ΦP；

其中，Ψ＝[ψ₁,ψ₂,...,ψ_N]，每一列ψ_i(i＝1,2,...,N)为一个无散光滑的标准正交基向量。

11.根据权利要求7所述装置，其特征在于，所述平滑参数计算模块计算速度场平滑参数包括：所述平滑参数计算模块通过求解GCV函数最小值的方法来确定权重系数中待定的参数s的最佳值。

12.根据权利要求7所述装置，其特征在于，所述速度场平滑模块根据所述无散光滑的标准正交基向量组以及速度场平滑参数对所述列向量速度场进行无散平滑包括：

所述速度场平滑模块将原始速度场u_exp在无散光滑的标准正交基向量组Ψ上投影获得初始投影系数矩阵a＝Ψ^Tu_exp；

对光滑性不同的基向量对应的初始投影系数乘以不同的权重，得到修订后的投影系数矩阵其中，(I+sΛ)^-1为初始投影系数对应的权重，I为单位阵，s为速度场平滑参数，Λ为矩阵Φ^TKΦ的特征值；

根据修订后的投影系数矩阵和无散光滑的标准正交基向量组Ψ，得到平滑后的速度场u_s：u_s＝Ψ(I+sΛ)^-1Ψ^Tu_exp。