CN101629965B - Multi-grid processing method in particle image velocimetry (PIV) - Google Patents

Abstract

The invention discloses a multi-grid processing method in particle image velocimitry (PIV). The method comprises the following steps: setting a first network of the image, obtaining an equations set under the first network (corresponding to a high-resolution image) corresponding to an equations set to be solved; iteratively solving the equations set and obtaining an equation error according to a practical solution; setting a second network with a size thicker than the first network of the image; mapping a solving scale of the equations set under the first network to a solving scale of the equations set under the second network; iteratively solving a high-frequency calculation error under the second network; mapping the scale of the high-frequency calculation error under the second network to obtain a low-frequency calculation error under the first network; and updating the practical solution with the low-frequency calculation error under the first network. The method can greatly accelerate the PIV processing process and can not lower the processing reliability and the processing precision.

Description

The multi-grid processing method that is used for particle image velocimetry

Technical field

The present invention relates to Flame Image Process, particularly relate to a kind of multi-grid processing method that is used for particle image velocimetry.

Background technology

The forms of motion of nature image varies, fluid motion is a kind of typical non-rigid motion, fluid motion image calculation and analysis are that particle image velocimetry (Particle ImageVelocimitry is called for short PIV) is a kind of novel contactless measuring technique.Usually PIV drops into trace particle in fluid, under the irradiation of laser sheet optical, taking particle picture with video camera with fluid motion with sheet light vertical direction, again image is carried out Treatment Analysis and calculate, finally obtain a kind of measuring method of two-dimension speed field on the square section, flow field.As a kind of whole flow field, contactless, no disturbance, high-precision flow visual method, complicate when being applicable to turbulent flow, nonstationary flow etc. the measurement in flow field of PIV.Nowadays, PIV has been a cross-synthesis technology interdisciplinary, its combine laser technology, video image processing technology, computer technology, fluid mechanics and modern age optical technology newest fruits, be widely used in many fields such as experimental fluid mechanics, biomedicine, industry manufacturing.PIV makes as a whole continuous function with whole motion vector field and estimates, its resolution can reach each pixel in the image, for microcosmic, accurately carry out topography's motion calculation and analysis provide may, overcome some inherent shortcomings in traditional correlation technique simultaneously.All there are important science and economic worth in the PIV system in fields such as experimental fluid mechanics, aerodynamics, manufacturing industry, medical science, industry, aircraft manufacturing, water conservancy and hydropower, meteorology.

Many graphical analysis problems need solve the oval partial differential equation that has boundary value problem, and finally by find the solution that linearity or Nonlinear System of Equations realize as: based on image non-linear diffusion, initiatively profile calculating, the image recovery etc. of variational method.Though optical flow computation has all obtained bigger leap at aspects such as computational accuracy, reliabilities, its numerical method and means do not have important breakthrough.

Many grid ideas are set forth in the sixties in last century, are applied to after the seventies in the actual science calculating.Nowadays many grid computings have been widely used in science such as physics, mechanics, electromagnetic field calculates, widely used many grid methods during engineering mathematics and finite element numerical are calculated.Grid is the basis of numerical solution of partial differential equations, and the quality of grid system directly influences precision of calculation results.The research of many grid methods has been experienced from the structuring to the destructuring, the process from single grid to hybrid grid, and the new grid generation technique at different situations constantly appears.

Summary of the invention

Fundamental purpose of the present invention is exactly at the deficiencies in the prior art, and a kind of multi-grid processing method is provided, and can effectively accelerate the image processing process of particle image velocimetry.

For achieving the above object, the present invention is by the following technical solutions:

A kind of multi-grid processing method that is used for particle image velocimetry, described particle image velocimetry comprises finds the solution system of linear equations AU=F to realize the step to the optimal approximation of sports ground, U finds the solution for waiting, A is a matrix of coefficients, it comprises the motion vector field parameter at least, F is a vector, in the image processing process of particle image velocimetry, for graphical analysis, the definition of grid adopts each pixel of original image to represent refined net, based on the value computing on such discrete grid block point, the multiple dimensioned strategy that approaches is step by step adopted in the image motion field, each dimensional variation all can carry out once upgrading based on the compensation of current estimated result, the variation of graphical rule is regarded as said method comprising the steps of the variation of different size grid:

A, first grid of image is set, obtains the system of equations A of system of equations AU=F correspondence under the first grid h ^hU ^h=F ^h

B, iterative system of equations A ^hU ^h=F ^h, separate according to gained is actual

Obtain error in equation $R^h=F^h-A^h{\stackrel{\‾}{U}}^h,$

C, the size of images second grid H thicker than first grid is set, will be to system of equations A under the first grid h ^he ^h=R ^hThe yardstick of finding the solution be mapped as under the second grid H system of equations A ^He ^H=R ^HFind the solution;

D, the iterative high frequency error of calculation e under the second grid H ^H

E, to the high frequency error of calculation e under the second grid H ^HCarry out the yardstick mapping, obtain the low frequency error of calculation e under the first grid h ^h

F, with the low frequency error of calculation e under the first grid h ^hUpgrade actual separating:

Preferably, also comprise:

After upgrading actual separating, iterative A under the first grid h of image again ^hU ^h=F ^h, eliminate

The high frequency error of calculation.

Preferably, the processing of described steps A to step F repeatedly carried out in circulation, and wherein the yardstick of second grid of image increases progressively one by one.

Beneficial technical effects of the present invention is:

The graphical analysis problem need solve the oval partial differential equation that has boundary value problem, and finally realize by finding the solution linearity or Nonlinear System of Equations, as: based on image non-linear diffusion, initiatively profile calculating, the image recovery etc. of variational method, the high frequency error of linear equation is easy to eliminate on the large scale grid, even opposite to the very difficult elimination of a large amount of iteration low frequency aberrations, but this low frequency aberration only needs just can obviously eliminate through small number of iterations on the small scale grid.Therefore, in the image processing process of particle image velocimetry, adopt many grid methods, carry out the processing of finding the solution of dependent linear equation group by the iteration of different grids, can effectively accelerate the convergence process of system of linear equations iteration, accelerate image processing speed greatly, can not reduce the reliability and the precision of processing simultaneously.

The method that adopts grid multiple dimensioned (might scale factor less than 2) to approach step by step can reduce the influence that the equation model error is brought, and improves computational accuracy.Each dimensional variation all can carry out once upgrading based on the compensation of current estimated result.Regard the variation of graphical rule the variation of different size grid as, multiple dimensioned in essence approaching calculated a kind of special case that can be regarded as multi-grid algorithm.

Description of drawings

Fig. 1 is the structured flowchart of particle image velocimetry disposal system;

Fig. 2 is for dividing the interior smooth synoptic diagram of discontinuous field, border continuously with level set function;

The process flow diagram of Fig. 3 for utilizing many grids to quicken to handle among a kind of embodiment;

Fig. 4 a and Fig. 4 b are respectively the synoptic diagram of V-type and the many networks of W type;

Fig. 5 is a kind of synoptic diagram of many networks of preferred embodiment.

Feature of the present invention and advantage will be elaborated in conjunction with the accompanying drawings by embodiment.

Embodiment

At first introduce a kind of particle image velocimetry method of the multi-grid processing method that can use the embodiment of the invention.

A complete particle image velocimetry disposal system key component as shown in Figure 1, mainly comprise data acquisition module, image pretreatment module, motion vector field (optical flow field) computing module and follow-up data analysis module, the follow-up data analysis module comprises curl computing module and divergence computing module.Motion vector field computing module wherein is the nucleus module of whole PIV system.

When calculating 2 dimension motion vector field, divergence and curl fields, adopted calculating thought, with the form calculating kinematical vector of whole function based on the light stream equation of motion.Can keep simultaneously the partial structurtes information of vector field, divergence and curl field again in order to restrain noise, also combine the thought of level set.

2 dimension motion vector fields are meant in size is N * M pixel coverage (as the PAL-system of knowing usually: 768 * 576 or 640 * 480 pixels) can describe the displacement vector of this pixel current time under each location of pixels (in N * M), this vector is a floating number two-dimensional vector, is expressed as: v=[v _x, v _y], v wherein _x, v _yBe x axle and axial two speed components of y.Describe for unified, represent vector with the capitalization runic, the small letter light face type is represented scalar.Velocity field v is actual to be coordinate position variable x=[x, y] function.The fluid motion image calculation has gone out outside the acquisition transient motion vector field v, also need calculate divergence field ξ and the curl field η of v, shown in the formula of (1) (2).

$\mathrm{\ξ}=\mathrm{div}\left(V\right)=\frac{{\∂v}_x}{\∂x}+\frac{{\∂v}_y}{\∂y}---\left(1\right)$

$\mathrm{\η}=\mathrm{curl}\left(V\right)=\frac{{\∂v}_x}{\∂y}-\frac{{\∂v}_y}{\∂x}---\left(2\right)$

Because actual fluid motion may have the phenomenon that many vector currents have a common boundary, be the moving boundaries and the problem of doing more physical exercises, in order to solve the moving boundaries and the problem of doing more physical exercises,, the strategy and the thought of level set is integrated in original optical flow computation framework simultaneously in order to improve the precision of fluid motion feature description.

As shown in Figure 2, motion vector, divergence and curl might be discontinuous in a certain regional Ω, Ω can be divided into two sub regions Ω for this reason ₁, Ω ₂, Ω ₁Sports ground, divergence field and curl field corresponding in the zone are: [v ₁, ξ ₁, η ₁]; Ω ₂Sports ground, divergence field and curl field corresponding in the zone are: [v ₂, ξ ₂, η ₂].Suppose at regional Ω ₁Interior [v ₁, ξ ₁, η ₁] smooth and continuous, at regional Ω ₂Interior [v ₂, ξ ₂, η ₂] smooth and continuous, but at Ω ₁, Ω ₂Intersection sports ground, divergence field and curl field and discontinuous.[v in order to guarantee to calculate ₁, ξ ₁, η ₁] and [v ₂, ξ ₂, η ₂] continuity, do not destroy simultaneously Ω again ₁With Ω ₂Between structural information and be convenient to calculate, make up level set function φ and indicator function H, φ is for being defined in the image magnitude range the (continuous function of N * M).H(x)＝1，x≥0and?H(x)＝0，x＜O。

φ divides Ω with level set function ₁, Ω ₂, the zone of φ 〉=0 belongs to Ω ₁, and the zone of φ＜0 belongs to Ω ₂Owing to adopted the thought of level set, therefore calculated [v ₁, ξ ₁, η ₁] and [v ₂, ξ ₂, η ₂] time will constantly calculate simultaneously and upgrade φ.

In order to estimate [v ₁, ξ ₁, η ₁] and [v ₂, ξ ₂, η ₂], construct following energy function.

E＝E ₁+λE ₂ (3)

Wherein E is the global energy functional, E ₁, E ₂Be respectively data constraint energy and smooth bound energy with following formula (4), (5) form.

$E_1=\underset{\mathrm{\&Omega;}}{\&Integral;}|I(X+V_1,t)-I(X,t){|}^{2}H\left(\mathrm{\&phi;}\right)+|I(X+{\mathrm{<figure-callout\; id="2"\; label="V"\; filenames="H2009101094293E00012.png"\; state="\{\{state\}\}">V</figure-callout>}}_2,t)-I(X,t){|}^{2}H(-\mathrm{\&phi;})\mathrm{dx}$

$+\underset{\mathrm{\&Omega;}}{\&Integral;}|\mathrm{div}\left(V_1\right)-{\mathrm{\&xi;}}_1{|}^{2}H\left(\mathrm{\&phi;}\right)+|\mathrm{div}\left(V_2\right)-{\mathrm{\&xi;}}_2{|}^{2}H(-\mathrm{\&phi;})\mathrm{dx}---\left(4\right)$

$+\underset{\mathrm{\&Omega;}}{\&Integral;}|\mathrm{curl}\left(V_1\right)-{\mathrm{\&eta;}}_1{|}^{2}H\left(\mathrm{\&phi;}\right)+|\mathrm{curl}\left(V_2\right)-{\mathrm{\&eta;}}_2{|}^{2}H(-\mathrm{\&phi;})\mathrm{dx}$

${\mathrm{<figure-callout\; id="2"\; label="E"\; filenames="H2009101094293E00012.png"\; state="\{\{state\}\}">E</figure-callout>}}_2=\mathrm{\&alpha;}\underset{\mathrm{\&Omega;}}{\&Integral;}|\&dtri;V_1{|}^{2}H\left(\mathrm{\&phi;}\right)+|\&dtri;{\mathrm{<figure-callout\; id="2"\; label="V"\; filenames="H2009101094293E00012.png"\; state="\{\{state\}\}">V</figure-callout>}}_2{|}^{2}H(-\mathrm{\&phi;})\mathrm{dx}$

$+\mathrm{\&alpha;}\underset{\mathrm{\&Omega;}}{\&Integral;}|\&dtri;{\mathrm{\&xi;}}_1{|}^{2}H\left(\mathrm{\&phi;}\right)+|\&dtri;{\mathrm{\&xi;}}_2{|}^{2}H(-\mathrm{\&phi;})\mathrm{dx}$

(5)

$+\mathrm{\&alpha;}\underset{\mathrm{\&Omega;}}{\&Integral;}|\&dtri;{\mathrm{\&eta;}}_1{|}^{2}H\left(\mathrm{\&phi;}\right)+|\&dtri;{\mathrm{\&eta;}}_2{|}^{2}H(-\mathrm{\&phi;})\mathrm{dx}$

$+\mathrm{\&beta;}\underset{\mathrm{\&Omega;}}{\&Integral;}|\&dtri;H\left(\mathrm{\&phi;}\right){|}^2\mathrm{dx}$

Here (X t) represents the particle picture sequence to I.The final optical flow field of estimating can be described as:

$V=\left\{\begin{array}{cc}{V}_1& \mathrm{\&phi;}\left(X\right)\&GreaterEqual;0\\ {\mathrm{<figure-callout\; id="2"\; label="V"\; filenames="H2009101094293E00012.png"\; state="\{\{state\}\}">V</figure-callout>}}_2& \mathrm{else}\end{array}\right.$

Realize the optimal approximation of optical flow field by minimizing of global energy functional.

Minimizing of energy functional (3) can realize by finding the solution system of linear equations (6).

AU＝F (6)

U=[V wherein ₁, V ₂, ξ ₁, ξ ₂, η ₁, η ₂, φ], A is 9 * 9 matrix of coefficients, F is the odd item of 9 * 1 vector representation linear equation.

The acquisition of system of linear equations (6) and calculating can realize by classical Variational Calculation method, not elaborate here.

For above-mentioned PIV method, can adopt many grids acceleration strategy to combine with it, particularly, promptly handle system of linear equations (6) by many grid methods, quicken PIV and calculate.

A central task of many grid computings is to find the solution various partial differential equation and system of linear equations from practical problems.Because many graphical analysis problems need solve the oval partial differential equation that has boundary value problem, and finally realize by finding the solution linearity or Nonlinear System of Equations, as: based on image non-linear diffusion, initiatively profile calculating, the image recovery etc. of variational method, the high frequency error of linear equation is easy to eliminate on the large scale grid, even opposite to the very difficult elimination of a large amount of iteration low frequency aberrations, but this low frequency aberration only needs just can obviously eliminate through small number of iterations on the small scale grid.Therefore many grid. policies are exactly to come the accelerating convergence process by the iterative computation of different grids.For graphical analysis, the definition of grid is fairly simple, and each pixel of original image is represented a grid, and promptly the thinnest grid calculates that the value computing all be based on such discrete grid block point finishes.

Please refer to Fig. 3, the multi-grid processing method of a kind of embodiment specifically comprises the steps:

At size of mesh opening is h _x* n _yThe first grid h under, formula (6) AU=F can be written as formula (7):

A ^hU ^h＝F ^h (7)

A ^h, U ^hBeing illustrated respectively in size of mesh opening is h _x* n _yEquation coefficient matrix under (image resolution ratios that different size of mesh opening are corresponding different) is conciliate.

Suppose that equation (7) is through n ₁Obtain after the inferior iteration separate for

, and make U ^hBe desirable truly separating, then the error of calculation of motion vector is:

$e^h=U^h-{\stackrel{\&OverBar;}{U}}^h---\left(8\right)$

Because the high frequency error of result of calculation just can be eliminated substantially through small number of iterations, so e ^hThe low frequency energy that mainly comprises the error of calculation here.If can know e ^hThen can be right Compensate and reach the purpose that reduces the error of calculation.Though directly obtain e ^hAnd be not easy, but can through type (9) and formula (10) find the solution indirectly.

By known

Can calculate the error R of equation (7) ^h, formula (9) is arranged.

$R^h=F^h-A^h{\stackrel{\&OverBar;}{U}}^h$

(9)

Because A ^hBe linear therefore having:

A ^he ^h＝R ^h (10)

The e here ^hIt mainly is the low frequency error of calculation under grid h.According to the analysis of front, only know to be difficult to obtain low frequency aberration e from equation (10) ^hThe ingenious part of many grid computings just is e ^hUnder original grid h, do not calculate, but on thicker grid, calculate.

For being of a size of H than coarse grid H _x* H _y, and H is arranged _x＞h _x, H _y＞h _y, on coarse grid, (10) can be written as:

A ^He ^H＝R ^H (11)

A wherein ^H, R ^HBe A ^h, R ^hThrough the new equation coefficient that obtains after the yardstick mapping, mapping relations are as follows:

A ^H＝P ^h→H(A ^h)，R ^H＝P ^h→H(R ^h)

P ^{H → H}Be mapping function, e ^HReflection be high frequency error of calculation grid H under, the while is again the low frequency error of calculation under the refined net h.Because e ^HAnd e ^hNot on a yardstick, therefore also need e ^HBe mapped as e ^h, as follows:

e ^h＝P ^H→h(e ^H)

High frequency error of calculation e under grid H ^HOnly need just can obtain, by mapping P by iteration seldom ^{H → h}(e ^H) finally can obtain e ^h, e has been arranged ^hJust can upgrade by (12)

${\stackrel{\&OverBar;}{U}}_{\mathrm{new}}^h=e^h+{\stackrel{\&OverBar;}{U}}^h---\left(12\right)$

Because the mapping between the yardstick also can bring certain high frequency error, therefore can also under grid h, carry out n more at last ₂Inferior iteration is eliminated

High frequency error.

Many networks can be combined with multiple dimensioned thought in the image calculation and be applied to PIV and calculate.The strategy that adopts multiple dimensioned (might scale factor less than 2) to approach step by step can reduce the influence that the equation model error is brought, and improves computational accuracy.Each dimensional variation all can carry out once upgrading based on the compensation of current estimated result.Regard the variation of graphical rule the variation of different size grid as, multiple dimensioned in essence approaching calculated a kind of special case that can be regarded as multi-grid algorithm.

Therefore, also can adopt the grid of more sizes according to the actual computation needs.Multi-grid algorithm is divided three classes in practice usually: the many grid computings of V-type, the many grid computings of W type and the many grid computings of perfect form, Fig. 4 a and Fig. 4 b have provided the synoptic diagram of V-type and the many grid computings of W type.

Preferably, employing many networks as shown in Figure 5.Round dot among Fig. 5 is represented the grid (size of mesh opening increases successively from bottom to up) of different scale, and the rising arrow represents from the small scale grid that to the transition of large scale grid arrow is then opposite downwards.Each just upgrades equation coefficient A and F in the formula (7) after finishing the many grid computings of V-type, carries out repeatedly circular treatment.The speed of convergence of system of linear equations iteration can be effectively accelerated by many grid computings, computation's reliability and precision can be do not reduced simultaneously.

Above content be in conjunction with concrete preferred implementation to further describing that the present invention did, can not assert that concrete enforcement of the present invention is confined to these explanations.For the general technical staff of the technical field of the invention, without departing from the inventive concept of the premise, can also make some simple deduction or replace, all should be considered as belonging to protection scope of the present invention.

Claims (3)

1. multi-grid processing method that is used for particle image velocimetry, described particle image velocimetry comprises finds the solution system of linear equations AU=F to realize the step to the optimal approximation of sports ground, U finds the solution for waiting, A is a matrix of coefficients, it comprises the motion vector field parameter at least, F is a vector, it is characterized in that, in the image processing process of particle image velocimetry, for graphical analysis, the definition of grid adopts each pixel of original image to represent refined net, based on the value computing on such discrete grid block point, the multiple dimensioned strategy that approaches is step by step adopted in the image motion field, each dimensional variation all can carry out once upgrading based on the compensation of current estimated result, the variation of graphical rule regarded as said method comprising the steps of the variation of different size grid:

A, first grid of image is set, obtains the system of equations A of system of equations AU=F correspondence under the first grid h ^hU ^h=F ^h

B, iterative system of equations A ^hU ^h=F ^h, separate according to gained is actual

Obtain error in equation $R^h=F^h-A^h{\stackrel{\&OverBar;}{U}}^h;$

C, the size of images second grid H thicker than first grid is set, will be to system of equations A under the first grid h ^he ^h=R ^hThe yardstick of finding the solution be mapped as under the second grid H system of equations A ^He ^H=R ^HFind the solution;

D, the iterative high frequency error of calculation e under the second grid H ^H

E, to the high frequency error of calculation e under the second grid H ^HCarry out the yardstick mapping, obtain the low frequency error of calculation e under the first grid h ^h

F, with the low frequency error of calculation e under the first grid h ^hUpgrade actual separating:

2. multi-grid processing method as claimed in claim 1 is characterized in that, also comprises:

After upgrading actual separating, iterative A under the first grid h of image again ^hU ^h=F ^h, eliminate

The high frequency error of calculation.

3. multi-grid processing method as claimed in claim 1 or 2 is characterized in that, the processing of described steps A to step F repeatedly carried out in circulation, and wherein the yardstick of second grid of image increases progressively one by one.

Images