CN103793924A - Flow field image self-adaption motion vector estimating method based on FHT-CC - Google Patents

Abstract

The invention discloses a flow field image self-adaption motion vector estimating method based on FHT-CC. Firstly, a flow field image is obtained, another flow field image is obtained after a certain time period and signals of the two flow field images are converted into Fourier frequency domain manners; secondly, the frequency domain correlation measure of the two flow filed images are calculated by means of two-dimension Hartley conversion cross-correlation and an air domain correlation curved face is obtained through inverse transformation; thirdly, a Gaussian curve equation is used for fitting sub-pixel coordinates of the peak value of the air domain correlation curve face which is obtained in the first step; finally, a window scale research zone is determined. In the zone, the self-adaption window selecting technology is adopted for searching for local peak values of the SNRs of the correlation curve faces to obtain the optimal window scale. Under the optimal window scale, the coordinates of the peak value of the correlation curve faces obtained in the second step serve as motion vectors of the flow field images. According to the flow field image self-adaption motion vector estimating method based on the FHT-CC, the FHT-CC serves as the correlation measure, the equivalence and the completeness of the correlation curve faces are maintained and at the same time, the calculated amount is greatly downsized.

Description

A kind of flow field figure based on FHT-CC is as adaptive motion vector method of estimation

Technical field

The present invention relates to a kind of flow field figure based on FHT-CC as adaptive motion vector method of estimation, belong to contactless instantaneous whole audience fluid-velocity survey technical field.

Background technology

Between two more than ten years in the past, the a collection of outstanding scientist such as Adrian and Merzkirch, realization and application in this contactless instantaneous whole audience fluid-velocity survey technology of particle image velocimetry (PIV) conduct in-depth research, make particle image velocimetry become a kind of applicable technology from principle, greatly improved the measurement capability of various Complex Flows under laboratory environment.Up to now, in two-dimentional whole-field velocity measurement technology, PIV is the most ripe a kind of new technology, becomes rapidly the standard method of testing the speed, and its product has also moved towards market (Technical Sourcing Internation of the U.S., Aerometrics company and Dantec company of Denmark etc. all have complete sets of products to release).

It is the developing milestone of experimental fluid mechanics now that particle is followed the tracks of speed-measuring method, it records the nowed forming in flow field by the method for shooting, adopt image analysis technology to calculate the velocity flow profile in whole flow field, be a kind of whole-field velocity measurement technology of instantaneous, unperturbed, in the Complex Flows such as research eddy current and turbulent flow, be widely used in recent years.Wherein the motion vector of particle picture estimating of fluid is core and the difficult point of PIV technology.

But for utilizing the natural water surface model such as the natural floating things such as leaf and whirlpool rolls, surface wave as for the large scale particle image velocimetry (LSPIV) of current tracer, most of natural water surface models result from Turbulence in Open channel Flow, they do not have clear and definite particle properties, and conventionally there is vibration, diffusion and distorted movement, thereby lack stable geometric properties, make particle follow the tracks of speed-measuring method and be difficult to directly be suitable for.Directly the spatial domain correlation matching algorithm of simple crosscorrelation (DCC) is computing method conventional in the motion vector of fluid, but for LSPIV, flow field area to be measured often covers hundreds of to thousands of square metres, image resolution ratio need reach 1,000,000 to ten million pixels, much larger than the order of magnitude of current laboratory 100,000 pixels, make the spatial domain correlation matching algorithm calculated amount based on DCC very large, be difficult to meet flow field, the flow requirement of real-time of monitoring continuously.Spectrum correlation matching method, based on fast fourier transform simple crosscorrelation (FFT-CC), in the time being input as real signal, is output as complex signal, and the calculating of imaginary part will consume a large amount of time and storage resources.In addition, because spectrum correlation coupling adopts fixing observation window, the tracer of the analyzed area of not considering to come in and go out, window size is comparatively responsive on the impact of measuring accuracy.At present the self-adaptation in PIV selects window method using related coefficient as iterated conditional, requires population in analyzed area more than 8～15, and the guarantee effectively related coefficient of coupling keeps stable, and this is conventionally difficult to meet in river Surface Picture.Therefore, the estimation of motion vectors in water surface flow field, river is faced with distinctive challenge, and research has important theory significance and using value for the method for estimating motion vector of this class non-rigid object of water surface pattern.

Summary of the invention

Technical matters to be solved by this invention is the deficiency for above-mentioned background technology, provide a kind of flow field figure based on FHT-CC as adaptive motion vector method of estimation, adopt FHT-CC as correlated measure, in keeping correlation surface equivalence and integrality, significantly simplify calculated amount.

In order to realize above-mentioned technical purpose, technical scheme of the present invention is:

Flow field figure based on FHT-CC is as an adaptive motion vector method of estimation, and the method comprises the following steps:

Step 1: first obtain a flow field figure picture, every after a time span, obtain another flow field figure picture, the signal of two flow field figure pictures is converted to Fourier's frequency domain form;

Step 2: the spectrum correlation that utilizes two-dimentional hartley conversion simple crosscorrelation model to calculate above-mentioned two width flow field figure pictures is estimated, then obtains spatial domain correlation surface through inverse transformation;

Step 3: adopt Gaussian curve equation model to go out the subpixel coordinates of above-mentioned spatial domain correlation surface peak value;

Step 4: determine the window size region of search of calculating above-mentioned spatial domain correlation surface according to sampling thheorem, in this interval, adopt self-adaptation to select the local peaking of window setting technique search correlation surface signal to noise ratio (S/N ratio) to obtain best window yardstick, and the subpixel coordinates of the correlation surface peak value that step 3 under best window yardstick is tried to achieve is as the motion vector of flow field figure picture.

Two-dimentional hartley conversion simple crosscorrelation (FHT-CC) model described in step 2:

C _H(u,v)＝F _He(u,v)G _He(u,v)+F _He(u,v)G _He(u,v)+F _He(u,-v)G _Ho(u,v)-F _Ho(u,v)G _He(u,-v)，

Wherein, F _h(u, v) and G _h(u, v) respectively corresponding F (u, v) converts with the two-dimentional hartley of G (u, v), F _he(u, v), F _he(u ,-v) be the odd function of former frame image through hartley transition structure, G _he(u, v), G _he(u ,-v) be the odd function of a rear two field picture through hartley transition structure, F _ho(u, v) is the even function of former frame image through hartley transition structure, G _ho(u, v) is the even function of a rear two field picture through hartley transition structure, C _h(u, v) estimates for spectrum correlation, and F (u, v) is Fourier's frequency domain form of former frame image, and G (u, v) is Fourier's frequency domain form of a rear two field picture.

Described in step 3, the method for asking of subpixel coordinates is:

Step 2 gained spatial domain correlation surface peak value distributes and can be expressed as at the one dimension of x direction:

$C\left(x\right)=\mathrm{Aexp}[-{(x-x_p)}^2/2{\mathrm{\σ}}_x^2],$

Wherein, x _pfor spatial domain correlation surface peak value is at the subpixel coordinates of x direction, A is the amplitude of correlation surface peak value, σ _xfor the standard deviation of x direction, set up a local coordinate system near of correlation surface peak value, obtain:

$x_p=x_m+\frac{\ln C(x_m-1)-\ln C(x_m+1)}{2\ln C(x_m-1)-4\ln C\left(x_m\right)+2\ln C(x_m+1)},$

$y_p=y_m+\frac{\ln C(y_m-1)-\ln C(y_m+1)}{2\ln C(y_m-1)-4\ln C\left(y_m\right)+2\ln C(y_m+1)},$

Wherein, note x _mand y _mfor whole pixel coordinate corresponding to maximum correlation coefficient, obtain thus the coordinate (x of the correlation surface peak value of sub-pixel precision _p, y _p).

Described in step 4, the method for asking of the window size region of search is:

For a fluid micellar in rectangular area, if adopt the stationary window that yardstick is M × N to observe its average displacement Δ x and Δ y in the x and y direction, the window size region of search is: 2 Δ x≤M≤3 Δ x, 2 Δ y≤N≤3 Δ y.

Self-adaptation described in step 4 selects window setting technique to comprise the following steps:

(1) choose initial gauges: choose one and be greater than in analyzed area in x and y direction the maximum square observation window for the treatment of 3 times of displacements as initial gauges:

M ₀＝N ₀≥3max(Δx,Δy)；

(2) calculate initial displacement: first with initial gauges M ₀× N ₀calculate FHT-CC, obtain observation window internal object x and y direction initial displacement Δ x ₀, Δ y ₀; Then adopt overall angular histogram to carry out verification to the displacement vector obtaining, and marking error vector is for flow field aftertreatment;

(3) iterative window yardstick: for correct vector, using 3 times of even numbers of displacement before this as current window size:

M _i＝3x _i-1,N _i＝3y _i-1,i＝1,2,3…

Recalculate FHT-CC; If gained vector is error vector, illustrate that current window size is not enough to the target information that provides enough, return to correct window size M corresponding to vector before this _i-1, N _i-1; If gained vector is correct vector, in the time that 3 times of displacement are greater than current window size, repeating step 3, the M that iteration makes new advances _i, N _i, until meet:

3Δx _i≤M _i,3Δy _i≤N _i；

(4) search SNR peak value: will approach the SNR peak value of 3 times of displacements most as the criterion of best window yardstick, search for successively the SNR peak value of large sense of displacement (x direction) and little sense of displacement (y direction) take 2pixel as the stepping of successively decreasing, until meet SNR _i>=SNR _i-1and SNR _i>=SNR _i+1in time, stops search, and returns to the best window yardstick M of analyzed area _i, N _iwith displacement x _i, Δ y _i; If there is not SNR peak value, in the time that window size is less than 2 times of initial displacement, i.e. M _i< 2 Δ x ₀, N _i>=2 Δ y ₀, stop search equally, return to window size M corresponding to maximum S/N R before this _i-1, N _i-1with displacement x _i-1, Δ y _i-1.

The beneficial effect that adopts technique scheme to bring is:

(1) there is the deficiency of redundant computation for FFT-CC in the method for estimating motion vector in the water surface flow field, river that the present invention proposes, adopt FHT-CC as correlated measure, in keeping correlation surface equivalence and integrality, calculated amount has reduced nearly half, and this advantage is along with the increase of window size is more obvious.

(2) the present invention is directed to peak value lockout issue, adopt 3 gaussian curve approximations to go out the subpixel coordinates of correlation peak, improved Displacement Estimation precision.

(3) the present invention is directed to the On The Choice of window size, best window yardstick district while having derived based on sampling thheorem that tracer meets followability and distribution occasion, adopt by slightly searching for the SNR peak value close to 3 times of displacements to smart Policy iteration, and then obtain adaptively best window yardstick, effectively improve the spatial resolution of flow field survey.

Accompanying drawing explanation

Fig. 1 is schematic flow sheet of the present invention.

Embodiment

Below with reference to accompanying drawing, technical scheme of the present invention is elaborated.

If Fig. 1 is schematic flow sheet of the present invention, a kind of flow field figure based on FHT-CC is as adaptive motion vector method of estimation, and the method comprises the following steps:

Step 1: first obtain a flow field figure picture, every after a time span, obtain another flow field figure picture, the signal of two flow field figure pictures is converted to Fourier's frequency domain form, the time span between former frame image and a rear two field picture is determined as the case may be.

Step 2: the spectrum correlation that utilizes two-dimentional hartley conversion simple crosscorrelation (FHT-CC) to calculate above-mentioned two width flow field figure pictures is estimated, and obtains spatial domain correlation surface through inverse transformation.First Fourier's frequency domain form F (u, v) of former frame image is carried out to hartley conversion and obtain F _h(u, v), carries out hartley conversion by Fourier's frequency domain form G (u, v) of a rear two field picture and obtains G _h(u, v), recycling formula (1) is tried to achieve spectrum correlation and is estimated C _h(u, v), to C _h(u, v) carries out inverse transformation and obtains spatial domain correlation surface C (d _x, d _y).

Wherein, two-dimentional hartley conversion simple crosscorrelation (FHT-CC) model is:

C _H(u,v)＝F _He(u,v)G _He(u,v)+F _He(u,v)G _He(u,v)+F _He(u,-v)G _Ho(u,v)-F _Ho(u,v)G _He(u,-v) （1）

In formula (1), F _he(u, v), F _he(u ,-v) be the odd function of former frame image through hartley transition structure, G _he(u, v), G _he(u ,-v) be the odd function of a rear two field picture through hartley transition structure, F _ho(u, v) is the even function of former frame image through hartley transition structure, G _ho(u, v) is the even function of a rear two field picture through hartley transition structure.

Step 3: adopt Gaussian curve equation model to go out the subpixel coordinates of above-mentioned gained spatial domain correlation surface peak value.

Step 2 gained spatial domain correlation surface peak value distributes and can be expressed as at the one dimension of x direction:

$C\left(x\right)=\mathrm{Aexp}[-{(x-x_p)}^2/2{\mathrm{\&sigma;}}_x^2]---\left(2\right)$

In formula (2), x _pfor spatial domain correlation surface peak value is at the subpixel coordinates of x direction, A is the amplitude of correlation surface peak value, σ _xfor the standard deviation of x direction, set up a local coordinate system near of correlation surface peak value, obtain:

$x_p=x_m+\frac{\ln C(x_m-1)-\ln C(x_m+1)}{2\ln C(x_m-1)-4\ln C\left(x_m\right)+2\ln C(x_m+1)}---\left(3\right)$

$y_p=y_m+\frac{\ln C(y_m-1)-\ln C(y_m+1)}{2\ln C(y_m-1)-4\ln C\left(y_m\right)+2\ln C(y_m+1)}---\left(4\right)$

In above formula, note x _mand y _mfor whole pixel coordinate corresponding to maximum correlation coefficient, obtain thus the coordinate (x of the correlation surface peak value of sub-pixel precision _p, y _p).

Step 4: determine the window size region of search of calculating correlation surface according to sampling thheorem, in this interval, adopt self-adaptation to select the local peaking of window setting technique search correlation surface signal to noise ratio (S/N ratio) to obtain best window yardstick, and the correlation surface peak coordinate that step 2 under best window yardstick is tried to achieve is as the motion vector of flow field figure picture.

The method of asking of the above-mentioned window size region of search is:

For a fluid micellar in rectangular area, observe it at x and y average displacement Δ x and Δ y if adopt the stationary window that yardstick is M × N, according to the sampling thheorem in information theory, the yardstick of observation window should be greater than fluid micellar 2 times of maximum displacement on correspondence direction:

M≥2ΔX _m,N≥2ΔY _m （5）

In formula (5), Δ X _m, Δ Y _mbe respectively fluid micellar maximum displacement in the x and y direction, in the time that the displacement of tracer meets consistance distributional assumption, have Δ X _m=Δ x, Δ Y _m=Δ _y, therefore the range scale of observation window is: M>=2 Δ x, N>=2 Δ y;

In the time that tracer meets the hypothesis being evenly distributed, between maximum displacement and average displacement, meet following relation:

$\mathrm{\&Delta;}{S}_m\&le;\sqrt{\mathrm{\&Delta;}{x}^2+\mathrm{\&Delta;}{y}^2}---\left(6\right)$

In formula (6), Δ S _mfor the tracer at center, rectangular area moves to the distance of angle point, for the situation of Δ x>=Δ y and Δ y>=Δ x, have respectively:

$\mathrm{\&Delta;}{X}_m=\mathrm{\&Delta;}{S}_m\&le;\sqrt{2}\mathrm{\&Delta;x},\mathrm{\&Delta;}{Y}_m=\mathrm{\&Delta;}{S}_m\&le;\sqrt{2}\mathrm{\&Delta;y}---\left(7\right)$

Can obtain in observation window best scale and be limited to:

this is that the tracer of rectangular area corner point moves ultimate range Δ S _mafter scope, obviously the tracer of other positions can not move out more than this scope composite type, the interval that can obtain best window yardstick place is:

2Δx≤M≤3Δx,2Δy≤N≤3Δy （8）

Above-mentioned self-adaptation selects window setting technique to comprise the following steps:

1, choose initial gauges: choose one and be greater than in analyzed area in x and y direction the maximum square observation window for the treatment of 3 times of displacements as initial gauges M ₀× N ₀:

M ₀＝N ₀≥3max(Δx,Δy) （9）

2, calculate initial displacement: first with initial gauges M ₀× N ₀calculate FHT-CC, obtain observation window internal object x and y direction initial displacement Δ x ₀, Δ y ₀; Then adopt overall angular histogram to carry out verification to the displacement vector obtaining, and marking error vector is for flow field aftertreatment;

3, iterative window yardstick: for correct vector, using 3 times of even numbers of displacement before this as current window size:

M _i＝3x _i-1,N _i＝3y _i-1,i＝1,2,3… （10）

Recalculate FHT-CC; If gained vector is error vector, illustrate that current window size is not enough to the target information that provides enough, return to correct window size M corresponding to vector before this _i-1, N _i-1; If gained vector is correct vector, in the time that 3 times of displacement are greater than current window size, repeating step 3, the M that iteration makes new advances _i, N _i, until meet:

3Δx _i≤M _i,3Δy _i≤N _i （11）

4, search SNR peak value: will approach the SNR peak value of 3 times of displacements most as the criterion of best window yardstick, search for successively the SNR peak value of x direction and y direction take 2pixel as the stepping of successively decreasing, until meet SNR _i>=SNR _i-1and SNR _i>=SNR _i+1in time, stops search, and returns to the best window yardstick M of analyzed area _i, N _iwith displacement x _i, Δ y _i; If there is not SNR peak value, in the time that window size is less than 2 times of initial displacement, i.e. M _i< 2 Δ x ₀, N _i>=2 Δ y ₀, stop search equally, return to window size M corresponding to maximum S/N R before this _i-1, N _i-1with displacement x _i-1, Δ y _i-1.

Above embodiment only, for explanation technological thought of the present invention, can not limit protection scope of the present invention with this, every technological thought proposing according to the present invention, and any change of doing on technical scheme basis, within all falling into protection domain of the present invention.

Claims (5)

1. the flow field figure based on FHT-CC, as an adaptive motion vector method of estimation, is characterized in that, the method comprises the following steps:

Step 1: first obtain a flow field figure picture, every after a time span, obtain another flow field figure picture, the signal of two flow field figure pictures is converted to Fourier's frequency domain form;

Step 2: the spectrum correlation that utilizes two-dimentional hartley conversion simple crosscorrelation model to calculate above-mentioned two width flow field figure pictures is estimated, then obtains spatial domain correlation surface through inverse transformation;

Step 3: adopt Gaussian curve equation model to go out the subpixel coordinates of above-mentioned spatial domain correlation surface peak value;

Step 4: determine the window size region of search of calculating above-mentioned spatial domain correlation surface according to sampling thheorem, in this interval, adopt self-adaptation to select the local peaking of window setting technique search correlation surface signal to noise ratio (S/N ratio) to obtain best window yardstick, and the subpixel coordinates of the correlation surface peak value that step 3 under best window yardstick is tried to achieve is as the motion vector of flow field figure picture.

2. a kind of flow field figure based on FHT-CC, as adaptive motion vector method of estimation, is characterized in that according to claim 1, and two-dimentional hartley conversion simple crosscorrelation model is described in step 2:

C _H(u,v)＝F _He(u,v)G _He(u,v)+F _He(u,v)G _He(u,v)+F _He(u,-v)G _Ho(u,v)-F _Ho(u,v)G _He(u,-v)，

Wherein, F _h(u, v) and G _h(u, v) respectively corresponding F (u, v) converts with the two-dimentional hartley of G (u, v), F _he(u, v), F _he(u ,-v) be the odd function of former frame image through hartley transition structure, G _he(u, v), G _he(u ,-v) be the odd function of a rear two field picture through hartley transition structure, F _ho(u, v) is the even function of former frame image through hartley transition structure, G _ho(u, v) is the even function of a rear two field picture through hartley transition structure, C _h(u, v) estimates for spectrum correlation, and F (u, v) is Fourier's frequency domain form of former frame image, and G (u, v) is Fourier's frequency domain form of a rear two field picture.

3. a kind of flow field figure based on FHT-CC, as adaptive motion vector method of estimation, is characterized in that according to claim 1, and the method for asking of subpixel coordinates is described in step 3:

Step 2 gained spatial domain correlation surface peak value distributes and can be expressed as at the one dimension of x direction:

$C\left(x\right)=\mathrm{Aexp}[-{(x-x_p)}^2/2{\mathrm{\&sigma;}}_x^2],$

Wherein, x _pfor spatial domain correlation surface peak value is at the subpixel coordinates of x direction, A is the amplitude of correlation surface peak value, σ _xfor the standard deviation of x direction, set up a local coordinate system near of correlation surface peak value, obtain:

$x_p=x_m+\frac{\ln C(x_m-1)-\ln C(x_m+1)}{2\ln C(x_m-1)-4\ln C\left(x_m\right)+2\ln C(x_m+1)},$

$y_p=y_m+\frac{\ln C(y_m-1)-\ln C(y_m+1)}{2\ln C(y_m-1)-4\ln C\left(y_m\right)+2\ln C(y_m+1)},$

Wherein, note x _mand y _mfor whole pixel coordinate corresponding to maximum correlation coefficient, obtain thus the coordinate (x of the correlation surface peak value of sub-pixel precision _p, y _p).

4. a kind of flow field figure based on FHT-CC, as adaptive motion vector method of estimation, is characterized in that according to claim 1, and the window size region of search is described in step 4:

For a fluid micellar in rectangular area, if adopt the stationary window that yardstick is M × N to observe its average displacement Δ x and Δ y in the x and y direction, the window size region of search is: 2 Δ x≤M≤3 Δ x, 2 Δ y≤N≤3 Δ y.

5. a kind of flow field figure based on FHT-CC, as adaptive motion vector method of estimation, is characterized in that according to claim 1, and self-adaptation selects window setting technique to comprise the following steps described in step 4:

(1) choose initial gauges: choose one and be greater than in analyzed area in x and y direction the maximum square observation window for the treatment of 3 times of displacements as initial gauges:

M ₀＝N ₀≥3max(Δx,Δy)；

(2) calculate initial displacement: first with initial gauges M ₀× N ₀calculate FHT-CC, obtain observation window internal object x and y direction initial displacement Δ x ₀, Δ y ₀, then adopt overall angular histogram to carry out verification to the displacement vector obtaining, and marking error vector is for flow field aftertreatment;

(3) iterative window yardstick: for correct vector, using 3 times of even numbers of displacement before this as current window size:

M _i＝3x _i-1,N _i＝3y _i-1,i＝1,2,3…

Recalculate FHT-CC; If gained vector is error vector, illustrate that current window size is not enough to the target information that provides enough, return to correct window size M corresponding to vector before this _i-1, N _i-1; If gained vector is correct vector, in the time that 3 times of displacement are greater than current window size, repeating step 3, the M that iteration makes new advances _i, N _i, until meet:

3Δx _i≤M _i,3Δy _i≤N _i；

(4) search SNR peak value: will approach the SNR peak value of 3 times of displacements most as the criterion of best window yardstick, search for successively the SNR peak value of x direction and y direction take 2pixel as the stepping of successively decreasing, until meet SNR _i>=SNR _i-1and SNR _i>=SNR _i+1in time, stops search, and returns to the best window yardstick M of analyzed area _i, N _iwith displacement x _i, Δ y _i; If there is not SNR peak value, in the time that window size is less than 2 times of initial displacement, i.e. M _i< 2 Δ x ₀, N _i>=2 Δ y ₀, stop search equally, return to window size M corresponding to maximum S/N R before this _i-1, N _i-1with displacement x _i-1, Δ y _i-1.

Images