CN102830121A - 一种软性磨粒流磨粒群实时检测方法 - Google Patents

Abstract

一种软性磨粒流磨粒群实时检测方法，包括：1）依照软性磨粒流磨粒群分布情况，建立相位场模型两相界面混合能量方程，同时建立相位场函数演化的控制方程，最后采用传统的二阶中心差分格式逼近粘性相和表面压力相，并用五阶WENO格式对得到的控制方程进行空间离散化重构，在时间方向采用具有TVD性质的Runge-Kutta方法进行离散，从而求解两相界面，得到磨粒群分布的理论模型；2）根据得到的磨粒群分布理论模型，实时地对软性磨粒流图像采集装置获得的磨粒群在约束流道内的分布区域进行图像处理，从而提取磨粒群分布区域轮廓特征，得到磨粒群动态边界的变化情况。本发明有效地对软性磨粒流磨粒群分布情况进行检测。

Description

一种软性磨粒流磨粒群实时检测方法

技术领域

本发明涉及软性磨粒流精密加工领域，尤其是一种软性磨粒流湍流图像实时检测方法。

背景技术

近年来，随着我国模具行业的迅猛发展，其应用范围已经涉及机械、轻工、电子、建材等多个行业。模具的传统机械加工方法大多采用数控铣或电火花加工方法进行粗加工或半精加工，但加工后仍存在加工痕迹，因此光整加工成为必不可少的重要工艺。液-固两相软性磨粒流可以在被加工工件的结构化表面形成湍流流动，它不是通过射流的形式强力冲击被加工表面，而是利用磨粒的微力微量切削的频繁作用实现表面的逐步光整，并配以约束模块，在加工工件表面形成磨粒流道，驱动磨粒对壁面的粗糙处进行切削，湍流流场中的磨粒运动的随机性实现了表面纹理无序化，直至实现结构化表面的无工具化精密加工。

软性磨粒流是一种具有弱粘性的液固两相流，在加工过程中，两相流在泵的驱动下进入约束流道，在流道内形成两相湍流，利用磨粒群的微力微量切削频繁作用于工件表面而实现逐步光整。因此的磨粒群的运动变化规律一直是研究的重点与难点，这对于实现湍流调控而言是必须克服的技术瓶颈。目前，这项技术的理论研究有了尚处在起步阶段。在实验测试方面，出现一种两相流粒子图像测速技术，以实现对多相流中流、固两相的同时测量。该技术通过在两相流场中散布能跟踪不同相流体流动的示踪粒子，用激光光源照射所测流场区域，通过在同一底片上连续曝光或在不同底片上曝光，从而获得记录示踪粒子移动的底片，最后利用多种特定的图像处理方法获得两相流场的速度分布。但是，软性磨粒流湍流流场具有结构拓扑变换特性，所得到的湍流图像单个示踪粒子无规律可循，因此需对示踪粒子分布较密集的磨粒群进行提取，分析其运动变化规律，现有的流场测试方法很难解决这一难题，更加无法快速准确地对软性磨粒流磨粒群分布情况进行实时检测与分析。

发明内容

为了克服现有的软性磨粒流磨粒群检测方法的快速性较差、准确性较差的不足，本发明提供了一种软性磨粒流磨粒群实时分析方法，可有效地对软性磨粒流磨粒群分布情况进行检测，快速准确地提取磨粒群在湍流流场中的运动变化图像。

一种软性磨粒流磨粒群实时检测方法，所述检测方法包括：

1）令Ω₁为一个二维物理空间模型，液固两相软性磨粒流为该空间的介质，磨粒群瞬时界面为分界面。假设相位场函数φ(x,t)在每个体积相为特定的常数，并且在分界面处经历快速且平滑的变化。两相流体及其交界面在t时刻的关系式由下式定义：

建立相位场模型两相界面混合能量方程：

$W(\mathrm{\φ},\▿\mathrm{\φ})={\∫}_{\mathrm{\Ω}}(\frac{1}{2}{|\▿\mathrm{\φ}|}^2+F\left(\mathrm{\φ}\right))\mathrm{dx}$

其中

$F\left(\mathrm{\φ}\right)=\frac{1}{4{\mathrm{\η}}^2}{({\left|\mathrm{\φ}\right|}^2-1)}^2$

式中F(φ)——体积能双势阱；η——交界面的宽度。

同时建立相位场函数演化的控制方程：

${\mathrm{\φ}}_t+\stackrel{\→}{u}\·\▿\mathrm{\φ}=-\mathrm{\γ}\frac{\mathrm{\δW}}{\mathrm{\δ\φ}}=\mathrm{\γ}(\mathrm{\Δ\φ}-f\left(\mathrm{\φ}\right))$

式中，f(φ)是φ的多项式，常数γ表示迁移率(m³s/kg)。

模型中两相流场的速度和压力通过Navier-Stokes方程控制，得到软性磨粒流两相流场控制方程：

$\▿\·\stackrel{\→}{u}=0$

$\stackrel{\→}{u_t}+\▿\·\left(\stackrel{\→}{u}\stackrel{\→}{u}\right)-\mathrm{\μ\Δ}\stackrel{\→}{u}+\▿p=-\mathrm{\λ\Δ\φ}\▿\mathrm{\φ}$

${\mathrm{\φ}}_t+\▿\·\left(\stackrel{\→}{u}\mathrm{\φ}\right)-\mathrm{\γ\Δ\φ}=\mathrm{\γ}(-f\left(\mathrm{\φ}\right)+\mathrm{\ξ}\left(t\right))$

$\frac{d}{\mathrm{dt}}{\∫}_{\mathrm{\Ω}}\mathrm{\φdx}=0$

式中，

——流体平均速度矢量，ρ——混合密度，p——压力，σ——应力张量（包括粘性张量和弹性应力张量），μ——流体混合粘度，λ——表面张力系数。

为张量积，其定义式为 $\▿(\▿\mathrm{\φ}\⊗\▿\mathrm{\φ})=\mathrm{\Δ\φ}\▿\mathrm{\φ}+\▿\left(\frac{1}{2}{|\▿\mathrm{\φ}|}^2\right),$

为

最后采用传统的二阶中心差分格式逼近粘性相和表面压力相，并用五阶WENO格式对得到的控制方程进行空间离散化重构，在时间方向采用具有TVD性质的Runge-Kutta方法进行离散，从而求解两相界面，得到磨粒群分布的理论模型；

2）根据得到的磨粒群分布理论模型，实时地对软性磨粒流图像采集装置获得的磨粒群在约束流道内的分布区域进行图像处理，从而提取磨粒群分布区域轮廓特征，得到磨粒群动态边界的变化情况；

所述图像处理流程如下：首先对采集的图像进行预处理，包括数字化处理，即采样和量化，图像剪切和图像增强；再利用水平集图像分割模型对磨粒群轮廓进行跟踪，分割出磨粒群形态特征；最后对图像进行后处理，提取磨粒群轮廓线。

进一步，所述对采集的图像进行预处理过程为：首先对图像进行数字化处理，即采样和量化，得到数字化图像；然后设置一个N*M的矩形剪切窗口，从数字化图像中剪切出包含整个磨粒群分布区域的子图像；再采用低通滤波增强的方法去除高频干扰噪声来平滑图像，并用Roberts算子锐化使滤波增强后图像的边缘、磨粒群轮廓线变的清晰。

再进一步，所述利用水平集图像分割模型对磨粒群轮廓进行跟踪的过程为：首先假设子图像区域为Ω，图像强度为u(x，y)：Ω→R，图像边界为水平集函数φ的零水平集Γ(t)曲线将输入图像划分为互不重叠的两个区域R₁和R₂，为防止水平集函数在其零水平集附近过于平缓或过陡，通常将水平集函数φ定义为符号距离函数：

其中，d(Γ(t),(x，y))为(x，y)到曲线Γ(t)的最短距离，从而φ(t,x，y)满足程函方程

构造水平集函数方程，使得在演化界面处水平集函数φ(X,t)保持为零，即φ(X(t),t)=0，引入正则化Heaviside函数

建立水平集图像分割能量模型：

$E(\mathrm{\φ},c_1,c_2)=v{\∫}_{\mathrm{\Ω}}|\▿H\left(\mathrm{\φ}\right)|\mathrm{dxdy}$

$+{\mathrm{\λ}}_1{\∫}_{\mathrm{\Ω}}{|u(x,y)-c_1|}^{2}H\left(\mathrm{\φ}\right)\mathrm{dxdy}$

$+{\mathrm{\λ}}_2{\∫}_{\mathrm{\Ω}}{|u(x,y)-c_2|}^2(1-H\left(\mathrm{\φ}\right))\mathrm{dxdy}$

运用变分法可得到Euler-Lagrange方程：

$\left\{\begin{array}{c}\frac{\∂\mathrm{\φ}}{\∂t}=\mathrm{\δ}\left(\mathrm{\φ}\right)[v\▿\·\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)+{\mathrm{\λ}}_2{(u-c_2)}^2-{\mathrm{\λ}}_1{(u-c_1)}^2]\mathrm{in\Ω}\\ \frac{\mathrm{\δ}\left(\mathrm{\φ}\right)}{|\▿\mathrm{\φ}|}\frac{\∂\mathrm{\φ}}{\∂n}=0\mathrm{on\Ω}\end{array}\right.$

其中Dirac函数δ(φ)＝H(φ)＇=ε/π(ε²+φ²)，c₁，c₂分别为{φ≥0}和{φ＜0}区域的灰度均值，可通过下式计算获得：

$\left\{\begin{array}{cc}{c}_1=\frac{\&Integral;\mathrm{uH}\left(\mathrm{\&phi;}\right)\mathrm{dxdy}}{\&Integral;H\left(\mathrm{\&phi;}\right)\mathrm{dxdy}}& \{\mathrm{\&phi;}\&GreaterEqual;0\}\\ c_2=\frac{\&Integral;u(1-H\left(\mathrm{\&phi;}\right))\mathrm{dxdy}}{\&Integral;(1-H\left(\mathrm{\&phi;}\right))\mathrm{dxdy}}& \{\mathrm{\&phi;}<0\}\end{array}\right.$

构造一个距离函数惩罚项，能够控制水平集演化函数保持为符号距离函数：

$P\left(\mathrm{\&phi;}\right)={\&Integral;}_{\mathrm{\&Omega;}}\frac{1}{2}{\left(\right|\&dtri;\mathrm{\&phi;}|-1)}^2\mathrm{dxdy}$

因此得到一个无需重新初始化的分割模型：

$E\left(\mathrm{\&phi;}\right)=E(\mathrm{\&phi;},c_1,c_2)+\mathrm{\&mu;P}\left(\mathrm{\&phi;}\right)$

$={\mathrm{\&lambda;}}_1{\&Integral;}_{\mathrm{\&Omega;}}{|u(x,y)-c_1|}^{2}H\left(\mathrm{\&phi;}\right)\mathrm{dxdy}+{\mathrm{\&lambda;}}_2{\&Integral;}_{\mathrm{\&Omega;}}{|u(x,y)-c_2|}^2(1-H\left(\mathrm{\&phi;}\right))\mathrm{dxdy}$

$+v{\&Integral;}_{\mathrm{\&Omega;}}|\&dtri;H\left(\mathrm{\&phi;}\right)|\mathrm{dxdy}+\mathrm{\&mu;}{\&Integral;}_{\mathrm{\&Omega;}}\frac{1}{2}{\left(\right|\&dtri;\mathrm{\&phi;}|-1)}^2\mathrm{dxdy}$

最小化能量函数E(φ)，得到相应的Euler-Lagrange方程，并由水平集方法得到曲线演化方程：

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}=\mathrm{\&delta;}\left(\mathrm{\&phi;}\right)[v\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)+{\mathrm{\&lambda;}}_2{(u-c_2)}^2-{\mathrm{\&lambda;}}_1{(u-c_1)}^2]+\mathrm{\&mu;}(\mathrm{\&Delta;\&phi;}-\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right))$

再采用水平集方法中常用的预先设置的办法人工选定初始化曲线，并初始化各参数，计算图像灰度值，将得到的曲线演化方程进行迭代求解，从而得到磨粒群分布形态特征。

更进一步，所述对图像进行后处理的过程为：采用合适的边缘检测方法提取磨粒群轮廓特征线，所述边缘检测方法为Marr边缘检测方法，先用高斯函数对图像进行平滑,然后再用拉普拉斯算子进行运算,得到Laplacian-Gauss算法,它使用一个墨西哥草帽函数形式：

$\mathrm{LOG}(x,y)=(\frac{{\&PartialD;}^2}{{\&PartialD;x}^2}+\frac{{\&PartialD;}^2}{{\&PartialD;y}^2})\frac{1}{2{\mathrm{\&pi;\&sigma;}}^2}\exp (-\frac{x^2+y^2}{{2\mathrm{\&sigma;}}^2})$

$=-\frac{1}{2{\mathrm{\&pi;\&sigma;}}^4}[2-\frac{x^2+y^2}{{\mathrm{\&sigma;}}^2}]\exp (-\frac{x^2+y^2}{{2\mathrm{\&sigma;}}^2})$

对于Log算子边缘检测的结果可以通过高斯函数标准偏差σ来进行调整，从而精确地确定磨粒群边缘轮廓线。

本发明的有益效果为：有效地对软性磨粒流磨粒群分布情况进行检测，快速准确地提取磨粒群在湍流流场中的运动变化图像。

附图说明

图1是软性磨粒流磨粒群图像轮廓特征提取流程图。

图2是Roberts算子梯度幅值计算示意图。

图3是本发明水平集图像分割流程图。

图4是水平集函数定义图。

图5是软性磨粒流图像检测装置结构示意图。

具体实施方式

下面结合附图对本发明做进一步说明。

参照图1~图5，一种软性磨粒流磨粒群实时检测方法，所述检测方法包括以下步骤：

1）令Ω₁为一个二维物理空间模型，液固两相软性磨粒流为该空间的介质，磨粒群瞬时界面为分界面。假设相位场函数φ(x,t)在每个体积相为特定的常数，并且在分界面处经历快速且平滑的变化。两相流体及其交界面在t时刻的关系式由下式定义：

建立相位场模型两相界面混合能量方程：

$W(\mathrm{\&phi;},\&dtri;\mathrm{\&phi;})={\&Integral;}_{{\mathrm{\&Omega;}}_1}(\frac{1}{2}{|\&dtri;\mathrm{\&phi;}|}^2+F\left(\mathrm{\&phi;}\right))\mathrm{dx}$

其中

$F\left(\mathrm{\&phi;}\right)=\frac{1}{4{\mathrm{\&eta;}}^2}{({\left|\mathrm{\&phi;}\right|}^2-1)}^2$

式中F(φ)——体积能双势阱；η——交界面的宽度。

同时，建立相位场函数演化的控制方程：

${\mathrm{\&phi;}}_t+\stackrel{\&RightArrow;}{u}\&CenterDot;\&dtri;\mathrm{\&phi;}=-\mathrm{\&gamma;}\frac{\mathrm{\&delta;W}}{\mathrm{\&delta;\&phi;}}=\mathrm{\&gamma;}(\mathrm{\&Delta;\&phi;}-f\left(\mathrm{\&phi;}\right))$

式中，f(φ)是φ的多项式，

常数γ表示迁移率(m³s/kg)。

模型中两相流场的速度和压力通过Navier-Stokes方程控制，其流体动量方程如下式：

$\mathrm{\&rho;}[\stackrel{\&RightArrow;}{u_t}+(\stackrel{\&RightArrow;}{u}\&CenterDot;\&dtri;)\stackrel{\&RightArrow;}{u}]=-\&dtri;p+\&dtri;\&CenterDot;\mathrm{\&sigma;}$

式中，

——流体平均速度矢量，ρ——混合密度，p——压力，σ——应力张量（包括粘性张量和弹性应力张量），由下式给出：

$\mathrm{\&sigma;}=\mathrm{\&mu;}[\&dtri;\stackrel{\&RightArrow;}{u}+{(\&dtri;\stackrel{\&RightArrow;}{u})}^T]-\mathrm{\&lambda;}\&dtri;\mathrm{\&phi;}\&CircleTimes;\&dtri;\mathrm{\&phi;}$

式中，μ——流体混合粘度，λ——表面张力系数。为张量积，其定义式为 $\&dtri;(\&dtri;\mathrm{\&phi;}\&CircleTimes;\&dtri;\mathrm{\&phi;})=\mathrm{\&Delta;\&phi;}\&dtri;\mathrm{\&phi;}+\&dtri;\left(\frac{1}{2}{|\&dtri;\mathrm{\&phi;}|}^2\right),$

为

将动量方程进一步简化，重新定义压力项：

$p=p+\frac{1}{2}\mathrm{\&lambda;}{|\&dtri;\mathrm{\&phi;}|}^2$

综上得到软性磨粒流两相流场控制方程：

$\&dtri;\&CenterDot;\stackrel{\&RightArrow;}{u}=0$

$\stackrel{\&RightArrow;}{u_t}+\&dtri;\&CenterDot;\left(\stackrel{\&RightArrow;}{u}\stackrel{\&RightArrow;}{u}\right)-\mathrm{\&mu;\&Delta;}\stackrel{\&RightArrow;}{u}+\&dtri;p=-\mathrm{\&lambda;\&Delta;\&phi;}\&dtri;\mathrm{\&phi;}$

${\mathrm{\&phi;}}_t+\&dtri;\&CenterDot;\left(\stackrel{\&RightArrow;}{u}\mathrm{\&phi;}\right)-\mathrm{\&gamma;\&Delta;\&phi;}=\mathrm{\&gamma;}(-f\left(\mathrm{\&phi;}\right)+\mathrm{\&xi;}\left(t\right))$

$\frac{d}{\mathrm{dt}}{\&Integral;}_{\mathrm{\&Omega;}}\mathrm{\&phi;dx}=0$

采用传统的二阶中心差分格式逼近粘性相和表面压力相，并用五阶WENO格式对得到的控制方程进行空间离散化重构：在时间方向采用具有TVD性质的Runge-Kutta方法进行离散，从而求解两相界面，得到磨粒群分布的理论模型。

2）根据得到的磨粒群分布理论模型，实时地对软性磨粒流图像采集装置获得的磨粒群在约束流道内的分布区域进行图像处理，从而提取磨粒群分布区域轮廓特征，得到磨粒群动态边界的变化情况。所述图像处理流程如图1所示：

2.1）首先对采集得到的图像进行预处理，包括数字化处理，即采样和量化，得到数字化图像，然后设置一个M×N的矩形剪切窗口，从数字化图像中剪切出包含整个磨粒群分布区域的子图像，再采用低通滤波增强的方法去除高频干扰噪声来平滑图像，并用Roberts算子锐化使滤波增强后图像的边缘、轮廓线以及图像的细节变的清晰。

所述采用低通滤波增强特征在于通过傅里叶变换将图像从空间域变换到频率域,然后在频率域对频谱进行操作处理,再将其反变换到空间域,从而得到增强后的图像，主要步骤为：

a)对原始图像f(x，y)进行傅里叶变换得到F(u,v)：假设图像以[f(x，y)]_M×N存储,则离散傅里叶变换[F(u,v)]_M×N可由如下公式得到：

$F(u,v)=\underset{x=0}{\overset{M-1}{\mathrm{\&Sigma;}}}\underset{y=0}{\overset{N-1}{\mathrm{\&Sigma;}}}f{(x,y)}^{-2\mathrm{i\&pi;}(\frac{\mathrm{xu}}{M}+\frac{\mathrm{yu}}{N})}$

b)将F(u,v)与传递函数H(u,v)进行卷积运算得到G(u,v)：

$G(u,v)=\underset{u=0}{\overset{M-1}{\mathrm{\&Sigma;}}}\underset{v=0}{\overset{N-1}{\mathrm{\&Sigma;}}}F(u,v)H(u,v)\frac{\mathrm{uv}}{\mathrm{MN}}$

c)将G(u,v)进行傅里叶逆变换得到增强图像g(x，y)：

$g(x,y)=\underset{u=0}{\overset{M-1}{\mathrm{\&Sigma;}}}\underset{v=0}{\overset{N-1}{\mathrm{\&Sigma;}}}f{(x,y)}^{2\mathrm{j\&pi;}(\frac{\mathrm{xu}}{M}+\frac{\mathrm{yu}}{N})}$

故低通滤波器的传递函数如下:

$H(u,v)=\left\{\begin{array}{c}1,D(u,v)\&le;D_0\\ 0,D(u,v)\&GreaterEqual;D_0\end{array}\right.$

式中,D₀是一个非负整数，D是从点(u,v)到频率平面原点的距离，即

$D(u,v)=\sqrt{u^2+v^2}$

在此低通滤波器中，小于D₀的频率,即以D₀为半径的圆内所有频率分量可以完全无损地通过,而圆外的频率,即大于D₀的频率分量则完全被除掉。

所述用Roberts算子锐化使滤波增强后图像的边缘、轮廓线以及图像的细节变的清晰，其梯度幅值近似计算方法如图2所示，(i，j)为当前像素的位置，f(i，j)为该点的灰度值,由如下锐化公式得到表示增强后的图像(i，j)位置处灰度值g(i，j)：

g(i，j)=|f(i，j)-f(i+1，j+1)|+|f(i+1，j)-f(i，j+1)|

2.2）利用水平集图像分割模型对磨粒群轮廓进行跟踪，分割出磨粒群信息，其流程如图3所示，假设输入软性磨粒流图像区域为Ω，图像强度为u(x，y)：Ω→R，图像边界为

水平集函数φ的零水平集Γ(t)曲线将输入图像划分为互不重叠的两个区域R₁和R₂，水平集函数的定义如图4所示，为防止水平集函数在其零水平集附近过于平缓或过陡，通常将水平集函数φ定义为符号距离函数：

其中，d(Γ(t),(x，y))为(x，y)到曲线Γ(t)的最短距离，从而φ(t,x，y)满足程函方程

构造水平集函数方程，使得在演化界面处水平集函数φ(X,t)保持为零，即φ(X(t),t)=0，引入正则化Heaviside函数建立水平集图像分割能量模型：

$E(\mathrm{\&phi;},c_1,c_2)=v{\&Integral;}_{\mathrm{\&Omega;}}|\&dtri;H\left(\mathrm{\&phi;}\right)|\mathrm{dxdy}$

$+{\mathrm{\&lambda;}}_1{\&Integral;}_{\mathrm{\&Omega;}}{|u(x,y)-c_1|}^{2}H\left(\mathrm{\&phi;}\right)\mathrm{dxdy}$

$+{\mathrm{\&lambda;}}_2{\&Integral;}_{\mathrm{\&Omega;}}{|u(x,y)-c_2|}^2(1-H\left(\mathrm{\&phi;}\right))\mathrm{dxdy}$

运用变分法可得到Euler-Lagrange方程：

$\left\{\begin{array}{c}\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}=\mathrm{\&delta;}\left(\mathrm{\&phi;}\right)[v\&dtri;\&CenterDot;\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)+{\mathrm{\&lambda;}}_2{(u-c_2)}^2-{\mathrm{\&lambda;}}_1{(u-c_1)}^2]\mathrm{in\&Omega;}\\ \frac{\mathrm{\&delta;}\left(\mathrm{\&phi;}\right)}{|\&dtri;\mathrm{\&phi;}|}\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;n}=0\mathrm{on\&Omega;}\end{array}\right.$

其中Dirac函数δ(φ)＝H(φ)＇=ε/π(ε²+φ²)，c₁，c₂分别为{φ≥0}和{φ＜0}区域的灰度均值，可通过下式计算获得：

$\left\{\begin{array}{cc}{c}_1=\frac{\&Integral;\mathrm{uH}\left(\mathrm{\&phi;}\right)\mathrm{dxdy}}{\&Integral;H\left(\mathrm{\&phi;}\right)\mathrm{dxdy}}& \{\mathrm{\&phi;}\&GreaterEqual;0\}\\ c_2=\frac{\&Integral;u(1-H\left(\mathrm{\&phi;}\right))\mathrm{dxdy}}{\&Integral;(1-H\left(\mathrm{\&phi;}\right))\mathrm{dxdy}}& \{\mathrm{\&phi;}<0\}\end{array}\right.$

为使水平集函数保持为符号距离函数，需要进行迭代中的初始曲线修正，即重新初始化符号距离函数，这样使计算量增大，降低了分割速度，也使分割准确度也难以保证。这里构造一个距离函数惩罚项，能够控制水平集演化函数保持为符号距离函数：

$P\left(\mathrm{\&phi;}\right)={\&Integral;}_{\mathrm{\&Omega;}}\frac{1}{2}{\left(\right|\&dtri;\mathrm{\&phi;}|-1)}^2\mathrm{dxdy}$

因此得到一个无需重新初始化的分割模型：

$E\left(\mathrm{\&phi;}\right)=E(\mathrm{\&phi;},c_1,c_2)+\mathrm{\&mu;P}\left(\mathrm{\&phi;}\right)$

$={\mathrm{\&lambda;}}_1{\&Integral;}_{\mathrm{\&Omega;}}{|u(x,y)-c_1|}^{2}H\left(\mathrm{\&phi;}\right)\mathrm{dxdy}+{\mathrm{\&lambda;}}_2{\&Integral;}_{\mathrm{\&Omega;}}{|u(x,y)-c_2|}^2(1-H\left(\mathrm{\&phi;}\right))\mathrm{dxdy}$

$+v{\&Integral;}_{\mathrm{\&Omega;}}|\&dtri;H\left(\mathrm{\&phi;}\right)|\mathrm{dxdy}+\mathrm{\&mu;}{\&Integral;}_{\mathrm{\&Omega;}}\frac{1}{2}{\left(\right|\&dtri;\mathrm{\&phi;}|-1)}^2\mathrm{dxdy}$

最小化能量函数E(φ)，得到相应的Euler-Lagrange方程，并由水平集方法得到曲线演化方程：

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}=\mathrm{\&delta;}\left(\mathrm{\&phi;}\right)[v\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)+{\mathrm{\&lambda;}}_2{(u-c_2)}^2-{\mathrm{\&lambda;}}_1{(u-c_1)}^2]+\mathrm{\&mu;}(\mathrm{\&Delta;\&phi;}-\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right))$

水平集函数在演化过程中始终保持为一个函数，利用有限差分的离散网格结构进行近似数值求解。对空间坐标变量和时间变量进行离散化，设离散网格的间隔为h，△t为时间步长，在n时刻，点(i，j)处的水平集函数为

则曲线演化方程的数值解的一种迭代形式为：

${\&dtri;}^{+}={[\max {(D_{\mathrm{ij}}^{-x},0)}^2+\min {(D_{\mathrm{ij}}^{+x},0)}^2+\max {(D_{\mathrm{ij}}^{-y},0)}^2+\min {(D_{\mathrm{ij}}^{+y},0)}^2]}^{1/2}$

${\&dtri;}^{-}={[\min {(D_{\mathrm{ij}}^{-x},0)}^2+\max {(D_{\mathrm{ij}}^{+x},0)}^2+\min {(D_{\mathrm{ij}}^{-y},0)}^2+\max {(D_{\mathrm{ij}}^{+y},0)}^2]}^{1/2}$

四个差分算子

可以表示为

$D_{\mathrm{ij}}^{-x}={\mathrm{\&phi;}}_{i,j}-{\mathrm{\&phi;}}_{i-1,j},$ $D_{\mathrm{ij}}^{+x}={\mathrm{\&phi;}}_{i+1,j}-{\mathrm{\&phi;}}_{i,j},$ $D_{\mathrm{ij}}^{-y}={\mathrm{\&phi;}}_{i,j}-{\mathrm{\&phi;}}_{i-1,j-1},$ $D_{\mathrm{ij}}^{+y}={\mathrm{\&phi;}}_{i,j+1}-{\mathrm{\&phi;}}_{i,j}$

在演化开始时，根据初始曲线计算符号距离函数，对每个网格点的水平集函数进行初始化得到零水平集，然后根据迭代来更新网格点的水平集函数，在t＝n△t时刻找到所有φ_i，j=0的点，便构成该时刻的曲线Γ(t)。因此初始曲线的选取非常重要，本发明中采用水平集方法中常用的预先设置的办法得到初始曲线位置。

为解决曲线演化方程中参数过多导致执行效率低的问题，需在模型应用前进行参数简化，首先是时间步长△t和离散网格间隔h的选取，这两个参数关系到使用有限差分方法进行偏微分方程数值求解，其次是能量函数中各参数v，λ₁，λ₂，μ的选取。时间步长△t与距离惩罚项的系数μ存在经验公式

且△t必须满足CFL条件：

能量函数的参数中v表示目标轮廓长度对拟合能量贡献的权值，λ₁，λ₂是灰度图中目标和背景同质所占的权值，μ则是距离惩罚项对能量函数贡献的权值。

图像灰度值的计算目前的技术成熟，这里不再赘述。本发明中还有一项水平集演化终止条件，假定Γ₀是当前曲线，N₀是Γ₀所包围的内部区域像素的集合，Γ_i是i次演化后的曲线，N_i是Γ_i包围的内部区域像素集合，则迭代n次后的停止准则为：

$\frac{|N_i-N_{n+i}|+|N_{n+i}-N_i|}{\frac{1}{n}\underset{j=0}{\overset{n}{\mathrm{\&Sigma;}}}{N}_{i+j}}<\mathrm{\&delta;}$

其中δ=b/N_row×N_col，N_row，N_col分别是图像的横向像素点和竖向像素点，b为某一常数。如果经过n次迭代前后，曲线位置变化很大则导致分子左边分式分子较大，准则不满足，迭代继续进行，曲线演化继续；若曲线位置基本不变，左边的分子近似为零，准则满足，迭代停止，曲线演化结束，得到磨粒群分布形态特征。

2.3）最后对图像进行后处理，其特征在于：采用合适的边缘检测方法提取磨粒群轮廓特征线，所述边缘检测方法为Marr边缘检测方法，先用高斯函数对图像进行平滑,然后再用拉普拉斯算子进行运算,得到Laplacian-Gauss算法,它使用一个墨西哥草帽函数形式：

$\mathrm{LOG}(x,y)=(\frac{{\&PartialD;}^2}{{\&PartialD;x}^2}+\frac{{\&PartialD;}^2}{{\&PartialD;y}^2})\frac{1}{2{\mathrm{\&pi;\&sigma;}}^2}\exp (-\frac{x^2+y^2}{{2\mathrm{\&sigma;}}^2})$

$=-\frac{1}{2{\mathrm{\&pi;\&sigma;}}^4}[2-\frac{x^2+y^2}{{\mathrm{\&sigma;}}^2}]\exp (-\frac{x^2+y^2}{{2\mathrm{\&sigma;}}^2})$

对于Log算子边缘检测的结果可以通过高斯函数标准偏差σ来进行调整，从而精确地确定磨粒群边缘轮廓线。

本发明所采用的技术方案是采用合适的示踪粒子跟踪磨粒群的运动，用高速CCD摄像机拍摄软性磨粒流两相流体的运动图像。由于示踪粒子和磨粒具有弱粘性，使得在两相流湍流过程中的某些时刻形成一个两相界面，因此可利用水平集图像分割方法追踪该两相动态界面的轮廓形态，保留它的拓扑变换特性，并采用合适的边缘检测方法提取磨粒群轮廓特征线。

本发明提供一种软性磨粒流图像采集装置，其结构示意图如图5所示。软性磨粒流图像检测装置依次具有氦氖激光器1、片光源2、被加工流道3、高速CCD摄像机4、图像采集卡5、控制及图像处理计算机6，所述被加工流道为透明拍摄背景一侧涂成黑色不透光材料，保证拍摄图片质量，高速CCD摄像机的光轴和片光源产生的光平面相垂直。

所述示踪粒子应具备两个基本要求：第一，粒子跟随流体的流动性好；第二，粒子是良好的散射体，其成像可见性好。

所述高速CCD摄像机的光轴和片光源产生的光平面相垂直，并且被加工流道为透明拍摄背景一侧涂成黑色不透光材料，保证拍摄质量。

在两相流场中散布合适大小和浓度的示踪粒子用于跟踪连续相流体的运动，示踪粒子在流场中分布均匀，能准确跟随流体的流动，且不改变流场特性。氦氖激光器1发出连续激光束经过片光源2形成一个扇形光片竖直向下入射到被加工流道3的上表面，流场中的示踪粒子和磨粒被激光打亮，高速CCD摄像机4拍摄流道内两相流场相变过程及其运动变化情况，通过图像采集卡5传输到图像处理计算机6进行后续的图像处理过程。

本说明书实施例所述的内容仅仅是对发明构思的实现形式的列举，本发明的保护范围的不应当被视为仅限于实施例所陈述的具体形式，本发明的保护范围也及于本领域技术人员根据本发明构思所能够想到的等同技术手段。

Claims (4)

1.一种软性磨粒流磨粒群实时检测方法，其特征在于：所述检测方法包括：

1）令Ω₁为一个二维物理空间模型，液固两相软性磨粒流为该空间的介质，磨粒群瞬时界面为分界面，假设相位场函数φ(x,t)在每个体积相为特定的常数，并且在分界面处经历快速且平滑的变化，两相流体及其交界面在t时刻的关系式由下式定义：

建立相位场模型两相界面混合能量方程： $W(\mathrm{\&phi;},\&dtri;\mathrm{\&phi;})={\&Integral;}_{\mathrm{\&Omega;}}(\frac{1}{2}{|\&dtri;\mathrm{\&phi;}|}^2+F\left(\mathrm{\&phi;}\right))\mathrm{dx}$

其中 $F\left(\mathrm{\&phi;}\right)=\frac{1}{4{\mathrm{\&eta;}}^2}{({\left|\mathrm{\&phi;}\right|}^2-1)}^2$

式中F(φ)——体积能双势阱；η——交界面的宽度；

同时建立相位场函数演化的控制方程：

${\mathrm{\&phi;}}_t+\stackrel{\&RightArrow;}{u}\&CenterDot;\&dtri;\mathrm{\&phi;}=-\mathrm{\&gamma;}\frac{\mathrm{\&delta;W}}{\mathrm{\&delta;\&phi;}}=\mathrm{\&gamma;}(\mathrm{\&Delta;\&phi;}-f\left(\mathrm{\&phi;}\right))$

式中，f(φ)是φ的多项式，

常数γ表示迁移率(m³s/kg)；

模型中两相流场的速度和压力通过Navier-Stokes方程控制，得到软性磨粒流两相流场控制方程：

$\&dtri;\&CenterDot;\stackrel{\&RightArrow;}{u}=0$

$\stackrel{\&RightArrow;}{u_t}+\&dtri;\&CenterDot;\left(\stackrel{\&RightArrow;}{u}\stackrel{\&RightArrow;}{u}\right)-\mathrm{\&mu;\&Delta;}\stackrel{\&RightArrow;}{u}+\&dtri;p=-\mathrm{\&lambda;\&Delta;\&phi;}\&dtri;\mathrm{\&phi;}$

${\mathrm{\&phi;}}_t+\&dtri;\&CenterDot;\left(\stackrel{\&RightArrow;}{u}\mathrm{\&phi;}\right)-\mathrm{\&gamma;\&Delta;\&phi;}=\mathrm{\&gamma;}(-f\left(\mathrm{\&phi;}\right)+\mathrm{\&xi;}\left(t\right))$

$\frac{d}{\mathrm{dt}}{\&Integral;}_{\mathrm{\&Omega;}}\mathrm{\&phi;dx}=0$

式中，

——流体平均速度矢量，ρ——混合密度，p——压力，σ——应力张量，μ——流体混合粘度，λ——表面张力系数，为张量积，其定义式为 $\&dtri;(\&dtri;\mathrm{\&phi;}\&CircleTimes;\&dtri;\mathrm{\&phi;})=\mathrm{\&Delta;\&phi;}\&dtri;\mathrm{\&phi;}+\&dtri;\left(\frac{1}{2}{|\&dtri;\mathrm{\&phi;}|}^2\right),$

为

最后采用传统的二阶中心差分格式逼近粘性相和表面压力相，并用五阶WENO格式对得到的控制方程进行空间离散化重构，在时间方向采用具有TVD性质的Runge-Kutta方法进行离散，从而求解两相界面，得到磨粒群分布的理论模型；

2）根据得到的磨粒群分布理论模型，实时地对软性磨粒流图像采集装置获得的磨粒群在约束流道内的分布区域进行图像处理，从而提取磨粒群分布区域轮廓特征，得到磨粒群动态边界的变化情况；

所述图像处理流程如下：首先对采集的图像进行预处理，包括数字化处理，即采样和量化，图像剪切和图像增强；再利用水平集图像分割模型对磨粒群轮廓进行跟踪，分割出磨粒群形态特征；最后对图像进行后处理，提取磨粒群轮廓线。

2.如权利要求1所述的一种软性磨粒流磨粒群实时检测方法，其特征在于：所述对采集的图像进行预处理过程为：首先对图像进行数字化处理，即采样和量化，得到数字化图像；然后设置一个N*M的矩形剪切窗口，从数字化图像中剪切出包含整个磨粒群分布区域的子图像；再采用低通滤波增强的方法去除高频干扰噪声来平滑图像，并用Roberts算子锐化使滤波增强后图像的边缘、磨粒群轮廓线变的清晰。

3.如权利要求1或2所述的一种软性磨粒流磨粒群实时检测方法，其特征在于：所述利用水平集图像分割模型对磨粒群轮廓进行跟踪的过程为：首先假设子图像区域为Ω，图像强度为u(x，y)：Ω→R，图像边界为

水平集函数φ的零水平集Γ(t)曲线将输入图像划分为互不重叠的两个区域R₁和R₂，为防止水平集函数在其零水平集附近过于平缓或过陡，通常将水平集函数φ定义为符号距离函数：

其中，d(Γ(t),(x，y))为(x，y)到曲线Γ(t)的最短距离，从而φ(t,x，y)满足程函方程

构造水平集函数方程，使得在演化界面处水平集函数φ(X,t)保持为零，即φ(X(t),t)=0，引入正则化Heaviside函数

建立水平集图像分割能量模型：

$E(\mathrm{\&phi;},c_1,c_2)=v{\&Integral;}_{\mathrm{\&Omega;}}|\&dtri;H\left(\mathrm{\&phi;}\right)|\mathrm{dxdy}$

$+{\mathrm{\&lambda;}}_1{\&Integral;}_{\mathrm{\&Omega;}}{|u(x,y)-c_1|}^{2}H\left(\mathrm{\&phi;}\right)\mathrm{dxdy}$

$+{\mathrm{\&lambda;}}_2{\&Integral;}_{\mathrm{\&Omega;}}{|u(x,y)-c_2|}^2(1-H\left(\mathrm{\&phi;}\right))\mathrm{dxdy}$

运用变分法可得到Eul er-Lagrange方程：

$\left\{\begin{array}{c}\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}=\mathrm{\&delta;}\left(\mathrm{\&phi;}\right)[v\&dtri;\&CenterDot;\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)+{\mathrm{\&lambda;}}_2{(u-c_2)}^2-{\mathrm{\&lambda;}}_1{(u-c_1)}^2]\mathrm{in\&Omega;}\\ \frac{\mathrm{\&delta;}\left(\mathrm{\&phi;}\right)}{|\&dtri;\mathrm{\&phi;}|}\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;n}=0\mathrm{on\&Omega;}\end{array}\right.$

其中Dirac函数δ(φ)＝H(φ)＇=ε/π(ε²+φ²)，c₁，c₂分别为{φ≥0}和{φ＜0}区域的灰度均值，可通过下式计算获得：

$\left\{\begin{array}{cc}{c}_1=\frac{\&Integral;\mathrm{uH}\left(\mathrm{\&phi;}\right)\mathrm{dxdy}}{\&Integral;H\left(\mathrm{\&phi;}\right)\mathrm{dxdy}}& \{\mathrm{\&phi;}\&GreaterEqual;0\}\\ c_2=\frac{\&Integral;u(1-H\left(\mathrm{\&phi;}\right))\mathrm{dxdy}}{\&Integral;(1-H\left(\mathrm{\&phi;}\right))\mathrm{dxdy}}& \{\mathrm{\&phi;}<0\}\end{array}\right.$

构造一个距离函数惩罚项，能够控制水平集演化函数保持为符号距离函数：

$P\left(\mathrm{\&phi;}\right)={\&Integral;}_{\mathrm{\&Omega;}}\frac{1}{2}{\left(\right|\&dtri;\mathrm{\&phi;}|-1)}^2\mathrm{dxdy}$

因此得到一个无需重新初始化的分割模型：

$E\left(\mathrm{\&phi;}\right)=E(\mathrm{\&phi;},c_1,c_2)+\mathrm{\&mu;P}\left(\mathrm{\&phi;}\right)$

$={\mathrm{\&lambda;}}_1{\&Integral;}_{\mathrm{\&Omega;}}{|u(x,y)-c_1|}^{2}H\left(\mathrm{\&phi;}\right)\mathrm{dxdy}+{\mathrm{\&lambda;}}_2{\&Integral;}_{\mathrm{\&Omega;}}{|u(x,y)-c_2|}^2(1-H\left(\mathrm{\&phi;}\right))\mathrm{dxdy}$

$+v{\&Integral;}_{\mathrm{\&Omega;}}|\&dtri;H\left(\mathrm{\&phi;}\right)|\mathrm{dxdy}+\mathrm{\&mu;}{\&Integral;}_{\mathrm{\&Omega;}}\frac{1}{2}{\left(\right|\&dtri;\mathrm{\&phi;}|-1)}^2\mathrm{dxdy}$

最小化能量函数E(φ)，得到相应的Euler-Lagrange方程，并由水平集方法得到曲线演化方程：

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}=\mathrm{\&delta;}\left(\mathrm{\&phi;}\right)[v\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)+{\mathrm{\&lambda;}}_2{(u-c_2)}^2-{\mathrm{\&lambda;}}_1{(u-c_1)}^2]+\mathrm{\&mu;}(\mathrm{\&Delta;\&phi;}-\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right))$

再采用水平集方法中常用的预先设置的办法人工选定初始化曲线，并初始化各参数，计算图像灰度值，将得到的曲线演化方程进行迭代求解，从而得到磨粒群分布形态特征。

4.如权利要求1或2所述的一种软性磨粒流磨粒群实时检测方法，其特征在于：所述对图像进行后处理的过程为：采用合适的边缘检测方法提取磨粒群轮廓特征线，所述边缘检测方法为Marr边缘检测方法，先用高斯函数对图像进行平滑,然后再用拉普拉斯算子进行运算,得到Laplacian-Gauss算法,它使用一个墨西哥草帽函数形式：

$\mathrm{LOG}(x,y)=(\frac{{\&PartialD;}^2}{{\&PartialD;x}^2}+\frac{{\&PartialD;}^2}{{\&PartialD;y}^2})\frac{1}{2{\mathrm{\&pi;\&sigma;}}^2}\exp (-\frac{x^2+y^2}{{2\mathrm{\&sigma;}}^2})$

$=-\frac{1}{2{\mathrm{\&pi;\&sigma;}}^4}[2-\frac{x^2+y^2}{{\mathrm{\&sigma;}}^2}]\exp (-\frac{x^2+y^2}{{2\mathrm{\&sigma;}}^2})$

对于Log算子边缘检测的结果可以通过高斯函数标准偏差σ来进行调整，从而精确地确定磨粒群边缘轮廓线。

Images