CN104657955A - Displacement field iteration smoothing method of kernel function based digital image correlation method - Google Patents

Abstract

The invention discloses a displacement field iteration smoothing method of a kernel function based digital image correlation method. A kernel function based correlation function rho* is adopted, and a displacement field obtained with the kernel function based digital image correlation method is processed with the iteration smoothing method. According to the displacement field iteration smoothing method of the kernel function based digital image correlation method, the measurement accuracy of the displacement field can be effectively improved, one new kernel function based correlation function rho* is defined, and image noise influence is considered, so that the displacement field iteration smoothing method has higher measurement accuracy compared with a traditional method; according to the displacement field iteration smoothing method, a rough penalty term is added on the basis of an original least square method, so that the accuracy and the roughness of a solution are balanced, a smoothing parameter alpha can be self-adaptively solved and does not require manual selection, convenience and practicability are achieved, output smooth displacement field is computed, and a result is more accurate.

Description

Based on the displacement field iteration smoothing method of the Digital Image Correlation Method of kernel function

Technical field

The invention belongs to digital image processing field, be specifically related to a kind of displacement field iteration smoothing method of the Digital Image Correlation Method based on kernel function.

Background technology

Digital picture is correlated with, proposed by Sutton et al. the earliest, be a kind of using the natural texture of body surface stochastic distribution or artificial speckle field as the carrier of deformation information, obtain the whole audience displacement of body structure surface under external load function and the measuring method of strain information by machine vision technique.Because DIC has plurality of advantages, such as measuring process is simple, measurement result accuracy is high, noncontact, can obtain full field data etc., is applied widely in recent years in Experimental Mechanics field.

As the effective measuring method of one, DIC technology obtains and develops rapidly in recent two decades, and at test mechanics, material and association area thereof obtain a large amount of successful Application.These are all owing to all advantages of DIC method itself, and such as contactless, pilot system is built facility, can be obtained full field data and relatively high measurement accuracy, etc.For this reason, large metering method is proposed for and improves the measuring accuracy of DIC and processing speed, promotes the method and to succeed more in a wider context application.Unfortunately, although obtain sufficient research in recent years based on the measure theory of DIC, still there is great number of issues in classic method, when in the face of practical application.The factors such as speckle pattern, image subsection size, related function, iterated conditional can to measurement result noise effect.For problems, more existing scholars carried out comparatively careful analysis to DIC performance, and proposed the bulking property principle that some instruct practical application.

Must be pointed out, in real world applications, picture noise is inevitable.And the measuring process of traditional DIC method depends critically upon the monochrome information of image.When noise exists, also can change along with the distribution of noise for the Luminance Distribution in the image subsection that mates before and after distortion, thus drastically influence the accuracy of measurement result.For reducing picture noise to the impact of measurement result, some methods eliminating measured deviation based on pre-treatment or post-processing technology are suggested.Although obtain certain effect, when being used for actual measurement, these methods still face the problem lacking smoothing parameter and how to choose.

The ultimate principle of DIC is very simple, the sequence image of interrecord structure distortion, and before being deformed after image subsection in search maximize certain correlation criterion, such as: zero-mean normalization crosscorrelation criterion, obtains the displacement data of measurement point.In order to improve the computational accuracy of displacement field, research is mainly paid close attention to and is exported high-precision Displacement by improving DIC algorithm.Higher-order gradients is introduced shape function to realize the description to complex deformation by Lu et al., thus improves the measurement effect of DIC to complex deformation field.Cofaru et al. proposes to revise DIC by irregular speckle pattern, and increases the output accuracy of displacement field in conjunction with regularization method.The research of Pan shows: employing size is that the gauss low frequency filter of 5*5 pixel carries out pre-service to speckle image, effectively can reduce the displacement field error that DIC exports.

In Experimental Mechanics field, compare simple displacement field data, strain field distribution information seems more valuable.Although DIC technology experienced by development for many years, due to current conditions restriction, still can there is various deviation, about the analysis of DIC system deviation refers to document in the displacement field that DIC exports.Because all information calculating strain field are included in the middle of displacement data, if directly utilize these displacement datas comprising noise to calculate strain, error will be exaggerated, to such an extent as to effective strain field distribution rule is difficult to obtain.Just because of this, the post-processing technology based on data smoothing or surface fitting is used to the noise eliminating displacement field.But the defect of these methods is the parameters needing artificial adjustment algorithm.In actual measurement process, owing to not having enough prioris to carry out guide parameters adjustment, thus hinder its practical application.

Therefore, need a kind of new Digital Image Correlation Method displacement field smoothing method to solve the problem.

Summary of the invention

The object of the invention is the deficiency for Digital Image Correlation Method in prior art and displacement field smoothing method thereof, a kind of displacement field iteration smoothing method of the Digital Image Correlation Method based on kernel function is provided.

For achieving the above object, the displacement field iteration smoothing method that the present invention is based on the Digital Image Correlation Method of kernel function can adopt following technical scheme:

Based on a displacement field iteration smoothing method for the Digital Image Correlation Method of kernel function, the Digital Image Correlation Method based on kernel function adopts the related function ρ based on kernel function ^*,

${\mathrm{\ρ}}^{*}=c\underset{s_p\∈S_0}{\mathrm{\Σ}}k\left({\left|\right|\frac{f\left(s_p\right)-g(s_p,P)}{h}\left|\right|}^2\right)$

Wherein, c is normaliztion constant, and k (.) is kernel function, and h is the bandwidth control parameter of kernel function, s _pfor being out of shape front image-region S ₀in pixel, f (s _p) be pixel s in image before distortion _pthe brightness of image at place, g (s _p, P) for after distortion in image with pixel s _pthe brightness of image at corresponding pixel place;

Wherein, iteration smoothing method comprises the following steps:

(1), the displacement field data d after distortion is measured based on the Digital Image Correlation Method of kernel function described in utilization ⁿ, wherein displacement field data d ⁿbe expressed from the next: d ⁿ=d+n;

Wherein, d be level and smooth after displacement field, n is displacement field noise;

(2), quadratic function is built $D(d,\mathrm{\α})=\mathrm{\α}{\left|\right|d-d^n\left|\right|}_2^2+(1-\mathrm{\α}){\left|\right|\mathrm{Cd}\left|\right|}_2^2$

Wherein, d ⁿfor the displacement field data after distortion, d be level and smooth after displacement field data, n is displacement field noise, and α is smoothing parameter, and C is high-order operator;

Solve the first-order partial derivative of quadratic function, obtain

$\begin{array}{c}\frac{\∂D(d,\mathrm{\α})}{\∂d}=2\mathrm{\α}(d-d^n)+2(1-\mathrm{\α})*4\left(\mathrm{Cd}\right),\\ \frac{\∂D(d,\mathrm{\α})}{\∂\mathrm{\α}}={\left|\right|d-d^n\left|\right|}_2^2-{\left|\right|\mathrm{Cd}\left|\right|}_2^2\end{array}$

Wherein, C is high-order operator,

D and α is become column vector and makes Q=[d1, d2, d3...dm, α], wherein, m is for counting; Then can obtain

$\▿D\left(Q\right)=\left[\begin{array}{c}2\mathrm{\α}(d_1-{d_1}^n)+2(1-\mathrm{\α})*4\left({\mathrm{Cd}}_1\right)\\ 2\mathrm{\α}(d_2-{d_2}^n)+2(1-\mathrm{\α})*4\left({\mathrm{Cd}}_2\right)\\ ......\\ 2\mathrm{\α}(d_m-{d_m}^n)+2(1-\mathrm{\α})*4\left({\mathrm{Cd}}_m\right)\\ {\left|\right|d-d^n\left|\right|}_2^2-{\left|\right|\mathrm{Cd}\left|\right|}_2^2\end{array}\right];$

(3), second-order partial differential coefficient ▽ ▽ D (Q) of quadratic function is solved according to the ▽ D (Q) of step (2),

$\frac{{\∂D}^2\left(Q\right)}{{\∂d}_i{\∂d}_j}=\left[\begin{array}{cccc}32-30\mathrm{\α}& & & \\ & 32-30\mathrm{\α}& & \\ & & ......& \\ & & & 32-30\mathrm{\α}\end{array}\right]$

$\frac{{\∂D}^2\left(Q\right)}{{\∂d}_i\∂\mathrm{\α}}=2(d-d^n)-8\left(\mathrm{Cd}\right)$

$\frac{{\∂D}^2\left(Q\right)}{{\∂\mathrm{\α}\∂d}_i}=2(d-d^n)-8\left(\mathrm{Cd}\right)$

$\frac{{\∂D}^2\left(Q\right)}{\∂\mathrm{\α}\∂\mathrm{\α}}=0$

Above-mentioned matrix is combined, obtains ▽ ▽ D (Q);

(4), the ▽ D (Q) that utilizes step (2) and (3) to obtain and ▽ ▽ D (Q), iterative computation is carried out to following formula:

$P_{k+1}=P_k-\frac{\▿D\left(Q\right)}{\▿\▿D\left(Q\right)},$

Wherein, the initial value of Q is [d ⁿ, α] ^t, wherein, d ⁿfor the displacement field data after distortion, when the d obtained after above formula convergence is desirable real displacement value.

Further, the scope of described α is (0,1).

Further, step 4) the condition of convergence be wherein, ε≤0.1.

Further, the described Digital Image Correlation Method based on kernel function comprises the following steps:

One), shape function is defined: the arbitrfary point (x in setting reference picture ₀, y ₀) and neighborhood S around, (x, y) is the coordinate of any pixel in neighborhood S in reference picture, for the coordinate of pixel corresponding with pixel (x, y) in target image, there is one group of mapping relations χ and following formula set up:

$\begin{array}{c}\mathrm{\χ}(x,y)\→(\tilde{x},\tilde{y})\\ f(x,y)=g(\tilde{x},\tilde{y})\end{array}$

Wherein, f (x, y) represents the brightness of image at pixel (x, y) place, represent pixel the brightness of image at place, mapping relations χ is shape function; Wherein, shape function χ is expressed from the next:

$\begin{array}{c}\tilde{x}=x+u(x_0,y_0)+\frac{\∂u}{\∂x}{|}_{x_0,y_0}\left({x-x}_0\right)+\frac{\∂u}{\∂y}{|}_{x_0,y_0}(y-y_0)\\ \tilde{y}=y+v(x_0,y_0)+\frac{\∂v}{\∂x}{|}_{x_0,y_0}(x-x_0)+\frac{\∂v}{\∂y}{|}_{x_0,y_0}(y-y_0)\end{array}$

Wherein, u and v is respectively the displacement of being out of shape x and the y direction caused, (x ₀, y ₀) be the center position coordinates of neighborhood S, for being out of shape the single order displacement gradient produced in the x and y direction;

Two), by step one) shape function parametrization, and to represent with vectorial P,

$P=(u,v,\frac{\∂u}{\∂x},\frac{\∂v}{\∂x},\frac{\∂u}{\∂y},\frac{\∂v}{\∂y}),$

Definition related function ρ ^*,

${\mathrm{\ρ}}^{*}=c\underset{s_p\∈S}{\mathrm{\Σ}}k\left({\left|\right|\frac{f\left(s_p\right)-g(s_p,P)}{h}\left|\right|}^2\right)$

Wherein, c is normaliztion constant, and k (.) is kernel function, and h is the bandwidth control parameter of kernel function, s _pfor the pixel of in neighborhood S, f (s _p) be pixel s in image before distortion _pthe brightness of image at place, g (s _p, P) for after distortion in image with pixel s _pthe brightness of image at corresponding pixel place;

Three), calculate and make related function ρ ^*shape function parameter P time minimum.

Further, step 3) in calculate related function ρ by following formula ^*minimum value, structure iterative equations:

$P=P_0-\frac{{\▿\mathrm{\ρ}}^{*}\left(P_0\right)}{{\▿\▿\mathrm{\ρ}}^{*}\left(P_0\right)}$

In formula, P ₀for deformation parameter initial value, ▽ ρ ^*with ▽ ▽ ρ ^*related function ρ ^*first-order Gradient and Hessian matrix, wherein ▽ ρ ^*with ▽ ▽ ρ ^*formula as follows:

$\begin{array}{c}{\▿\mathrm{\ρ}}^{*}={\left(\frac{{\∂\mathrm{\ρ}}^{*}}{{\∂P}_i}\right)}_{i=1,...,6}\\ =\frac{-2}{\underset{s_p\∈S}{\mathrm{\Σ}}{f}^2\left(s_p\right)}{\left\{\underset{s_p\∈S}{\mathrm{\Σ}}\right[f\left(s_p\right)-g(s_p,P)\left]\frac{\∂g(s_p,P)}{{\∂P}_i}\right\}}_{i=1,...,6}\end{array}$

$\begin{array}{c}{\▿\▿\mathrm{\ρ}}^{*}={\left(\frac{{{\∂}^2\mathrm{\ρ}}^{*}}{{\∂P}_i{\∂P}_j}\right)}_{\begin{array}{c}i=1,...,6\\ j=1,...,6\end{array}}\\ =\frac{-2}{\underset{s_p\∈S}{\mathrm{\Σ}}{f}^2\left(s_p\right)}{\left\{\underset{s_p\∈S}{\mathrm{\Σ}}\frac{\∂g(s_p,P)}{{\∂P}_i}\frac{\∂g(s_p,P)}{{\∂P}_j}\right\}}_{\begin{array}{c}i=1,...,6\\ j=1,...,6\end{array}}\end{array};$

Utilize iterative equations calculate and make related function ρ ^*minimum time solution.

Further, the described Digital Image Correlation Method based on kernel function comprises the following steps:

1), gather image before deformation of body and after distortion, the image before distortion is reference picture, and the image after distortion is target image;

2), region-of-interest S is chosen on a reference ₁, region-of-interest S ₁for reference picture subarea;

3), on target image, region of search S is set up ₂, region of search S ₂for target image subarea, target image subarea comprises the reference picture subarea after distortion;

4) gray-scale value of the point in target image subarea and reference picture subarea, is obtained;

5), build shape function, determine the position relationship of corresponding point in reference picture subarea and target image subarea, wherein, (x ₁, y ₁) be the coordinate of pixel any in reference picture subarea, (x ₂, y ₂) in target image subarea with pixel (x ₁, y ₁) coordinate of corresponding pixel, there is one group of mapping relations χ and following formula set up:

χ(x ₁,y ₁)→(x ₂,y ₂)

f(x ₁,y ₁)＝g(x ₂,y ₂)

Wherein, f (x ₁, y ₁) represent pixel (x ₁, y ₁) brightness of image at place, g (x ₂, y ₂) represent pixel (x ₂, y ₂) brightness of image at place, mapping relations χ is shape function; Wherein, shape function χ is expressed from the next:

$\begin{array}{c}{x}_2=x_1+u^{\′}(x^{\′},y^{\′})+\frac{{\∂u}^{\′}}{{\∂x}_1}{|}_{x^{\′},y^{\′}}(x_1-x^{\′})+\frac{{\∂u}^{\′}}{{\∂y}_1}{|}_{x^{\′},y^{\′}}(y_1-y^{\′})\\ y_2=y_1+v^{\′}(x^{\′},y^{\′})+\frac{{\∂v}^{\′}}{{\∂x}_1}{|}_{x^{\′},y^{\′}}(x_1-x^{\′})+\frac{{\∂v}^{\′}}{{\∂y}_1}{|}_{x^{\′},y^{\′}}(y_1-y^{\′})\end{array}$

Wherein, u and v is respectively the displacement of being out of shape x and the y direction caused, and (x', y') is the center position coordinates in reference picture subarea, for distortion is at x ₁and y ₁the single order displacement gradient that direction produces;

6), by step 5) shape function parametrization, and to represent with vectorial P',

$P^{\′}=(u^{\′},v^{\′},\frac{\∂u^{\′}}{\∂x_1},\frac{\∂v^{\′}}{\∂x_1},\frac{\∂u^{\′}}{\∂y_1},\frac{\∂v^{\′}}{\∂y_1}),$

Definition related function ρ ^*,

${\mathrm{\ρ}}^{*\′}=c\underset{{s_p}^{\′}\∈S_1}{\mathrm{\Σ}}k\left({\left|\right|\frac{f\left({s_p}^{\′}\right)-g({s_p}^{\′},P^{\′})}{h}\left|\right|}^2\right)$

Wherein, c is normaliztion constant, and k (.) is kernel function, and h is the bandwidth control parameter of kernel function, s _p' be the pixel of in reference picture subarea, f (s _p') be pixel s in image before distortion _p' place brightness of image, g (s _p', P') for after distortion in image with pixel s _p' the brightness of image at corresponding pixel place;

7), setting vector initial value P ₀', and set maximum iteration time n,

8), by P ₀' bring iterative equations into:

${P_1}^{\′}={P_0}^{\′}-\frac{{\▿\mathrm{\ρ}}^{*\′}\left({P_0}^{\′}\right)}{{\▿\▿\mathrm{\ρ}}^{*\′}\left({P_0}^{\′}\right)}$

In formula, P ₀' be deformation parameter initial value, ▽ ρ ^*' and ▽ ▽ ρ ^*' be related function ρ ^*first-order Gradient and Hessian matrix, wherein ▽ ρ ^*' and ▽ ▽ ρ ^*' formula as follows:

${\▿\mathrm{\ρ}}^{*\′}=\frac{-2}{\underset{s_p\∈S_1}{\mathrm{\Σ}}{f}^2\left({s_p}^{\′}\right)}{\left\{\underset{{s_p}^{\′}\∈S_1}{\mathrm{\Σ}}\right[f\left({s_p}^{\′}\right)-g({s_p}^{\′},P^{\′})\left]\frac{\∂g({s_p}^{\′},P^{\′})}{{{\∂P}_i}^{\′}}\right\}}_{i=1,...,6}$

${\▿\▿\mathrm{\ρ}}^{*\′}=\frac{-2}{\underset{s_p\∈S_1}{\mathrm{\Σ}}{f}^2\left({s_p}^{\′}\right)}{\left\{\underset{{s_p}^{\′}\∈S_1}{\mathrm{\Σ}}\frac{\∂g({s_p}^{\′},P^{\′})}{{{\∂P}_i}^{\′}}\frac{\∂g({s_p}^{\′},P^{\′})}{{{\∂P}_j}^{\′}}\right\}}_{\begin{array}{c}i=1,...,6\\ j=1,...,6\end{array}};$

In formula, P _i' and P _j' be respectively i-th element and a jth element of vectorial P';

P is calculated according to iterative equations ₁';

9), according to following formula judge whether iterative equations restrains, if convergence, then preserve region-of-interest S ₁shift value u' and v'; If convergence, does not make P ₀'=P ₁', and repeat step 8) and 9); If iterations reaches maximum iteration time n, then finishing iteration;

10) shift value of region-of-interest, is exported.

Further, step 4) in obtain the gray-scale value of the point in target image subarea by carrying out bicubic spline interpolation to target image subarea, the gray-scale value of sub-pixel location in acquisition target image subarea.

Further, step 4) in obtain the gray-scale value of the point in reference picture subarea by carrying out bicubic spline interpolation to reference picture subarea, the gray-scale value of sub-pixel location in acquisition reference picture subarea.

Further, step 2) in choose region-of-interest S on a reference ₁, region-of-interest S ₁for reference picture subarea is obtained by following steps: divide the first grid respectively on a reference, the number of the first grid be M capable × N row, spacing in the length and Width of reference picture between the first grid is respectively l and w, sets up the region-of-interest S centered by a the first grid on a reference ₁, region-of-interest S ₁for reference picture subarea, length and the width in reference picture subarea are respectively L and W, wherein, and a=1,2,3 ..., M × N.Utilize and choose region-of-interest by the mode of grid division on a reference, simple and convenient, easily realize.

Further, step 3) on target image, set up region of search S ₂, region of search S ₂for target image subarea, reference picture subarea after target image subarea comprises distortion is obtained by following steps: on target image, divide the second grid, the number of the second grid be M capable × N row, the spacing in the length and Width of target image between the second grid is respectively l and w; Target image is set up the region of search S centered by b the second grid ₂, region of search S ₂for target image subarea, length and the width in target image subarea are respectively k*L and k*W, wherein, and k>1, b=a.The mode by grid division on target image is utilized to choose region of search, simple and convenient, easily realize.

Beneficial effect: the displacement field iteration smoothing method of the Digital Image Correlation Method based on kernel function of the present invention effectively can improve the measuring accuracy of displacement field.Define a kind of likeness coefficient based on kernel function newly, consider the impact of picture noise, make the present invention's comparatively classic method, there is better measuring accuracy.Displacement field iteration smoothing method is on the basis in former least square method, add coarse penalty term, a balance is reached between the roughening that the accuracy of solution is conciliate, wherein, smoothing parameter α can self-adaptation solve, and need not select this parameter artificially, convenient, practical, calculate the level and smooth displacement field exported, result is more accurate.

Accompanying drawing explanation

Fig. 1, be the process flow diagram of the displacement field iteration smoothing method of the Digital Image Correlation Method based on kernel function of the present invention;

Reference picture before Fig. 2, distortion;

Image (adding 4% noise) after Fig. 3, distortion;

Four groups of displacement error comparison diagrams that Fig. 4, NR method calculates;

Four groups of shift standards difference comparison diagrams that Fig. 5, NR method calculates;

Fig. 6, the four groups of displacement error comparison diagrams calculated based on the displacement field iteration smoothing method of the Digital Image Correlation Method of kernel function;

Fig. 7, the four groups of shift standards difference comparison diagrams calculated based on the displacement field iteration smoothing method of the Digital Image Correlation Method of kernel function;

Fig. 8, be test specimen figure;

Fig. 9, be the treatment effect figure of the displacement field iteration smoothing method of Digital Image Correlation Method based on kernel function;

Figure 10, be test specimen handling principle figure;

Figure 11 is the process flow diagram of displacement field iteration smoothing method of the present invention;

Figure 12 is under different α initial condition, α iteration situation of change;

Figure 13 is the comparison diagram of N-R iterative computation shift value, level and smooth rear shift value and real displacement value;

Figure 14 is the result utilizing conventional method process DIC to export data in inhomogeneous deformation situation;

Figure 15 is the result utilizing displacement field iteration smoothing method process DIC output data to obtain in inhomogeneous deformation situation;

Figure 16 is simulated speckle pattern;

Figure 17 is speckle pattern;

Figure 18 is speckle schematic diagram;

Figure 19 is the tensile strain field deformation measurement result that Multiple-Hole Specimen utilizes the process of conventional DIC method to obtain;

Figure 20 is the shear strain field deformation measurement that Multiple-Hole Specimen utilizes the process of conventional DIC method and obtains;

Figure 21 is the tensile strain field deformation measurement result that Multiple-Hole Specimen utilizes method process of the present invention to obtain;

Figure 22 is the shear strain field deformation measurement that Multiple-Hole Specimen utilizes method process of the present invention and obtains.

Embodiment

Below in conjunction with the drawings and specific embodiments, illustrate the present invention further, these embodiments should be understood only be not used in for illustration of the present invention and limit the scope of the invention, after having read the present invention, the amendment of those skilled in the art to the various equivalent form of value of the present invention has all fallen within the application's claims limited range.

Tradition DIC method

The ultimate principle of DIC method is very simple, and before being converted into test piece deformation by the measurement problem of malformation, the relevant matches problem of (reference picture), rear (target image) image solves.Therefore, be the deformation on description scheme surface, first need to define shape function.Assuming that arbitrfary point (x, y) in reference picture and the little neighborhood S of around thereof, there is one group of mapping relations χ and meet

$\mathrm{\χ}(x,y)\→(\tilde{x},\tilde{y}),f(x,y)=g(\tilde{x},\tilde{y})$

Wherein, f (x, y) represents the brightness of image at point (x, y) place, represent the rear point of distortion the brightness of image at place.

Map χ and be called as so-called shape function.If neighborhood S and deflection are enough little, shape function χ can be described by formula (1),

$\begin{array}{c}\tilde{x}=x+u(x_0,y_0)+\frac{\∂u}{\∂x}{|}_{x_0,y_0}\left({x-x}_0\right)+\frac{\∂u}{\∂y}{|}_{x_0,y_0}(y-y_0)\\ \tilde{y}=y+v(x_0,y_0)+\frac{\∂v}{\∂x}{|}_{x_0,y_0}(x-x_0)+\frac{\∂v}{\∂y}{|}_{x_0,y_0}(y-y_0)\end{array}---\left(1\right)$

Wherein, u and v is respectively the in-plane displacement being out of shape x and the y direction caused, (x ₀, y ₀) be the center position coordinates of region S.

Write shape function as vector form,

$P=(u,v,\frac{\∂u}{\∂x},\frac{\∂v}{\∂x},\frac{\∂u}{\∂y},\frac{\∂v}{\∂y}),$

And define related coefficient,

$\mathrm{\ρ}=\frac{\underset{s_p\∈S}{\mathrm{\Σ}}{[f\left(s_p\right)-g(s_p,P)]}^2}{\underset{s_p\∈S}{\mathrm{\Σ}}{f}^2\left(s_p\right)}---\left(2\right)$

Then obtain by means of nonlinear optimization method and make formula (2) minimized optimum solution, problem just achieves a solution.

As can be seen from formula (2), when related function minimalization, before and after distortion, the similarity of image subsection reaches maximal value.Now, displacement parameter u and v that parameter vector P comprises represents the optimum estimate to displacement after distortion, calculates in the same way, can obtain whole audience displacement to all measurement points.

Minimize ρ and have multiple method for solving, known to author, Newton-Raphson method because computational accuracy is higher adopt by numerous document, be namely constructed as follows iterative equations

$P=P_0-\frac{\▿\mathrm{\ρ}\left(P_0\right)}{\▿\▿\mathrm{\ρ}\left(P_0\right)}---\left(3\right)$

In formula, P ₀for deformation parameter initial value, ▽ ρ and ▽ ▽ ρ is First-order Gradient and the Hessian matrix of related function ρ, meets

$\begin{array}{c}\▿\mathrm{\ρ}={\left(\frac{\∂\mathrm{\ρ}}{{\∂P}_i}\right)}_{i=1,...,6}\\ =\frac{-2}{\underset{s_p\∈S}{\mathrm{\Σ}}{f}^2\left(s_p\right)}{\left\{\underset{s_p\∈S}{\mathrm{\Σ}}\right[f\left(s_p\right)-g(s_p,P)\left]\frac{\∂g(s_p,P)}{{\∂P}_i}\right\}}_{i=1,...,6}\end{array}---\left(4\right)$

With

$\begin{array}{c}\▿\▿\mathrm{\ρ}={\left(\frac{{\∂}^2\mathrm{\ρ}}{{\∂P}_i{\∂P}_j}\right)}_{\begin{array}{c}i=1,...,6\\ j=1,...,6\end{array}}\\ \≈\frac{2}{\underset{s_p\∈S}{\mathrm{\Σ}}{f}^2\left(s_p\right)}{\left\{\underset{s_p\∈S}{\mathrm{\Σ}}\frac{\∂g(s_p,P)}{{\∂P}_i}\frac{\∂g(s_p,P)}{{\∂P}_j}\right\}}_{\begin{array}{c}i=1,...,6\\ j=1,...,6\end{array}}\end{array}---\left(5\right)$

Carefully analyze formula (3-5) can find out, traditional DIC in the process minimizing ρ, f (s _p)-g (s _p, P) and whether two very large on iteration result impact, even can affect iteration and restrain.Ideally (do not have noise), once, the parameter of shape function all can along gradient direction to true strain parameter convergence one step for every iteration.When iteration meets end condition, Ying You trend towards zero.But, when picture noise exists, not only f (s _p)-g (s _p, P) and item is not equal to zero, item also can because of the impact of noise magnification distortion.Therefore, traditional DIC can reduce measurement effect because of the existence of noise when practical application.

DIC method of the present invention

For improving the adaptability of traditional DIC to noise, kernel function is incorporated into the definition of similarity function herein.Kernel function is normally defined certain radially symmetrical scalar function.Such as, gaussian kernel function, i.e. so-called radial basis function (Radial Basis Function, RBF), to be defined as in space any point x to a certain center x _cbetween the monotonic quantity of Euclidean distance, can be denoted as k (|| x-x _c||), form is

$k\left(\right||x-x_c|\left|\right)=\exp (-\frac{{\left|\right|x-x_c\left|\right|}^2}{{2\mathrm{\σ}}^2})$

Wherein, x _cfor kernel function center, σ is the width parameter of function, the radial effect scope of control function, and its effect is local, when x is away from x _ctime function value very little.

For utilizing the control ability of kernel function in radial extension, in given neighborhood S, by the brightness of image before and after distortion is done poor mode, the similarity measurement of image subsection before and after distortion is transformed to the feature space centered by zero, and represent with the form of error of sum square, construct the related coefficient based on kernel function

${\mathrm{\ρ}}^{*}=c\underset{s_p\∈S}{\mathrm{\Σ}}k\left({\left|\right|\frac{f\left(s_p\right)-g(s_p,P)}{h}\left|\right|}^2\right)---\left(6\right)$

Wherein, c is normaliztion constant, and k (.) is kernel function, and h is the bandwidth control parameter of kernel function.

Express for simplifying, order

r(s _p,p)＝f(s _p)-g(s _p,p)

Formula (4) is reduced to,

${\mathrm{\ρ}}^{*}=c\underset{s_p\∈s}{\mathrm{\Σ}}k\left({\left|\right|\frac{r(s_p,p)}{h}\left|\right|}^2\right)---\left(7\right)$

First-order Gradient and the Hessian matrix of ρ be,

$\begin{array}{c}\▿\mathrm{\ρ}\left(p\right)={\left(\frac{\∂\mathrm{\ρ}}{{\∂p}_i}\right)}_{i=1,...,7}\\ =-\frac{2c}{h}{\{\underset{s_p\∈s}{\mathrm{\Σ}}{k}^{\′}\left({\left|\right|\frac{r(s_p,p)}{h}\left|\right|}^2\right)\·r(s_p,p)\·\frac{\∂g(s_p,p)}{{\∂p}_i}\}}_{i=1,...,6}\end{array}---\left(8\right)$

$\begin{array}{c}\▿\▿\mathrm{\ρ}\left(p\right)\≈{\left(\frac{{\∂}^2\mathrm{\ρ}}{{\∂p}_i{\∂p}_j}\right)}_{\begin{array}{c}i=1,...,6\\ j=1,...,6\end{array}}\\ =-\frac{2c}{h}{\{\underset{s_p\∈s}{\mathrm{\Σ}}{k}^{\′}\left({\left|\right|\frac{r(s_p,p)}{h}\left|\right|}^2\right)\·\frac{\∂g(s_p,p)}{{\∂p}_i}\frac{\∂g(s_p,p)}{{\∂p}_j}\}}_{\begin{array}{c}i=1,...,6\\ j=1,...,6\end{array}}\end{array}---\left(9\right)$

In like manner, can formula (8) and formula (9) be brought into formula (3), and obtain the numerical solution minimizing formula (7) by iterative process.

The contrast of two kinds of methods

The First-order Gradient and the Hessian matrix that compare two kinds of method similarity functions can be found out, formal KDIC is the weighted version of TDIC, and weight function is equivalent to the first order derivative of kernel function.Therefore, the weighted version that different IPs CWinInetConnection type will be corresponding different.For intuitive analysis, discuss for the most frequently used two kinds of kernel functions.One is Epanechnikov core, and another kind is gaussian kernel.

Case1:Epanechnikov core

$k_E\left(x\right)=\left\{\begin{array}{c}1-x,\mathrm{if}0<x<1\\ 0,\mathrm{otherwise}\end{array}\right.$

That is,

$k_E^{\&prime;}\left(x\right)=\left\{\begin{array}{c}1,\mathrm{if}0<x<1\\ 0,\mathrm{otherwise}\end{array}\right.$

If Epanechnikov nucleus band to be entered formula (8) and (9), from expression-form, KDIC and TDIC is equivalent.That is, when adopting Epanechnikov core, TDIC is a special case of KDIC.

Case2: gaussian kernel

$k_G\left(x\right)=\exp (-\frac{1}{2}x),x>0$

And

$k_G^{\&prime;}\left(x\right)=\frac{1}{2}\exp (-\frac{1}{2}x),x>0$

Because the derivative of Gaussian function has identical profile with original function, we can utilize the mathematical characteristic of Gaussian function and and the physical meaning of reality between relation to realize weighting, thus obtain immunity to picture noise.Suppose that the noise in image is the white noise of stochastic distribution, and variance is σ.When accordingly, the bandwidth parameter (h in formula (6)) in kernel function can be set as the multiple of σ by us.Showing when deviate exceedes this multiple the to change the time pixel at place is very serious by noise pollution, distributes almost nil weights and directly can eliminate larger noise spot and affect.On the contrary, when deviate drops within the scope of this multiple, be weighted according to the deviate of noise pollution degree to difference by the mode of Gauss's weighting, distribute contribution intensity as each pixel in zoning, thus obtain measurement result more accurately.

Embodiment 1:

Based on the displacement field iteration smoothing method of the Digital Image Correlation Method of kernel function, comprise the following steps:

1), gather image before deformation of body and after distortion, the image before distortion is reference picture, and the image after distortion is target image;

2), region-of-interest S is chosen on a reference ₁, region-of-interest S ₁for reference picture subarea; Obtained by following steps: divide the first grid respectively on a reference, the number of the first grid be M capable × N row, spacing in the length and Width of reference picture between the first grid is respectively l and w, sets up the region-of-interest S centered by i-th the first grid on a reference ₁, region-of-interest S ₁for reference picture subarea, length and the width in reference picture subarea are respectively L and W, wherein, and a=1,2,3 ..., M × N.Utilize and choose region-of-interest by the mode of grid division on a reference, simple and convenient, easily realize.

3), on target image, region of search S is set up ₂, region of search S ₂for target image subarea, target image subarea comprises the reference picture subarea after distortion; Obtained by following steps: on target image, divide the second grid, the number of the second grid be M capable × N row, the spacing in the length and Width of target image between the second grid is respectively l and w; Target image is set up the region of search S centered by b the second grid ₂, region of search S ₂for target image subarea, length and the width in target image subarea are respectively k*L and k*W, wherein, and k>1, b=a.The mode by grid division on target image is utilized to choose region of search, simple and convenient, easily realize.Wherein, l and w is 3-5 pixel.

4) gray-scale value of the point in target image subarea and reference picture subarea, is obtained.Obtaining the gray-scale value of the point in target image subarea by carrying out bicubic spline interpolation to target image subarea, obtaining the gray-scale value of sub-pixel location in target image subarea.Obtaining the gray-scale value of the point in reference picture subarea also by carrying out bicubic spline interpolation to reference picture subarea, obtaining the gray-scale value of sub-pixel location in reference picture subarea.

5), build shape function, determine the position relationship of corresponding point in reference picture subarea and target image subarea, wherein, (x ₁, y ₁) be the coordinate of pixel any in reference picture subarea, (x ₂, y ₂) in target image subarea with pixel (x ₁, y ₁) coordinate of corresponding pixel, there is one group of mapping relations χ and following formula set up:

χ(x ₁,y ₁)→(x ₂,y ₂)

f(x ₁,y ₁)＝g(x ₂,y ₂)

Wherein, f (x ₁, y ₁) represent pixel (x ₁, y ₁) brightness of image at place, g (x ₂, y ₂) represent pixel (x ₂, y ₂) brightness of image at place, mapping relations χ is shape function; Wherein, shape function χ is expressed from the next:

$\begin{array}{c}{x}_2=x_1+u^{\&prime;}(x^{\&prime;},y^{\&prime;})+\frac{{\&PartialD;u}^{\&prime;}}{{\&PartialD;x}_1}{|}_{x^{\&prime;},y^{\&prime;}}(x_1-x^{\&prime;})+\frac{{\&PartialD;u}^{\&prime;}}{{\&PartialD;y}_1}{|}_{x^{\&prime;},y^{\&prime;}}(y_1-y^{\&prime;})\\ y_2=y_1+v^{\&prime;}(x^{\&prime;},y^{\&prime;})+\frac{{\&PartialD;v}^{\&prime;}}{{\&PartialD;x}_1}{|}_{x^{\&prime;},y^{\&prime;}}(x_1-x^{\&prime;})+\frac{{\&PartialD;v}^{\&prime;}}{{\&PartialD;y}_1}{|}_{x^{\&prime;},y^{\&prime;}}(y_1-y^{\&prime;})\end{array}$

Wherein, u and v is respectively the displacement of being out of shape x and the y direction caused, and (x', y') is the center position coordinates in reference picture subarea, for distortion is at x ₁and y ₁the single order displacement gradient that direction produces;

6), by step 5) shape function parametrization, and to represent with vectorial P',

$P^{\&prime;}=(u^{\&prime;},v^{\&prime;},\frac{\&PartialD;u^{\&prime;}}{\&PartialD;x_1},\frac{\&PartialD;v^{\&prime;}}{\&PartialD;x_1},\frac{\&PartialD;u^{\&prime;}}{\&PartialD;y_1},\frac{\&PartialD;v^{\&prime;}}{\&PartialD;y_1}),$

Definition related function ρ ^*,

${\mathrm{\&rho;}}^{*\&prime;}=c\underset{{s_p}^{\&prime;}\&Element;S_1}{\mathrm{\&Sigma;}}k\left({\left|\right|\frac{f\left({s_p}^{\&prime;}\right)-g({s_p}^{\&prime;},P^{\&prime;})}{h}\left|\right|}^2\right)$

Wherein, c is normaliztion constant, and k (.) is kernel function, and h is the bandwidth control parameter of kernel function, s _p' be the pixel of in reference picture subarea, f (s _p') be pixel s in image before distortion _p' place brightness of image, g (s _p', P') for after distortion in image with pixel s _p' the brightness of image at corresponding pixel place;

7), setting vector initial value P ₀', and set maximum iteration time n,

8), by P ₀' bring iterative equations into:

${P_1}^{\&prime;}={P_0}^{\&prime;}-\frac{{\&dtri;\mathrm{\&rho;}}^{*\&prime;}\left({P_0}^{\&prime;}\right)}{{\&dtri;\&dtri;\mathrm{\&rho;}}^{*\&prime;}\left({P_0}^{\&prime;}\right)}$

In formula, P ₀' be deformation parameter initial value, ▽ ρ ^*' and ▽ ▽ ρ ^*' be related function ρ ^*first-order Gradient and Hessian matrix, wherein ▽ ρ ^*' and ▽ ▽ ρ ^*' formula as follows:

${\&dtri;\mathrm{\&rho;}}^{*\&prime;}=\frac{-2}{\underset{s_p\&Element;S_1}{\mathrm{\&Sigma;}}{f}^2\left({s_p}^{\&prime;}\right)}{\left\{\underset{{s_p}^{\&prime;}\&Element;S_1}{\mathrm{\&Sigma;}}\right[f\left({s_p}^{\&prime;}\right)-g({s_p}^{\&prime;},P^{\&prime;})\left]\frac{\&PartialD;g({s_p}^{\&prime;},P^{\&prime;})}{{{\&PartialD;P}_i}^{\&prime;}}\right\}}_{i=1,...,6}$

${\&dtri;\&dtri;\mathrm{\&rho;}}^{*\&prime;}=\frac{-2}{\underset{s_p\&Element;S_1}{\mathrm{\&Sigma;}}{f}^2\left({s_p}^{\&prime;}\right)}{\left\{\underset{{s_p}^{\&prime;}\&Element;S_1}{\mathrm{\&Sigma;}}\frac{\&PartialD;g({s_p}^{\&prime;},P^{\&prime;})}{{{\&PartialD;P}_i}^{\&prime;}}\frac{\&PartialD;g({s_p}^{\&prime;},P^{\&prime;})}{{{\&PartialD;P}_j}^{\&prime;}}\right\}}_{\begin{array}{c}i=1,...,6\\ j=1,...,6\end{array}};$

In formula, P _i' and P _j' be respectively i-th element and a jth element of vectorial P';

P is calculated according to iterative equations ₁';

9), according to following formula judge whether iterative equations restrains, if convergence, then preserve region-of-interest S ₁shift value u' and v'; If convergence, does not make P ₀'=P ₁', and repeat step 8) and 9); If iterations reaches maximum iteration time n, then finishing iteration;

10) shift value of region-of-interest, is exported.

The input of this algorithm be speckle pattern before and after the test piece deformation that gathers, utilize the gray feature of image, final what export is the displacement field of test piece deformation; During calculating, the area-of-interest grid division point of test specimen before being deformed, object is exactly the position of asking for these net points image after deformation, thus can do the displacement that difference obtains these net point places.

Centered by the net point divided, divide reference picture subarea, to divide target image subarea in same central point after deformation image, size is k times of reference picture subarea, k>1.In order to the position at net point place after finding distortion, carry out point by point search, find the most similar subarea with reference to image subsection in target image subarea, so the center position in this subarea is the position of net point after deformation in image.

Utilize gradation of image feature construction similarity function, as minimum squared distance function etc.For improving the noise immunity of this similarity function, for this function adds weighted value items, in formula, k can select Gaussian function.

In formula, f is the gray-scale value of each point in reference picture subarea, and g is the gray-scale value of each point after distortion in image subsection.The position (x', y') of the point after distortion in image subsection is described by shape function, and (x', y') is sub-pixel location value, therefore needs to carry out interpolation to target subarea, as bicubic spline interpolation, obtains the gray-scale value of sub-pixel location.P=[U, V, Ux, Vx, Uy, Vy] ^tfor the unknown quantity in shape function, in fact, this similar function is the function of P, minimizes this function, can solve the value of P.

Similar function be a multivariable function, by nonlinear optimization method, as Newton-Raphson method, iterative can be carried out.

The iteration convergence condition of iterative equations: desirable ε < 0.1, under the prerequisite ensureing convergence, ε is less, and output displacement field is more accurate.Or convergence number of times is set, as 50 times, under the prerequisite ensureing convergence, this value is larger, and output displacement field is more accurate.

Sunykatuib analysis: for assessing the validity of said method, the method that digital speckle image simulation measures test piece deformation is used.Because deformation parameter is known, the quality for the treatment of effect can be evaluated objectively.Reference picture is the digital speckle image of stochastic generation, and resolution is 256 × 256pixels, and speckle particle number is 2000, and speckle particle radius is 4, as shown in Figure 2.By the deformation parameter preset, in the y-direction, generate an amplitude variation shape speckle pattern every 0.5pixel, symbiosis becomes 20.Be the noise that the image of 0-1pixel adds 1%, 2% and 4% respectively to this group displacement, add muting one group of speckle pattern, obtain 4 groups of speckle patterns altogether.Wherein, the speckle pattern of 4% noise level is added as shown in Figure 3.

For investigating the impact of different stage noise, first utilize traditional DIC method process four groups of speckle patterns, calculate shift value and employing two kinds of methods evaluated, be i.e. displacement average error and standard deviation, expression formula is as follows:

$e_v=\stackrel{\&OverBar;}{v}-v_{\mathrm{true}}$

${\mathrm{\&sigma;}}_v=\sqrt{\frac{1}{N-1}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(v_i-\stackrel{\&OverBar;}{v})}^2}$

Wherein, represent the mean value of all calculation level displacements in same width figure, namely v _truerepresentation theory deformation parameter; N represents that all calculating is counted.

Figure 4 and 5 show respectively the error and standard deviation that are obtained by NR method displacement calculating field post analysis.As seen from Figure 4, random noise calculates impact significantly to displacement field.Along with the increase of noise grade, displacement error entirety is in rising trend, from maximum error value, is increased to 0.0994pixel from 0.0077pixel.

Equally, as seen from Figure 5, the standard deviation of gained displacement field is also increase along with the increase of noise grade, and maximum standard deviation is increased to 0.0133pixel from 0.0029pixel.

For checking improves the validity of DIC method herein, KDIC method is used to calculating four groups emulation speckle pattern.Equally, displacement average error and standard deviation two kinds of evaluation indexes are used to the validity of appraisal procedure.Result of calculation as shown in Figure 6 and Figure 7.

As can see from Figure 7, contrast the displacement field error containing 4% noise level that traditional DIC calculates, the error of all noise levels calculated by KDIC method is all less, wherein for 4% noise level error of calculation analysis, maximum error drops to 0.0271 from 0.0094pixel, new method Be very effective.

For ease of comparing, maximum error and the average standard deviation of 4 groups of displacement result of calculations are listed in table 1.Result proves, the DIC method based on kernel function can effectively reduce the error of calculation of displacement field, and the larger effect of noise level is more remarkable.But from the result of calculation of standard deviation, treatment effect is not clearly.For 1% noise level, namely signal noise ratio (snr) of image is 40dB, the error of calculation of displacement field is 6 ‰, and the signal to noise ratio (S/N ratio) that the industrial camera being actually used in collection image can reach 40dB is substantially even higher, therefore, during the image of the method for the treatment of actual acquisition, the precision of thousand points of positions can meet the demands completely.

Table 1 TDIC and KDIC method result of calculation maximum error and the contrast of average standard deviation

Above simulation analysis demonstrates and improves the validity that DIC method reduces the displacement field error of calculation herein.For proving that the method effectively can process the image of actual acquisition, computational analysis will be carried out with KDIC method below.As shown in Figure 8, test specimen is aluminum material to actual picture, carries out the stretching of y direction.Treatment effect as shown in Figure 9.

Can see from result, the displacement field calculated through KDIC method is comparatively level and smooth, effectively can reduce the error that strain field that the later stage obtains by displacement field difference calculates, this result demonstrates improves the feasibility and validity that DIC method analyzes for actual experiment herein.

Refer to shown in Figure 11, the displacement field adaptive smooth method being applicable to DIC of the present invention.True according to the slickness of body surface, set up penalized least-squares regressive object function, and from noise data, estimate that penalty factor is to punish the roughening of solution by means of GCV method, thus realize effective smoothing processing of displacement field.The method have realize simple, calculated amount is little and full automatic advantage.

The deformation field utilizing DIC to carry out body structure surface is measured:

Deformation field DIC being used for body structure surface measures problem, mainly comprises three important steps.First, the suitable distortion of shape function description scheme under external load function is constructed; Then, set up certain correlation metric, be used for the degree of similarity of the image brightness distribution before and after the tested malformation of quantitative evaluation; Finally, solved by Multi-variables optimum design method and make the maximized shape function parameter of similarity criterion, thus indirect inspection go out distortion after displacement field data.

, there is one group of mapping relations χ and meet in the arbitrfary point (x, y) in given non-deformation pattern and the little neighborhood S of around one thereof and f (x, y) represents the brightness of image at point (x, y) place, represent the brightness of image at the rear corresponding coordinate place of distortion.If neighborhood S is enough little, mapping relations χ can be described by formula (1),

$\begin{array}{c}\tilde{x}=x+u(x_0,y_0)+\frac{\&PartialD;u}{\&PartialD;x}{|}_{x_0,y_0}\left({x-x}_0\right)+\frac{\&PartialD;u}{\&PartialD;y}{|}_{x_0,y_0}(y-y_0)\\ \tilde{y}=y+v(x_0,y_0)+\frac{\&PartialD;v}{\&PartialD;x}{|}_{x_0,y_0}(x-x_0)+\frac{\&PartialD;v}{\&PartialD;y}{|}_{x_0,y_0}(y-y_0)\end{array}---\left(10\right)$

Wherein, u and v is respectively the in-plane displacement in x and y direction, (x ₀, y ₀) be the center of region S.

Make parameter vector $P=(u,v,\frac{\&PartialD;u}{\&PartialD;x},\frac{\&PartialD;v}{\&PartialD;x},\frac{\&PartialD;u}{\&PartialD;y},\frac{\&PartialD;v}{\&PartialD;y}),$ And define related coefficient

$\mathrm{\&rho;}=\frac{\underset{(x,y)\&Element;S}{\mathrm{\&Sigma;}}{[f(x,y)-g(x,y,P)]}^2}{\underset{(x,y)\&Element;S}{\mathrm{\&Sigma;}}{f}^2(x,y)}---\left(11\right)$

As can be seen from the above equation, when related function minimalization, before and after distortion, the similarity of image subsection reaches maximal value, now, displacement parameter u and v that parameter vector P comprises represents the optimum estimate to displacement after distortion, in the same way all measurement points are calculated, whole audience displacement can be obtained.

For minimizing ρ, be the solution of zero by solving formula (10) First-order Gradient,

$\begin{array}{c}\&dtri;\mathrm{\&rho;}={\left(\frac{\&PartialD;\mathrm{\&rho;}}{{\&PartialD;P}_i}\right)}_{i=1,...,6}\\ =\frac{2}{\underset{(x,y)\&Element;S}{\mathrm{\&Sigma;}}{f}^2(x,y)}{\left\{\underset{(x,y)\&Element;S}{\mathrm{\&Sigma;}}\right[f(x,y)-g(x,y,P)\left]\frac{\&PartialD;g(x,y,P)}{\&PartialD;P}\right\}}_{i=1,...,6}=0\end{array}---\left(12\right)$

There is a lot of method to can be used to the formula that solves (12), adopt Newton-Raphson method iterative herein, have

$P=P_0-\frac{\&dtri;\mathrm{\&rho;}\left(P_0\right)}{\&dtri;\&dtri;\mathrm{\&rho;}\left(P_0\right)}---\left(13\right)$

In formula, P ₀for deformation parameter initial value, ▽ ▽ C (P ₀) be the Hessian matrix of related function ρ, meet

$\begin{array}{c}\&dtri;\&dtri;\mathrm{\&rho;}={\left(\frac{{\&PartialD;}^2\mathrm{\&rho;}}{{\&PartialD;P}_i{\&PartialD;P}_j}\right)}_{\begin{array}{c}i=1,...,6\\ j=1,...,6\end{array}}\\ =\frac{2}{\underset{(x,y)\&Element;S}{\mathrm{\&Sigma;}}{f}^2(x,y)}{\left\{\underset{(x,y)\&Element;S}{\mathrm{\&Sigma;}}\frac{\&PartialD;g(x,y,P)}{{\&PartialD;P}_i}\frac{\&PartialD;g(x,y,P)}{{\&PartialD;P}_j}\right\}}_{\begin{array}{c}i=1,...,6\\ j=1,...,6\end{array}}\end{array}---\left(14\right)$

As can be seen from above process, many factors can impact the accuracy of displacement field measurement result, such as subarea size, interpretational criteria, iteration convergence condition etc.Although DIC method obtains large quantifier elimination, the effective ways eliminating various systematic survey deviation are instructed also to lack very much.As hope obtains valuable strain field distribution data, need meticulously to process the displacement field data containing stochastic error.

Utilize N-R alternative manner to the smoothing process of displacement field:

Solve the shape function parameter making correlation coefficient ρ minimum, go out the displacement field data d after distortion according to the shape function parameter measurement obtained ⁿ, wherein displacement field data d ⁿbe expressed from the next: d ⁿ=d+n;

Wherein, d is desirable real displacement value, and n is displacement field noise;

Eliminating the stochastic error of the displacement field data U after distortion, eliminating stochastic error by minimizing following formula

$D(d,\mathrm{\&alpha;})=\mathrm{\&alpha;}{\left|\right|d-d^n\left|\right|}_2^2+(1-\mathrm{\&alpha;}){\left|\right|\mathrm{Cd}\left|\right|}_2^2$

Wherein, d ⁿfor the displacement field data after distortion, d is desirable real displacement value, and n is displacement field noise, and α is smoothing parameter; To above formula differentiate, obtain

$\begin{array}{c}\frac{\&PartialD;D(d,\mathrm{\&alpha;})}{\&PartialD;d}=2\mathrm{\&alpha;}(d-d^n)+2(1-\mathrm{\&alpha;})*4\left(\mathrm{Cd}\right),\\ \frac{\&PartialD;D(d,\mathrm{\&alpha;})}{\&PartialD;\mathrm{\&alpha;}}={\left|\right|d-d^n\left|\right|}_2^2-{\left|\right|\mathrm{Cd}\left|\right|}_2^2\end{array}$

Wherein, C is Laplace operator,

D and α is become column vector and makes Q=[d1, d2, d3...dm, α], wherein, m is the number of element in d; Then can obtain

$\&dtri;D\left(Q\right)=\left[\begin{array}{c}2\mathrm{\&alpha;}(d_i-d^n)+2(1-\mathrm{\&alpha;})*4\left({\mathrm{Cd}}_i\right)\\ ......\\ 2\mathrm{\&alpha;}(d_m-d^n)+2(1-\mathrm{\&alpha;})*4\left({\mathrm{Cd}}_m\right)\\ \mathrm{\&alpha;}\end{array}\right];$

Second-order partial differential coefficient ▽ ▽ D (Q) is tried to achieve according to ▽ D (Q),

$\frac{{\&PartialD;D}^2\left(Q\right)}{{\&PartialD;d}_i{\&PartialD;d}_j}=\left[\begin{array}{cccc}32-30\mathrm{\&alpha;}& & & \\ & 32-30\mathrm{\&alpha;}& & \\ & & ......& \\ & & & 32-30\mathrm{\&alpha;}\end{array}\right]$

$\frac{{\&PartialD;D}^2\left(Q\right)}{{\&PartialD;d}_i\&PartialD;\mathrm{\&alpha;}}=2(d-d^n)-8\left(\mathrm{Cd}\right)$

$\frac{{\&PartialD;D}^2\left(Q\right)}{{\&PartialD;\mathrm{\&alpha;}\&PartialD;d}_i}=2(d-d^n)-8\left(\mathrm{Cd}\right)$

$\frac{{\&PartialD;D}^2\left(Q\right)}{\&PartialD;\mathrm{\&alpha;}\&PartialD;\mathrm{\&alpha;}}=0$

Above-mentioned matrix is combined, obtains ▽ ▽ D (Q);

Utilize the ▽ D (Q) and ▽ ▽ D (Q) that obtain, iterative computation carried out to following formula:

$P_{k+1}=P_k-\frac{\&dtri;D\left(Q\right)}{\&dtri;\&dtri;D\left(Q\right)},$

Wherein, the initial value of Q is [d ⁿ, α] ^t, the condition of convergence is when the d obtained after above formula convergence is desirable real displacement value.

Simulation analysis: emulation speckle pattern parameter, 500pixel × 500pixel, speckle particle radius 4, speckle particle number 4000

(1) homogeneous deformation, x direction 1000 microstrain, adds noise in y direction displacement 0.3 image.In two kinds of situation:

The first situation: signal to noise ratio (S/N ratio) 40db

Table 1 displacement, strain by N-R compare with this paper error calculated

The first situation: signal to noise ratio (S/N ratio) 35db

Table 2 displacement, strain by N-R compare with this paper error calculated

Change the initial value of α, carry out iterative computation, α changes shown in Figure 12, and Figure 12 is under different α initial condition, α iteration situation of change

(2) inhomogeneous deformation: emulation speckle pattern parameter, 500pixel × 500pixel, speckle particle radius 4, speckle particle number 4000

Deformation: u (x, y)=0.1sin (2*pi*x/200) is in two kinds of situation: the first, signal to noise ratio (S/N ratio) 30; The second, signal to noise ratio (S/N ratio) 35

The first situation: signal to noise ratio (S/N ratio) 30db

Refer to shown in Figure 13, Figure 14 and Figure 15,

Under table 1, signal to noise ratio (S/N ratio) 30db condition, two kinds of method Output rusults application condition

Method	Std_v/pixel	Std_v y /με
			N-R alternative manner	0.02488	1227.8411
Distinguish level and smooth displacement t	0.02299	770.0738

The second situation, signal to noise ratio (S/N ratio) 35db

Refer to Figure 16, Figure 17 and Figure 18,

Under table 2, signal to noise ratio (S/N ratio) 35db condition, two kinds of method Output rusults application condition

Experimental analysis: rivet plate test specimen strain calculation

During experiment, rivet plate lower end is fixed, and upper end stretches, and theoretical analysis should be that rivet lower end is stressed maximum.

What Figure 19 and Figure 20 provided is rivet peripheral region linear deformation situation (can supplement other Direction distortion situations again)

(2) Multiple-Hole Specimen experiment

Refer to shown in Figure 21 and Figure 22,

Conclusion

Known by analysis, strain field accurately be obtained, two conditions need be met:

(1) displacement field of N-R calculating output is more accurate.If the displacement field deviation that N-R exports is too large, that strain result calculated must be inaccurate.Displacement field deviation is very large, and the data peaks after level and smooth can not be fallen a lot, because there is a balance between data accuracy and flatness; If displacement field data are less than normal, the value after that is level and smooth can be less, and error is also large.Therefore, N-R is necessary compared with the output of exact shift field.

(2) suitable displacement noise-reduction method.Obtain because strain is divided by displacement difference, the fluctuation of therefore displacement can bring comparatively big error to strain calculation.

Claims (10)

1. based on a displacement field iteration smoothing method for the Digital Image Correlation Method of kernel function, it is characterized in that: the Digital Image Correlation Method based on kernel function adopts the related function ρ based on kernel function ^*,

${\mathrm{\&rho;}}^{*}=c\underset{s_p\&Element;S_0}{\mathrm{\&Sigma;}}k\left({\left|\right|\frac{f\left(s_p\right)-g(s_p,P)}{h}\left|\right|}^2\right)$

Wherein, c is normaliztion constant, and k (.) is kernel function, and h is the bandwidth control parameter of kernel function, s _pfor being out of shape front image-region S ₀in pixel, f (s _p) be pixel s in image before distortion _pthe brightness of image at place, g (s _p, P) for after distortion in image with pixel s _pthe brightness of image at corresponding pixel place;

Wherein, iteration smoothing method comprises the following steps:

(1), the displacement field data d after distortion is measured based on the Digital Image Correlation Method of kernel function described in utilization ⁿ, wherein displacement field data d ⁿbe expressed from the next: d ⁿ=d+n;

Wherein, d be level and smooth after displacement field, n is displacement field noise;

(2), quadratic function is built $D(d,\mathrm{\&alpha;})=\mathrm{\&alpha;}{\left|\right|d-d^n\left|\right|}_2^2+(1-\mathrm{\&alpha;}){\left|\right|\mathrm{Cd}\left|\right|}_2^2$

Wherein, d ⁿfor the displacement field data after distortion, d be level and smooth after displacement field data, n is displacement field noise, and α is smoothing parameter, and C is high-order operator;

Solve the first-order partial derivative of quadratic function, obtain

$\frac{\&PartialD;D(d,\mathrm{\&alpha;})}{\&PartialD;d}=2\mathrm{\&alpha;}(d-d^n)+2(1-\mathrm{\&alpha;})*4\left(\mathrm{Cd}\right),$

$\frac{\&PartialD;D(d,\mathrm{\&alpha;})}{\&PartialD;\mathrm{\&alpha;}}={\left|\right|d-d^n\left|\right|}_2^2-{\left|\right|\mathrm{Cd}\left|\right|}_2^2$

Wherein, C is high-order operator,

D and α is become column vector and makes Q=[d1, d2, d3...dm, α], wherein, m is for counting; Then can obtain

$\&dtri;D\left(Q\right)=\left[\begin{array}{c}2\mathrm{\&alpha;}(d_1-{d_1}^n)+2(1-\mathrm{\&alpha;})*4\left(C{d}_1\right)\\ 2\mathrm{\&alpha;}(d_2-{d_2}^n)+2(1-\mathrm{\&alpha;})*4\left(C{d}_2\right)\\ ......\\ 2\mathrm{\&alpha;}(d_m-{d_m}^n)+2(1-\mathrm{\&alpha;})*4\left(C{d}_m\right)\\ {\left|\right|d-d^n\left|\right|}_2^2-{\left|\right|\mathrm{Cd}\left|\right|}_2^2\end{array}\right];$

(3), second-order partial differential coefficient ▽ ▽ D (Q) of quadratic function is solved according to the ▽ D (Q) of step (2),

$\frac{\&PartialD;D^2\left(Q\right)}{\&PartialD;d_i\&PartialD;d_j}=\left[\begin{array}{cccc}32-30\mathrm{\&alpha;}& & & \\ & 32-30\mathrm{\&alpha;}& & \\ & & ......& \\ & & & 32-30\mathrm{\&alpha;}\end{array}\right]$

$\frac{\&PartialD;D^2\left(Q\right)}{\&PartialD;d_i\&PartialD;\mathrm{\&alpha;}}=2(d-d^n)-8\left(\mathrm{Cd}\right)$

$\frac{\&PartialD;D^2\left(Q\right)}{\&PartialD;\mathrm{\&alpha;}\&PartialD;d_i}=2(d-d^n)-8\left(\mathrm{Cd}\right)$

$\frac{\&PartialD;D^2\left(Q\right)}{\&PartialD;\mathrm{\&alpha;}\&PartialD;\mathrm{\&alpha;}}=0$

Above-mentioned matrix is combined, obtains ▽ ▽ D (Q);

(4), the ▽ D (Q) that utilizes step (2) and (3) to obtain and ▽ ▽ D (Q), iterative computation is carried out to following formula:

$P_{k+1}=P_k-\frac{\&dtri;D\left(Q\right)}{\&dtri;\&dtri;D\left(Q\right)},$

Wherein, the initial value of Q is [d ⁿ, α] ^t, wherein, d ⁿfor the displacement field data after distortion, when the d obtained after above formula convergence is desirable real displacement value.

2. be applicable to the displacement field iteration smoothing method that digital picture is relevant as claimed in claim 1, it is characterized in that: the scope of described α is (0,1).

3. be applicable to the displacement field iteration smoothing method that digital picture is relevant as claimed in claim 1, it is characterized in that: the condition of convergence of step (4) is wherein, ε≤0.1.

4., as claimed in claim 1 based on the displacement field iteration smoothing method of the Digital Image Correlation Method of kernel function, it is characterized in that: the described Digital Image Correlation Method based on kernel function comprises the following steps:

One), shape function is defined: the arbitrfary point (x in setting reference picture ₀, y ₀) and neighborhood S around, (x, y) is the coordinate of any pixel in neighborhood S in reference picture, for the coordinate of pixel corresponding with pixel (x, y) in target image, there is one group of mapping relations χ and following formula set up:

$\mathrm{\&chi;}(x,y)\&RightArrow;(\tilde{x},\tilde{y})$

$f(x,y)=g(\tilde{x},\tilde{y})$

Wherein, f (x, y) represents the brightness of image at pixel (x, y) place, represent pixel the brightness of image at place, mapping relations χ is shape function; Wherein, shape function χ is expressed from the next:

$\tilde{x}=x+u(x_0,y_0)+\frac{\&PartialD;u}{\&PartialD;x}{|}_{x_0,y_0}(x-x_0)+\frac{\&PartialD;u}{\&PartialD;y}{|}_{x_0,y_0}(y-y_0)$

$\tilde{y}=y+v(x_0,y_0)+\frac{\&PartialD;v}{\&PartialD;x}{|}_{x_0,y_0}(x-x_0)+\frac{\&PartialD;v}{\&PartialD;y}{|}_{x_0,y_0}(y-y_0)$

Wherein, u and v is respectively the displacement of being out of shape x and the y direction caused, (x ₀, y ₀) be the center position coordinates of neighborhood S, for being out of shape the single order displacement gradient produced in the x and y direction;

Two), by step one) shape function parametrization, and to represent with vectorial P,

$P=(u,v,\frac{\&PartialD;u}{\&PartialD;x},\frac{\&PartialD;v}{\&PartialD;x},\frac{\&PartialD;u}{\&PartialD;y},\frac{\&PartialD;v}{\&PartialD;y}),$

Definition related function ρ ^*,

${\mathrm{\&rho;}}^{*}=c\underset{s_p\&Element;S}{\mathrm{\&Sigma;}}k\left({\left|\right|\frac{f\left(s_p\right)-g(s_p,P)}{h}\left|\right|}^2\right)$

Wherein, c is normaliztion constant, and k (.) is kernel function, and h is the bandwidth control parameter of kernel function, s _pfor the pixel of in neighborhood S, f (s _p) be pixel s in image before distortion _pthe brightness of image at place, g (s _p, P) for after distortion in image with pixel s _pthe brightness of image at corresponding pixel place;

Three), calculate and make related function ρ ^*shape function parameter P time minimum.

5., as claimed in claim 4 based on the displacement field iteration smoothing method of the Digital Image Correlation Method of kernel function, it is characterized in that: step 3) in calculate related function ρ by following formula ^*minimum value, structure iterative equations:

$P=P_0-\frac{\&dtri;{\mathrm{\&rho;}}^{*}\left(P_0\right)}{\&dtri;\&dtri;{\mathrm{\&rho;}}^{*}\left(P_0\right)}$

In formula, P ₀for deformation parameter initial value, ▽ ρ ^*with ▽ ▽ ρ ^*related function ρ ^*first-order Gradient and Hessian matrix, wherein ▽ ρ ^*with ▽ ▽ ρ ^*formula as follows:

$\begin{array}{c}\&dtri;{\mathrm{\&rho;}}^{*}={\left(\frac{\&PartialD;{\mathrm{\&rho;}}^{*}}{\&PartialD;P_i}\right)}_{i=1,...,6}\\ =\frac{-2}{\underset{s_p\&Element;S}{\mathrm{\&Sigma;}}{f}^2\left(s_p\right)}{\left\{\underset{s_p\&Element;S}{\mathrm{\&Sigma;}}\right[f\left(s_p\right)-g(s_p,P)\left]\frac{\&PartialD;g(s_p,P)}{\&PartialD;P_i}\right\}}_{i=1,...,6}\end{array}$

$\begin{array}{c}\&dtri;\&dtri;{\mathrm{\&rho;}}^{*}={\left(\frac{{\&PartialD;}^2{\mathrm{\&rho;}}^{*}}{\&PartialD;P_i\&PartialD;P_j}\right)}_{\begin{array}{c}i=1,...,6\\ j=1,...,6\end{array}}\\ =\frac{2}{\underset{s_p\&Element;S}{\mathrm{\&Sigma;}}{f}^2\left(s_p\right)}{\left\{\underset{s_p\&Element;S}{\mathrm{\&Sigma;}}\frac{\&PartialD;g(s_p,P)}{\&PartialD;P_i}\frac{\&PartialD;g(s_p,P)}{\&PartialD;P_j}\right\}}_{\begin{array}{c}i=1,...,6\\ j=1,...,6\end{array}}\end{array};$

Utilize iterative equations calculate and make related function ρ ^*minimum time solution.

6., as claimed in claim 1 based on the displacement field iteration smoothing method of the Digital Image Correlation Method of kernel function, it is characterized in that: the described Digital Image Correlation Method based on kernel function comprises the following steps:

1), gather image before deformation of body and after distortion, the image before distortion is reference picture, and the image after distortion is target image;

2), region-of-interest S is chosen on a reference ₁, region-of-interest S ₁for reference picture subarea;

3), on target image, region of search S is set up ₂, region of search S ₂for target image subarea, target image subarea comprises the reference picture subarea after distortion;

4) gray-scale value of the point in target image subarea and reference picture subarea, is obtained;

5), build shape function, determine the position relationship of corresponding point in reference picture subarea and target image subarea, wherein, (x ₁, y ₁) be the coordinate of pixel any in reference picture subarea, (x ₂, y ₂) in target image subarea with pixel (x ₁, y ₁) coordinate of corresponding pixel, there is one group of mapping relations χ and following formula set up:

χ(x ₁,y ₁)→(x ₂,y ₂)

f(x ₁,y ₁)＝g(x ₂,y ₂)

Wherein, f (x ₁, y ₁) represent pixel (x ₁, y ₁) brightness of image at place, g (x ₂, y ₂) represent pixel (x ₂, y ₂) brightness of image at place, mapping relations χ is shape function; Wherein, shape function χ is expressed from the next:

$x_2=x_1+u^{\&prime;}(x^{\&prime;},y^{\&prime;})+\frac{\&PartialD;u^{\&prime;}}{\&PartialD;x_1}{|}_{x^{\&prime;},y^{\&prime;}}(x_1-x^{\&prime;})+\frac{\&PartialD;u^{\&prime;}}{\&PartialD;y_1}{|}_{x^{\&prime;},y^{\&prime;}}(y_1-y^{\&prime;})$

$y_2=y_1+v^{\&prime;}(x^{\&prime;},y^{\&prime;})+\frac{\&PartialD;v^{\&prime;}}{\&PartialD;x_1}{|}_{x^{\&prime;},y^{\&prime;}}(x_1-x^{\&prime;})+\frac{\&PartialD;v^{\&prime;}}{\&PartialD;y_1}{|}_{x^{\&prime;},y^{\&prime;}}(y_1-y^{\&prime;})$

Wherein, u and v is respectively the displacement of being out of shape x and the y direction caused, and (x', y') is the center position coordinates in reference picture subarea, for distortion is at x ₁and y ₁the single order displacement gradient that direction produces;

6), by step 5) shape function parametrization, and to represent with vectorial P',

$P^{\&prime;}=(u^{\&prime;},v^{\&prime;},\frac{\&PartialD;u^{\&prime;}}{\&PartialD;x_1},\frac{\&PartialD;v^{\&prime;}}{\&PartialD;x_1},\frac{\&PartialD;u^{\&prime;}}{\&PartialD;y_1},\frac{\&PartialD;v^{\&prime;}}{\&PartialD;y_1}),$

Definition related function ρ ^*,

${\mathrm{\&rho;}}^{*\&prime;}=c\underset{{s_p}^{\&prime;}\&Element;S_1}{\mathrm{\&Sigma;}}k\left({\left|\right|\frac{f\left({s_p}^{\&prime;}\right)-g({s_p}^{\&prime;},P^{\&prime;})}{h}\left|\right|}^2\right)$

Wherein, c is normaliztion constant, and k (.) is kernel function, and h is the bandwidth control parameter of kernel function, s _p' be the pixel of in reference picture subarea, f (s _p') be pixel s in image before distortion _p' place brightness of image, g (s _p', P') for after distortion in image with pixel s _p' the brightness of image at corresponding pixel place;

7), setting vector initial value P ₀', and set maximum iteration time n,

8), by P ₀' bring iterative equations into:

${P_1}^{\&prime;}={P_0}^{\&prime;}-\frac{\&dtri;{\mathrm{\&rho;}}^{*\&prime;}\left({P_0}^{\&prime;}\right)}{\&dtri;\&dtri;{\mathrm{\&rho;}}^{*\&prime;}\left({P_0}^{\&prime;}\right)}$

In formula, P ₀' be deformation parameter initial value, ▽ ρ ^*' and ▽ ▽ ρ ^*' be related function ρ ^*first-order Gradient and Hessian matrix, wherein ▽ ρ ^*' and ▽ ▽ ρ ^*' formula as follows:

$\&dtri;{\mathrm{\&rho;}}^{*\&prime;}=\frac{-2}{\underset{s_p\&Element;S_1}{\mathrm{\&Sigma;}}{f}^2\left({s_p}^{\&prime;}\right)}{\left\{\underset{{s_p}^{\&prime;}\&Element;S_1}{\mathrm{\&Sigma;}}\right[f\left({s_p}^{\&prime;}\right)-g({s_p}^{\&prime;},P^{\&prime;})\left]\frac{\&PartialD;g({s_p}^{\&prime;},P^{\&prime;})}{\&PartialD;{P_i}^{\&prime;}}\right\}}_{i=1,...,6}$

$\&dtri;\&dtri;{\mathrm{\&rho;}}^{*\&prime;}=\frac{2}{\underset{s_p\&Element;S_1}{\mathrm{\&Sigma;}}{f}^2\left({s_p}^{\&prime;}\right)}{\left\{\underset{{s_p}^{\&prime;}\&Element;S_1}{\mathrm{\&Sigma;}}\frac{\&PartialD;g({s_p}^{\&prime;},P^{\&prime;})}{\&PartialD;{P_i}^{\&prime;}}\frac{\&PartialD;g({s_p}^{\&prime;},P^{\&prime;})}{\&PartialD;{P_j}^{\&prime;}}\right\}}_{\begin{array}{c}i=1,...,6\\ j=1,...,6\end{array}};$

In formula, P _i' and P _j' be respectively i-th element and a jth element of vectorial P';

P is calculated according to iterative equations ₁';

9), according to following formula judge whether iterative equations restrains, if convergence, then preserve region-of-interest S ₁shift value u' and v'; If convergence, does not make P ₀'=P ₁', and repeat step 8) and 9); If iterations reaches maximum iteration time n, then finishing iteration;

10) shift value of region-of-interest, is exported.

7. as claimed in claim 6 based on the displacement field iteration smoothing method of the Digital Image Correlation Method of kernel function, it is characterized in that: step 4) in obtain the gray-scale value of the point in target image subarea by carrying out bicubic spline interpolation to target image subarea, the gray-scale value of sub-pixel location in acquisition target image subarea.

8. as claimed in claim 6 based on the displacement field iteration smoothing method of the Digital Image Correlation Method of kernel function, it is characterized in that: step 4) in obtain the gray-scale value of the point in reference picture subarea by carrying out bicubic spline interpolation to reference picture subarea, the gray-scale value of sub-pixel location in acquisition reference picture subarea.

9., as claimed in claim 6 based on the displacement field iteration smoothing method of the Digital Image Correlation Method of kernel function, it is characterized in that: step 2) in choose region-of-interest S on a reference ₁, region-of-interest S ₁for reference picture subarea is obtained by following steps: divide the first grid respectively on a reference, the number of the first grid be M capable × N row, spacing in the length and Width of reference picture between the first grid is respectively l and w, sets up the region-of-interest S centered by a the first grid on a reference ₁, region-of-interest S ₁for reference picture subarea, length and the width in reference picture subarea are respectively L and W, wherein, and a=1,2,3 ..., M × N.

10., as claimed in claim 6 based on the displacement field iteration smoothing method of the Digital Image Correlation Method of kernel function, it is characterized in that: step 3) on target image, set up region of search S ₂, region of search S ₂for target image subarea, reference picture subarea after target image subarea comprises distortion is obtained by following steps: on target image, divide the second grid, the number of the second grid be M capable × N row, the spacing in the length and Width of target image between the second grid is respectively l and w; Target image is set up the region of search S centered by b the second grid ₂, region of search S ₂for target image subarea, length and the width in target image subarea are respectively k*L and k*W, wherein, and k>1, b=a.