CN100495022C - Concrete ultrasound tomography algorithm - Google Patents

Abstract

The disclosed concrete ultrasonic chromatography imaging algorithm comprises: providing a probabilistic ART algorithm for real-time 2D inversion imaging; wherein, preliminary determining the probability of ray to pass through the defect unit according to the projection value of every ray; then further determining the probability of every mesh as a defect unit to select the initial iterative value for ART algorithm; assigning the projection error during the iteration till stop. This invention can improve computation precision and image reconstruction quality to efficient take inversion the inner strength distribution and defect position and size for concrete.

Description

Concrete ultrasound tomography algorithm

Technical field

The invention belongs to concrete Ultrasonic NDT field, particularly concrete ultrasound tomography algorithm.

Technical background

Tomography (Computerized Tomography) is under the condition of not damaging research " object " inner structure, utilize certain radiographic source, according to outside with data for projection that checkout equipment obtained from " object ", according to certain physical and mathematical relation, utilize the distribution of certain physical quantity of the inner the unknown of computing machine inverting " object ", generate two dimension, 3-D view, reappear " object " internal feature.The main target that chromatography imaging technique is applied to the concrete Non-Destructive Testing is under the condition of not damaging concrete inner structure, determines the fine structure and the local unevenness of interior of building.

At present, ray path tomography inversion method relatively more commonly used mainly contains backprojection reconstruction algorithm, transform reconstruction class algorithm, discrete picture algebraic reconstruction class algorithm.Wherein, backprojection algorithm is fastest up to now a kind of algorithm, but computational accuracy is not high, is easy to generate " artefact "; With Fourier transform reconstruction algorithm and anti-(contrary) projection algorithm of filtering is the transform reconstruction method of representative, the noise resistance interference capability is poor, if and data for projection is not the simple integral along straight line, so have to be less than the closed form of resolving inversion formula, under these circumstances, it is invalid that converter technique just becomes, and therefore is not suitable for the imaging of concrete chromatography; Comparatively commonly used in the discrete picture algebraic reconstruction class algorithm have ART, SIRT, a constraint least square class algorithm (containing maximum entropy algorithm and optimization image reconstruction algorithm) etc., be applicable to that data for projection is incomplete, skewness, raypath is the occasion of curve, and be convenient to computer realization, therefore be used widely.Wherein, the SIRT algorithm only just demonstrates its superiority on reconstruction quality when measurement data is inaccurate especially, and other advantages and not obvious.And the ART algorithm is with respect to constraint least square class algorithm, and step is simple, is easy to programming and realizes, therefore, the most frequently used algorithm is the ART algorithm in concrete ultrasound tomography at present.The ART algorithm is through constantly improving, although simulation result and test findings are effectively, the precision of image reconstruction and rapidity still well do not solve.

Summary of the invention

The objective of the invention is to overcome above-mentioned prior art deficiency, a kind of concrete ultrasound tomography algorithm is provided, and this method can improve computational accuracy and computing velocity, by the result of two dimensional inversion imaging, effectively reflect concrete inner structure, thereby determine feature, size and the position of defective.

Technical scheme of the present invention is achieved in that concrete ultrasound tomography algorithm, adopts probability ART algorithm to carry out according to the following steps:

Step 1: by formula (1) obtains the velocity of wave of every ray

$v_i=\underset{j}{\mathrm{\Σ}}{a}_{\mathrm{ij}}/{\mathrm{\τ}}_i---\left(1\right)$

Wherein, v _iBe the velocity of wave of i bar ray, a _IjBe the ray length that i bar ray passes j grid, τ _iDuring for the i bar ray of measuring gained walking from the shot point to the acceptance point (during sound);

Step 2: according to ray velocity of wave Normal Distribution, then $\frac{\stackrel{\‾}{v}-V_t}{s/\sqrt{n}}~N\left(\mathrm{0,1}\right),$ Utilize the upside fractile table of standardized normal distribution, find fiducial probability P and be 100%, 90% ..., 10% o'clock pairing upside fractile λ _t(t=1,2 ... 10), obtain the lower limit of normal region velocity of wave then according to formula (2)

$\begin{array}{cc}{V}_t=\stackrel{\‾}{v}-{\mathrm{\λ}}_t\·s_v/\sqrt{n}& (t=\mathrm{1,2},\·\·\·,10)\end{array}---\left(2\right)$

Wherein, n is the ray sum by the survey district, and v is a velocity of wave mean value, s _vBe the velocity of wave standard deviation, v and s _vBy formula try to achieve (3) (4) respectively;

$\stackrel{\‾}{v}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{v}_i/n---\left(3\right)$

$s_v=\sqrt{\mathrm{\Σ}{(v_i-\stackrel{\‾}{v})}^2/(n-1)}---\left(4\right)$

Step 3: determine that according to formula (5) every ray passes the probability size of defective unit

$\begin{array}{cc}{p}_i=P_t{|}_{v_i<V_t}& (i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n)\end{array}---\left(5\right)$

Promptly working as the ray velocity of wave is P less than fiducial probability _t(t=1,2 ... 10) the normal velocity of wave lower limit of correspondence the time, then to pass the probability of defective unit be P to this ray _t

Step 4: determine that according to formula (6) each grid is the probability size of defective

Wherein, α is parameter factors (0＜α≤1),

=∑

_(i)/ lj, lj are the number of lines of penetrating that passes j grid,

Be the line segment length a that i bar ray passes j grid _IjRatio with the total length of i bar ray;

Step 5: the initial velocity of wave of choosing iteration according to formula (7)

$\begin{array}{cc}{v}_j^{\left(0\right)}={\mathrm{<figure-callout\; id="1"\; label="V"\; filenames="C200610104463E00132.png"\; state="\{\{state\}\}">V</figure-callout>}}_1\&CenterDot;q_j+{\mathrm{<figure-callout\; id="10"\; label="V"\; filenames="C200610104463E00141.png"\; state="\{\{state\}\}">V</figure-callout>}}_{10}\&CenterDot;(1-q_j)& (j=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,m)\end{array}---\left(7\right)$

Promptly obtaining the slow initial value of ripple is

${{\hat{f}}_j}^{\left(0\right)}=1/{v_j}^{\left(0\right)}---\left(8\right)$

Step 6: i bar ray to j the slow estimated value of grid ripple is when remembering q wheel iteration

Use (9) formula, one by one ray i (i=1,2 ..., n) to the correction of ripple slow-motion row, wherein b _Ij=a _IjQ _j, distributing to walk time error as weights, μ is the Poisson factor, 0＜μ≤1, a _IjBe the line segment length of i bar ray in j image-generating unit, τ _iDuring for the i bar ray of measuring gained walking from the shot point to the acceptance point (during sound),

For utilizing

And a _IjDuring the i bar ray that calculates gained walking from the shot point to the acceptance point (during sound);

${\hat{f}}_j^{q,i+1}={\hat{f}}_j^{q,i}+\mathrm{\&Delta;}{f}_j={\hat{f}}_j^{q,i}+\mathrm{\&mu;}\frac{b_{\mathrm{ij}}\&CenterDot;q_j}{\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{b}_{\mathrm{ij}}^2}({\mathrm{\&tau;}}_i-{\widehat{\mathrm{\&tau;}}}_i^q)---\left(9\right)$

Step 7: to the slowness vector of being tried to achieve

Carrying out degree of convergence judges:

‖f ^q-f ^q-1‖ _∞<ε (10)

Wherein,

Be the slowness vector that q wheel iteration obtains, the error bound of ε for setting are positive numbers, then stop iteration if following formula is set up, and change step 6, carry out q+1 wheel iteration, till satisfying convergence criterion.

The present invention can directly apply to concrete Non-Destructive Testing scene, and object to be detected is carried out real-time two dimensional inversion imaging, thereby accurately reflects concrete inner structure.

It is probability ART algorithm that the present invention adopts inversion algorithm.Determine tentatively that according to the projection value (when walking) of every ray every ray passes the probability of defective unit, determine further that then each grid is the probability of defective unit, utilize this probable value to choose the iteration initial value of ART algorithm, and in iterative process, reasonably distribute projection (when walking) error according to the probability size, shortage probability is big more, distribute to walk time error many more, the less number of times of iteration just can obtain the imaging effect of better quality, judgement to defective locations, shape and size is more accurate, thereby has effectively improved the speed and the precision of image reconstruction.

Description of drawings

Fig. 1 is a probability ART algorithm flow chart of the present invention;

Fig. 2 is that the present invention singly surveys the detection mode synoptic diagram;

Fig. 3 is a computer simulation experiment model sectional view of the present invention;

Fig. 4 is the velocity of wave 3-D display figure of computer simulation experiment of the present invention, wherein, figure (a) is model velocity of wave 3-D display figure, and figure (b) is the velocity of wave 3-D display figure that obtains for 100 times with traditional ART algorithm iteration, and figure (c) is the velocity of wave 3-D display figure that obtains for 50 times with probability ART algorithm iteration;

Fig. 5 is the section velocity of wave distribution plan that concrete sample obtains for 100 times with traditional ART algorithm iteration;

Fig. 6 is the section velocity of wave distribution plan that concrete sample obtains for 50 times with probability ART algorithm iteration.

Below in conjunction with accompanying drawing content of the present invention is described in further detail.

Embodiment

The iteration initial value of each grid cell of the present invention no longer be rely on priori compose to identical value, but determine according to the shortage probability of each grid; In the process of iteration, each bar ray walk all unit that time error is distributed to the ray process no longer fifty-fifty, but come reasonable distribution projection (when walking) error according to the shortage probability size of grid, shortage probability is big more, distribute to walk time error many more, thereby improved the speed and the precision of image reconstruction.

With reference to shown in Figure 1, its concrete steps are as follows:

1) ray passes determining of defective unit probability

$v_i=\underset{j}{\mathrm{\&Sigma;}}{a}_{\mathrm{ij}}/{\mathrm{\&tau;}}_i;---\left(1\right)$

$\stackrel{\&OverBar;}{v}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{v}_i/n;---\left(2\right)$

$s_v=\sqrt{\mathrm{\&Sigma;}{(v_i-\stackrel{\&OverBar;}{v})}^2/(n-1)};---\left(3\right)$

v _iThe velocity of wave of-the i bar ray; V-velocity of wave mean value; s _v-velocity of wave standard deviation;

In same survey district, suppose each bar ray velocity of wave Normal Distribution, then $\frac{\stackrel{\&OverBar;}{v}-V_t}{s/\sqrt{n}}~N\left(\mathrm{0,1}\right),$ Utilize the upside fractile table of standardized normal distribution, find fiducial probability P and be 100%, 90% ..., 10% o'clock pairing upside fractile λ _t(t=1,2 ... 10), and then to obtain fiducial probability be P _t(t=1,2 ... 10) lower limit of normal region velocity of wave the time, $V_t=\stackrel{\&OverBar;}{v}+{\mathrm{\&lambda;}}_t\&CenterDot;s_v/\sqrt{n}$ (t＝1，2，…，10) (4)

If the ray velocity of wave is P less than fiducial probability _tThe time normal velocity of wave lower limit, then to pass the probability of defective unit be P to this ray _t, can obtain the probability that every ray passes defective unit, be calculated as follows

$\begin{array}{cc}{p}_i=P_t{|}_{v_i<V_t}& (i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n)\end{array}---\left(5\right)$

2) grid is determining of defective unit probability

The probability that passes defective unit as if a certain ray is p _i, think that then all grids that it passes are that the probability of defective unit is p _i, and with P _iAssignment is given its all grid that pass through, owing to have different rays to pass in each grid, so each grid will obtain one group of probability data, and problem is converted into how according to a different set of probable value (p ₍₁₎, p ₍₂₎..., p _(l)) determine that grid is the probability q of defective unit _jp _iHour, illustrate that all grids that i bar ray passes are the probability q of defective _j(j=1,2 ... m) all less, otherwise and, if q _iBigger, can only illustrate that then indivedual grids that it passes are the probability q of defective _jBigger, so grid is the probability q of defective _jShould choose (p ₍₁₎, p ₍₂₎..., p _(l)) in less relatively data, but since grid to the difference that influences of different rays, q _jChoose the line segment length a of every ray in j grid that should consider to pass j grid _IjRatio with the length overall of ray

Size.So with

Be weight coefficient, adopt the method for filtering to choose in the probability data than fractional value as q _j

To pass the p of its each bar ray for each grid _iAnd it is to each bar ray

Corresponding one by one assignment is in it, and presses p _iSize ordering, obtain one group of data { (p ₍₁₎,

), (p ₍₂₎,

) ..., (p _(lj),

), p wherein ₍₁₎≤ p ₍₂₎≤ ... ≤ p _(lj), l _jFor passing the number of lines of penetrating of j grid.Choose the q of each grid by formula (7) _j, wherein, α is a parameter factors.

(6)

3) the iteration initial value chooses

Choose the velocity of wave initial value of each grid according to the probability size

$\begin{array}{cc}{v}_j^{\left(0\right)}={\mathrm{<figure-callout\; id="1"\; label="V"\; filenames="C200610104463E00132.png"\; state="\{\{state\}\}">V</figure-callout>}}_1\&CenterDot;q_j+{\mathrm{<figure-callout\; id="10"\; label="V"\; filenames="C200610104463E00141.png"\; state="\{\{state\}\}">V</figure-callout>}}_{10}\&CenterDot;(1-q_j)& (j=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,m)\end{array}---\left(8\right)$

Then the slow initial value of ripple is

${{\hat{f}}_j}^{\left(0\right)}=1/{v_j}^{\left(0\right)}---\left(9\right)$

4) ripple is revised determining of increment slowly

With each grid is in the probability introducing error distribution of defective unit, with a _IjWith q _jProduct be that weights distribute to walk time error, obtain the slow modification increment of following formula ripple:

$\begin{array}{cc}\mathrm{\&Delta;}{f}_j=\mathrm{\&mu;}\frac{b_{\mathrm{ij}}\&CenterDot;q_j}{\underset{j=1}{\overset{\mathrm{<figure-callout\; id="2"\; label="m\; b\; ij"\; filenames="C200610104463E00132.png"\; state="\{\{state\}\}">m</figure-callout>}}{\mathrm{\&Sigma;}}}{b}_{\mathrm{ij}}^{}{}_2^{}}{\mathrm{\&Delta;\&tau;}}_i& (i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n;j=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,m)\end{array}---\left(10\right)$

Wherein, b _Ij=a _IjQ _j

I bar ray to j the slow estimated value of grid ripple is when remembering q wheel iteration

Use (11) formula, one by one ray i (i=1,2 ..., n) ripple is modified as follows slowly:

${\hat{f}}_j^{q,i+1}={\hat{f}}_j^{q,i}+\mathrm{\&Delta;}{f}_j={\hat{f}}_j^{q,i}+\mathrm{\&mu;}\frac{b_{\mathrm{ij}}\&CenterDot;q_j}{\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{b}_{\mathrm{ij}}^2}({\mathrm{\&tau;}}_i-{\widehat{\mathrm{\&tau;}}}_i^q)---\left(11\right)$

5) the slowness vector to being tried to achieve

Carrying out degree of convergence judges:

‖f ^q-f ^q-1‖ _∞<ε (12)

Wherein,

Be the slowness vector that q wheel iteration obtains, the error bound of ε for setting are positive numbers, then stop iteration if following formula is set up, otherwise, carry out q+1 wheel iteration, till satisfying convergence criterion.

The present invention has the following advantages: probability ART algorithm has iterative convergence speed and higher computational accuracy faster, has improved the image reconstruction quality, effectively is finally inversed by the size and the position of the intensity distributions and the defective of concrete inner structure.

With reference to shown in Figure 2, arrange transmitting transducer (R respectively in the both sides of timing length direction ₁, R ₂... R _m) and receiving transducer (T ₁, T ₂... T _m), each transducer is placed on the midpoint on grid cell border, wherein, and T _iRepresent i transmitting transducer, R _jRepresent j receiving transducer.

Embodiments of the invention one:

With reference to shown in Figure 3, survey the district and be the rectangular region of 100cm * 60cm, be divided into 10 * 6 grid (image-generating unit), its size is 10cm * 10cm, with grid according to from top to bottom, from left to right serial number, wherein 22,23,28,29,34, No. 35 is the defective unit lattice, the defective velocity of wave is 4050m/s, normal velocity of wave is 4500m/s, arranges emission and receiving transducer in the both sides of length direction.

Shown in Fig. 4 a-c, use ART algorithm iteration 100 times and probability ART algorithm iteration 50 times, the result of calculation contrast is as shown in table 1, and among the figure, the x axle is represented timing length, and unit is a rice, and sector width is surveyed in the representative of y axle, and unit is a rice, and the z axle is represented velocity of wave, and unit is a meter per second.As seen, the result of calculation of probability ART algorithm is more near the computer-experiment model, and normal cell velocity of wave difference is very little, and defective locations is more outstanding, makes computational accuracy improve greatly.

Table 1 probability ART and the contrast of traditional ART result of calculation

Embodiments of the invention two:

Do simulation test with concrete sample, the test specimen section is surveyed head of district 44cm, dark 50cm, be divided into 11 * 10 grids, wherein, core is defective (a concrete folder mud), and size is 15 * 15cm, arrange emission and receiving transducer in the both sides of length direction respectively, when obtaining 11 * 11 sound, utilize ART algorithm and the inverting of probability ART algorithm respectively, the section velocity of wave that obtains distributes shown in accompanying drawing 5 and accompanying drawing 6.

With reference to shown in Figure 5, the cross-sectional imaging of traditional ART algorithm is deviation to some extent, and two defect areas have appearred in the test specimen center, are that an even defective conforms to not to the utmost with the concrete sample center.

With reference to shown in Figure 6, probability ART algorithm image quality improves greatly, and the position of defective and shape with outstanding, are that an even defective is consistent with the concrete center more accurately.

Claims (1)

1, concrete ultrasound tomography algorithm is characterized in that, adopts probability ART algorithm to carry out according to the following steps:

Step 1: by formula (1) obtains the velocity of wave of every ray

$v_i=\underset{j}{\mathrm{\&Sigma;}}{a}_{\mathrm{ij}}/{\mathrm{\&tau;}}_i---\left(1\right)$

Wherein, v _iBe the velocity of wave of i bar ray, a _IjBe the ray length that i bar ray passes j grid, τ _iDuring for the i bar ray of measuring gained walking from the shot point to the acceptance point (during sound);

Step 2: according to ray velocity of wave Normal Distribution, then $\frac{\stackrel{\&OverBar;}{v}-V_t}{s/\sqrt{n}}~N\left(\mathrm{0,1}\right),$ Utilize the upside fractile table of standardized normal distribution, find fiducial probability P and be 100%, 90% ..., 10% o'clock pairing upside fractile λ _t(t=1,2 ... 10), obtain the lower limit of normal region velocity of wave then according to formula (2);

$\begin{array}{cc}{V}_t=\stackrel{\&OverBar;}{v}-{\mathrm{\&lambda;}}_t\&CenterDot;s_v/\sqrt{n}& (t=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,10)\end{array}---\left(2\right)$

Wherein, n is the ray sum by the survey district, and v is a velocity of wave mean value, s _vBe the velocity of wave standard deviation, v and s _vBy formula try to achieve (3) (4) respectively;

$\stackrel{\&OverBar;}{v}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{v}_i/n---\left(3\right)$

$s_v=\sqrt{\mathrm{\&Sigma;}{(v_i-\stackrel{\&OverBar;}{v})}^2/(n-1)}---\left(4\right)$

Step 3: determine that according to formula (5) every ray passes the probability size of defective unit

$\begin{array}{cc}{p}_i=P_i{|}_{v_i<V_t}& (i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n)\end{array}---\left(5\right)$

Promptly working as the ray velocity of wave is P less than fiducial probability _t(t=1,2 ... 10) the normal velocity of wave lower limit of correspondence the time, then to pass the probability of defective unit be P to this ray _t

Step 4: determine that according to formula (6) each grid is the probability size of defective

Wherein, α is parameter factors (0＜α≤1),

Lj is the number of lines of penetrating that passes j grid,

Be the line segment length a that i bar ray passes j grid _IjRatio with the total length of i bar ray;

Step 5: the initial velocity of wave of choosing iteration according to formula (7)

$\begin{array}{cc}{v}_j^{\left(0\right)}=V_1\&CenterDot;q_j+V_{10}\&CenterDot;(1-q_j)& (j=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,m)\end{array}---\left(7\right)$

Promptly obtaining the slow initial value of ripple is

${\hat{f}}_j^{\left(0\right)}=1/{v_j}^{\left(0\right)}---\left(8\right)$

Step 6: i bar ray to j the slow estimated value of grid ripple is when remembering q wheel iteration

Use (9) formula, one by one ray i (i=1,2 ..., n) to the correction of ripple slow-motion row, wherein b _Ij=a _IjQ _j, distributing to walk time error as weights, μ is the Poisson factor, 0＜μ≤1, a _IjBe the line segment length of i bar ray in j image-generating unit, τ _iDuring for the i bar ray of measuring gained walking from the shot point to the acceptance point or during sound, For utilizing

And a _IjDuring the i bar ray that calculates gained walking from the shot point to the acceptance point or during sound;

${\hat{f}}_j^{q,i+1}={\hat{f}}_j^{q,i}+\mathrm{\&Delta;}{f}_j={\hat{f}}_j^{q,i}+\mathrm{\&mu;}\frac{b_{\mathrm{ij}}\&CenterDot;q_j}{\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{b}_{\mathrm{ij}}^2}({\mathrm{\&tau;}}_i-{\widehat{\mathrm{\&tau;}}}_i^q)---\left(9\right)$

Step 7: to the slowness vector of being tried to achieve

Carrying out degree of convergence judges:

‖f ^q-f ^q-1‖ _∞<ε (10)

Wherein,

Be the slowness vector that q wheel iteration obtains, the error bound of ε for setting are positive numbers, then stop iteration if following formula is set up, and change step 6, carry out q+1 wheel iteration, till satisfying convergence criterion.

Images