1. An ultrasonic image restoration method based on point spread function parameter optimization is characterized by comprising the following steps:

s1: acquiring an ultrasonic C-scan image, and removing extreme value noise points on the ultrasonic C-scan image to obtain an ultrasonic image q to be restored;

s2: iterating the ultrasonic image f and the point spread function h to obtain the optimal point spread function h_s ；

S3: according to the point spread function h with the optimal parameters_sA corresponding restored ultrasound image is obtained.

The iterative solution formula of the ultrasonic image f and the point spread function h is expressed as follows:

k represents the number of alternating iterations, and k is 0,1,2,3 …; f. of_kAnd h_kRespectively representing the ultrasonic image and the point spread function, f, obtained in the k-th iterative solution_k+1And h_k+1Respectively representing an ultrasonic image and a point spread function obtained in the (k + 1) th iterative solution; h is_sThe point spread function with the optimal parameters is represented, s represents a parameter vector of the ultrasonic imaging system, and g represents an ultrasonic image to be restored; equation (1) adopts RL-TV algorithm of total variation regularization constraint to solve f_k+1(ii) a The formula (2) is based on a multivariate Gaussian sound beam model MGB to obtain the parameter optimization expression of the point spread function, and h is carried out by adopting a particle swarm optimization algorithm_k+1And (4) carrying out optimization solution on the parameters.

2. The method of claim 1, wherein the step S2 includes:

s21: initialization f₀Initializing h as g₀(x，y)：

Wherein, (x, y) is the pixel coordinate in the ultrasonic image g to be restored;

s22: solving the intermediate ultrasonic image f by adopting RL-TV algorithm_k+1The iterative format of the RL-TV algorithm is:

where n denotes the number of iterations, and n is 0,1,2,3 …, f_n(x, y) is the ultrasound image obtained in the nth iteration, f_n+1(x, y) is an ultrasonic image obtained by the (n + 1) th iteration, h (x, y) is a given point spread function, and h (-x, -y) represents a result obtained after the h (x, y) is turned along the x axis and the y axis; g (x, y) represents an ultrasound image to be restored; lambda [ alpha ]_TVFor regularization coefficients, div is the divergence calculation;

step S23: performing alternate solving on the formula (1) and the formula (2) through the formula (4), judging whether the maximum iteration number is reached according to the result of the alternate solving, if the maximum iteration number is reached, turning to the step S25, otherwise, turning to the step S24;

s24: solving intermediate point diffusion function h by adopting particle swarm optimization algorithm_k+1And go to step S22 to continue the iteration;

s25: and obtaining the point spread function with the optimal parameters and the corresponding restored ultrasonic image.

3. The method of claim 2, wherein the step S24 includes:

s241: for the sound field of water immersion ultrasonic detection, when the incident direction of ultrasonic waves is vertical to a water-solid interface, the sound field distribution in a solid workpieceComprises the following steps:

wherein the content of the first and second substances,representing particles (x) in a solid workpiece₂，y₂，z₂) ω ═ 2 pi f, which represents the ultrasonic probe circle frequency; j represents an imaginary number;

v₀(ω) -surface particle vibration velocity of the ultrasonic probe;

gamma 2-wave form in solid; p-longitudinal wave;

ρ₁water density, p₂-a solid workpiece density;

the speed of sound of longitudinal waves in water,-longitudinal acoustic velocity in solids;

a-probe radius;

A_r，B_r-multivariate gaussian superposition complex coefficients;

z₁-water layer depth between the ultrasonic probe and the workpiece surface;

s242: after obtaining the sound field distribution of the ultrasonic probe, calculating the inherent z of the sound field₀Point spread function at depth h (x, y):

whereinIs the sound field distribution obtained in step S341, v₀Representing the vibration velocity of particle on the surface of the ultrasonic probe, p representing the density of the solid workpiece,

s243: for the radius a of the ultrasonic probe and the center frequency f of the ultrasonic probe₀And depth z of investigation of sound field in workpiece₂Optimizing, comprising: let the characteristic parameter vector s ═ z₂，a，f₀]^TThe point spread function with the optimal parameters determined by the feature parameter vector S according to step S242 is recorded as h_s＝MGB(s)；

The optimization objective function when the particle swarm optimization algorithm is adopted for parameter optimization is as follows:

if the number of particle populations is I and the maximum number of iterations is T, the position of the I-th (I-1, 2,3, …, I) particle at the T (T-1, 2,3, …, T) iteration is:

the moving speed of the ith particle at the t iteration is:

the optimal position searched by the ith particle at present is as follows:

the optimal position currently searched by the whole particle population is as follows:

G^t＝(p_G1，p_G2，p_G3)； (11)；

the particles update speed and position according to equation (12):

where eta represents the inertia factor, lambda₁、λ₂Are self-learning factor and social learning factor, respectively, delta₁、δ₂Are all [0,1]A random number within a range;

s244: the parameter optimizing step of the point spread function based on the particle swarm optimization algorithm comprises the following steps:

s2441: setting the number I of particle populations, the maximum iteration number T and the initial positions of particlesInitial velocity of particlesParticle position boundary [ P ]_L，P_U]And particle velocity boundary [ V ]_L，V_U]；

S2442: calculating the fitness values of all the particles of the current iteration according to an objective function F(s)

S2443: comparing the current position fitness value of each particleA fitness value corresponding to the optimal position of the particleIf it is notUpdating the optimal position of the particle to be the current position, otherwise, keeping the optimal position unchanged;

s2444: comparing the current position fitness value of each particleFitness value F (G) of optimal position of particle swarm^t) If, ifUpdating the optimal position of the particle swarm to be the current position of the particle, otherwise, keeping the optimal position unchanged;

s2445: updating the velocity of individual particlesAnd position

S2446: judging whether the maximum iteration times is reached, if so, performing step S2447, otherwise, turning to step S2442;

s2447: outputting the current optimal result h_s。

Background

With the continuous development of ultrasonic imaging technology, ultrasonic C-scan imaging technology is increasingly used for measuring the size of defects inside workpieces. However, due to the resolution of the ultrasound imaging system, the ultrasound C-scan image is very blurred, from which it is difficult to determine the exact defect shape and size. The ultrasound image g (x, y) is generally regarded as the result of a two-dimensional convolution of the true defect distribution f (x, y) with the Point Spread Function (PSF) h (x, y) of the ultrasound imaging system: g (x, y) ═ f (x, y) × h (x, y) + n (x, y) #; where x represents a two-dimensional convolution operation and n (x, y) is the image noise incorporated during imaging. The essence of image restoration is the process of solving for the true defect distribution (i.e., sharp image) f (x, y) by deconvolution.

An image restoration method using a Point Spread Function (PSF) known in advance is called an image non-blind restoration method, but in such an image restoration method, it is also a difficult point that an accurate point spread function of an ultrasound imaging system is obtained in advance as the most important step. Therefore, when the point spread function of the ultrasound imaging system is not accurately obtained in advance, how to restore the ultrasound image is an urgent problem to be solved.

Disclosure of Invention

The utility model provides an ultrasonic image restoration method based on point spread function parameter optimization, which aims to directly obtain a clear image and a point spread function of an ultrasonic imaging system through a fuzzy ultrasonic image and restore the ultrasonic image under the condition of not presetting the point spread function of the ultrasonic imaging system.

The technical purpose of the present disclosure is achieved by the following technical solutions:

an ultrasonic image restoration method based on point spread function parameter optimization comprises the following steps:

s1: acquiring an ultrasonic C-scan image, and removing extreme value noise points on the ultrasonic C-scan image to obtain an ultrasonic image q to be restored;

s2: iterating the ultrasonic image f and the point spread function h to obtain the optimal point spread function h_s ；

S3: according to the point spread function h with the optimal parameters_sA corresponding restored ultrasound image is obtained.

The iterative solution formula of the ultrasonic image f and the point spread function h is expressed as follows:

k represents the number of alternating iterations, and k is 0,1,2,3 …; f. of_kAnd h_kRespectively representing the ultrasonic image and the point spread function, f, obtained in the k-th iterative solution_k+1And h_k+1Respectively representing an ultrasonic image and a point spread function obtained in the (k + 1) th iterative solution; h is_sThe point spread function with the optimal parameters is represented, s represents a parameter vector of the ultrasonic imaging system, and g represents an ultrasonic image to be restored; equation (1) adopts RL-TV algorithm of total variation regularization constraint to solve f_k+1(ii) a The formula (2) is based on a multivariate Gaussian sound beam model MGB to obtain the parameter optimization expression of the point spread function, and h is carried out by adopting a particle swarm optimization algorithm_k+1And (4) carrying out optimization solution on the parameters.

The beneficial effect of this disclosure lies in: the ultrasonic image and the point spread function are simultaneously solved through the alternate minimization framework, and the RL-TV algorithm is adopted to complete the solving part of the intermediate image in the alternate minimization; and then, a multi-parameter optimization solving part of the point spread function in the middle of the alternate minimization is completed by adopting a particle swarm optimization algorithm, and finally, a restored ultrasonic image and a corresponding point spread function are obtained. The method does not need to preset the point spread function of the ultrasonic imaging system, and can directly obtain the clear image and the point spread function of the ultrasonic imaging system only through the fuzzy ultrasonic image.

The method has the beneficial effects that:

(1) the method overcomes the defect that the image restoration effect of the traditional non-blind restoration method depends on the point spread function which is estimated accurately in advance, and has higher practicability.

(2) According to the method, the sound field distribution of the ultrasonic probe is obtained based on the multivariate Gaussian sound beam model, and then the parametric representation of the point spread function of the ultrasonic imaging system is obtained, so that the problem of insufficient theoretical basis of the point spread function obtained in the traditional image blind restoration and non-blind restoration methods is solved.

(3) The point spread function is optimized through the parameters of the point spread function instead of directly substituting the actual parameters into the model to obtain the point spread function, and the problem that the accuracy of the point spread function is not high due to the deviation of an MGB theoretical sound field model and the actual sound field is solved.

(4) The method and the device have the advantages that parameter optimization is carried out on the point spread function in the limited freedom degree space, the freedom degree of solving the point spread function is obviously reduced, the solving complexity of an optimization problem is reduced, meanwhile, overfitting can be avoided, and the overall stability of the solving process is improved.

(5) The method and the device have higher reference and reference meanings for restoring images in other related fields which can be represented by point spread function parameterization of the ultrasonic imaging system.

Drawings

FIG. 1 is a flow chart of a method described herein;

FIG. 2 is a flow chart of an embodiment of a method described herein;

FIG. 3 is a C-scan image of a stainless steel flat bottom hole acquired in the experiments of the present application;

FIG. 4 is a recovery result of a flat-bottom hole C-scan image obtained by a non-blind recovery method;

FIG. 5 is a graph of a flat bottom hole C scan image restored by the method of the present application;

FIG. 6 is a C-scan line for No. 1 flat bottom hole;

figure 7 is a C-scan line for No. 2 flat bottom hole.

Detailed Description

The technical scheme of the disclosure will be described in detail with reference to the accompanying drawings.

Fig. 1 is a flow chart of a method according to the present application, as shown in fig. 1, the method comprising:

step S1: and acquiring an ultrasonic C-scan image, and removing extreme value noise points on the ultrasonic C-scan image to obtain an ultrasonic image q to be restored.

Step S2: iterating the ultrasonic image f and the point spread function h to obtain the optimal point spread functionPoint spread function h_s 。

Specifically, the iterative solution formula of the ultrasound image f and the point spread function h is expressed as:

k represents the number of alternating iterations, and k is 0,1,2,3 …; f. of_kAnd h_kRespectively representing the ultrasonic image and the point spread function, f, obtained in the k-th iterative solution_k+1And h_k+1Respectively representing an ultrasonic image and a point spread function obtained in the (k + 1) th iterative solution; h is_sThe point spread function with the optimal parameters is represented, s represents a parameter vector of the ultrasonic imaging system, and g represents an ultrasonic image to be restored; equation (1) adopts RL-TV algorithm of total variation regularization constraint to solve f_k+1(ii) a The formula (2) is based on a multivariate Gaussian sound beam model MGB to obtain the parameter optimization expression of the point spread function, and h is carried out by adopting a particle swarm optimization algorithm_k+1And (4) carrying out optimization solution on the parameters.

Step S3: according to the point spread function h with the optimal parameters_sA corresponding restored ultrasound image is obtained.

Fig. 2 is a flowchart of an embodiment of the method of the present application, and it can be seen from fig. 2 that step S2 further includes:

s21: initialization f₀Initializing h as g₀(x,y)：

Wherein, (x, y) is the pixel coordinate in the ultrasonic image g to be restored.

S22: solving the intermediate ultrasonic image f by adopting RL-TV algorithm_k+1The iterative format of the RL-TV algorithm is:

wherein n represents the number of iterationsAnd n is 0,1,2,3 …, f_n(x, y) is the ultrasound image obtained in the nth iteration, f_n+1(x, y) is an ultrasonic image obtained by the (n + 1) th iteration, h (x, y) is a given point spread function, and h (-x, -y) represents a result obtained after the h (x, y) is turned along the x axis and the y axis; g (x, y) represents an ultrasound image to be restored; lambda [ alpha ]_TVFor regularization coefficients, div is the divergence calculation.

Step S23: and (4) performing alternate solving on the formula (1) and the formula (2) through the formula (4), judging whether the maximum iteration number is reached according to the alternate solving result, if the maximum iteration number is reached, turning to the step S25, and if not, turning to the step S24.

S24: solving intermediate point diffusion function h by adopting particle swarm optimization algorithm_k+1And go to step S22 to continue the iteration.

S25: and obtaining the point spread function with the optimal parameters and the corresponding restored ultrasonic image.

Step S24 further includes:

s241: for the sound field of water immersion ultrasonic detection, when the incident direction of ultrasonic waves is vertical to a water-solid interface, the sound field distribution in a solid workpieceComprises the following steps:

wherein the content of the first and second substances,representing particles (x) in a solid workpiece₂,y₂,z₂) ω ═ 2 pi f, which represents the ultrasonic probe circle frequency; j represents an imaginary number; v. of₀(ω) -surface particle vibration velocity of the ultrasound probe.

-transmission coefficient of gamma 2 wave in solid to p wave in water at liquid-solid interface; gamma 2-wave form in solid; p-longitudinal wave.

ρ₁Water density, p₂-a solid workpiece density;-the speed of longitudinal waves in water,-longitudinal acoustic velocity in solids;-a rayleigh distance;-the wave number of longitudinal waves in water;-number of longitudinal waves in the solid workpiece; a-probe radius; a. the_r,B_r-multivariate gaussian superposition complex coefficients; z is a radical of₁-water layer depth between the ultrasound probe and the workpiece surface.

S242: after obtaining the sound field distribution of the ultrasonic probe, calculating the inherent z of the sound field₀Point spread function at depth h (x, y):

whereinIs the sound field distribution obtained in step S341, v₀Representing the vibration velocity of particle on the surface of the ultrasonic probe, p representing the density of the solid workpiece,

s243: for the radius a of the ultrasonic probe and the center frequency f of the ultrasonic probe₀And depth z of investigation of sound field in workpiece₂Optimizing, comprising: let the characteristic parameter vector s ═ z₂,a,f₀]^TAccording to a characteristic parameter vector sThe point spread function with the optimal parameters determined in step S242 is denoted as h_s＝MGB(s)；

The optimization objective function when the particle swarm optimization algorithm is adopted for parameter optimization is as follows:

if the number of particle populations is I and the maximum number of iterations is T, the position of the I-th (I-1, 2,3, …, I) particle at the T (T-1, 2,3, …, T) iteration is:

the moving speed of the ith particle at the t iteration is:

the optimal position searched by the ith particle at present is as follows:

the optimal position currently searched by the whole particle population is as follows:

G^t＝(p_G1,p_G2,p_G3)； (11)；

the particles update speed and position according to equation (12):

where eta represents the inertia factor, lambda₁、λ₂Are self-learning factor and social learning factor, respectively, delta₁、δ₂Are all [0,1]Random numbers within a range.

S244: the parameter optimizing step of the point spread function based on the particle swarm optimization algorithm comprises the following steps:

s2441: setting the number I of particle populations, the maximum iteration number T and the initial positions of particlesInitial velocity of particlesParticle position boundary [ P ]_L,P_U]And particle velocity boundary [ V ]_L,V_U]。

S2442: calculating the fitness values of all the particles of the current iteration according to an objective function F(s)

S2443: comparing the current position fitness value of each particleA fitness value corresponding to the optimal position of the particleIf it is notThe optimal position of the particle is updated to the current position, otherwise, the optimal position of the particle is not changed.

S2444: comparing the current position fitness value of each particleFitness value F (G) of optimal position of particle swarm^t) If, ifThe optimal position of the particle swarm is updated toThe current position of the particle, otherwise, it is unchanged.

S2445: updating the velocity of individual particlesAnd position

S2446: judging whether the maximum iteration times is reached, if so, performing step S2447, otherwise, turning to step S2442;

s2447: outputting the current optimal result h_s。

As a specific example, fig. 3 is a C-scan image of a stainless steel flat bottom hole acquired in the experiments of the present application. The test block is made of 304 stainless steel, the diameters of No. 1 and No. 2 flat-bottom holes are 4mm and 6mm respectively, and the buried depth is 70 mm. Water immersion type ultrasonic C scanning imaging is adopted, the center frequency of an ultrasonic probe is 2.5MHz, the water range is 50mm, and the scanning step is 0.2 mm.

Fig. 4 shows the restoration result of the flat-bottom hole C-scan image obtained by the non-blind restoration method. The non-blind restoration method improves the size resolution of the original defect image, but because the accuracy of the initial point diffusion function is insufficient, the non-blind restoration cannot obtain an effective restoration result, and the original fuzzy and adhesive defects on the image are still adhered after restoration.

FIG. 5 is a graph of the recovery of a C-scan flat-bottom hole image by the method of the present application. The method and the device remarkably improve the definition of the original image and successfully separate two flat-bottom hole images.

Fig. 6 and 7 show the C-scan lines longitudinally through the center of the hole on two images of the flat-bottom hole, respectively, and the diameter size of each flat-bottom hole is obtained by the-6 dB method. Due to the non-blind restoration failure, there is no meaningful C-scan line in the restoration results of fig. 6 and 7, and accordingly no dimensional information for the two flat-bottom holes can be obtained. The diameters obtained on the original images of the No. 1 and No. 2 flat-bottom holes are respectively 19.6mm (relative error is 226.7%) and 24.8mm (relative error is 520.0%), the diameters obtained after the No. 1 and No. 2 flat-bottom holes are restored by the method are respectively 4.2mm (relative error is 30.0%) and 4.6mm (relative error is 15.0%), and the size accuracy is obviously improved.

The above description is only a preferred embodiment of the present invention, and is not intended to limit the present invention in any way, but any modifications or equivalent variations made according to the technical spirit of the present invention are within the scope of the present invention as claimed.