A post-processing method for particle field reconstruction using tomographic particle image velocimetry

Abstract

The present invention discloses a post-processing method for reconstructing a particle field by tomographic particle image velocimetry, comprising the following steps: (1) obtaining a reconstructed particle field of a flow field to be measured; (2) performing a first deconvolution operation on the reconstructed particle field to obtain a first predicted particle field; (3) performing an anisotropic convolution kernel estimation based on the first predicted particle field and the reconstructed particle field; (4) performing a deconvolution operation on the reconstructed particle field using the anisotropic convolution kernel obtained by the latest estimation to obtain a second predicted particle field; (5) performing an anisotropic convolution kernel estimation based on the second predicted particle field and the reconstructed particle field; (6) repeating steps (4) and (5) until an iteration stop condition is met; (7) performing a deconvolution operation on the reconstructed particle field using the anisotropic convolution kernel estimated in step (6) to obtain a repaired particle field. Compared with the prior art, the present invention has the characteristics of effectively restoring the spherical shape of the original particles and significantly suppressing ghost particles in the reconstructed particle field.

Description

Technical Field

[0001] The present invention belongs to the field of particle image velocimetry, and more particularly, relates to a post-processing method for reconstructing a particle field using tomographic particle image velocimetry. Background Art

[0002] Tomographic Particle Image Velocimetry (Tomo-PIV) is commonly used to measure three-dimensional, three-component flow fields. Tomo-PIV uses multiple cameras to simultaneously capture images of particles in a flow field and then analyzes the captured images to measure three-dimensional flow velocity. Its fundamental assumption is that well-traceable particles can reveal the structure of the flow field. Tracer particles are dispersed in the flow field, illuminated by a volume light source, and simultaneously recorded by cameras positioned at different angles. The three-dimensional particle intensity distribution is reconstructed tomographically from these recorded two-dimensional images using previously obtained calibration results. Finally, the two reconstructed particle fields are analyzed using a cross-correlation algorithm or a variational optical flow algorithm to obtain the particle displacement field. Based on Tomo-PIV's fundamental assumption, the particle displacement field can be assumed to be the displacement field of the fluid in which the particles reside. Dividing the obtained displacement field by the time interval between camera captures yields the fluid velocity. Tomo-PIV is now widely used in a variety of fields, including experimental fluid mechanics, biomedicine, and aerospace, as a powerful, full-field, instantaneous, and non-contact quantitative measurement method.

[0003] The accuracy of Tomo-PIV reconstruction significantly impacts the accuracy of subsequent velocity measurements. Tomo-PIV employs the classic MART algorithm for tomographic reconstruction of the particle field, which can lead to two errors: In the reconstructed particle field, particle elongation along the depth direction leads to uncertainty in the velocity measurement in the depth direction; Furthermore, the reconstructed ghost particles significantly affect the peak signal-to-noise ratio (PSNR) of the cross-correlation measurement and the grayscale conservation assumption of the optical flow measurement, resulting in a sharp decrease in the accuracy of the velocity measurement algorithm. Reducing the effects of particle elongation and ghost particles can be achieved by reducing the tracer particle concentration. However, a lack of sufficient tracer particles can result in the loss of fine flow field structure, reducing the spatial resolution of the velocity measurement.

[0004] Many researchers are currently studying the reconstruction characteristics of Tomo-PIV, seeking to improve the traditional MART algorithm to reduce particle elongation or suppress ghost particles in the reconstructed field, thereby achieving high-precision reconstruction results under conditions of high tracer particle concentrations. Classic improvements include reconstruction algorithms such as SFIT-MART (spatial filtering improved tomography MART) and IntE-MART (intensity enhanced MART). However, SFIT-MART primarily suppresses the intensity of ghost particles to improve reconstruction accuracy, without considering the errors introduced by the reconstructed particles. Furthermore, because SFIT-MART filters the particle field, suppressing ghost particles simultaneously weakens the intensity of real particles, resulting in limited accuracy improvement. IntE-MART's one-dimensional inverse diffusion equation only considers particle elongation along the depth direction, but not deformation in other directions. While IntE-MART can effectively restore the spherical shape of particles and suppress ghost particles, its accuracy needs further improvement.

[0005] In summary, particle elongation and ghost particles seriously reduce the accuracy of reconstructed particle fields, hindering the application of Tomo-PIV at high concentrations of tracer particles. It is necessary to improve the existing algorithm to address these two problems. Summary of the Invention

[0006] This invention addresses the problems of existing particle reconstruction methods, which cannot effectively restore particle shape and suppress ghost particles. It provides a post-processing method for reconstructing particle fields using tomographic particle image velocimetry. This method aims to estimate the anisotropic convolution kernel acting on the real particle field based on a given reconstructed particle field, and to design constraints for ghost particles. This method effectively suppresses ghost particles in the reconstructed particle field while deconvolving the reconstructed particle field to restore the original particle shape, resulting in a more accurate reconstruction. This method enables the application of Tomo-PIV at higher particle concentrations and improves the accuracy and resolution of the final measured velocity field.

[0007] To achieve the above objectives, in a first aspect, the present invention provides a post-processing method for reconstructing a particle field by tomographic particle image velocimetry, comprising the following steps:

[0008] (1) Obtaining the reconstructed particle field of the flow field to be measured;

[0009] (2) performing a first deconvolution operation on the reconstructed particle field to obtain a first predicted particle field;

[0010] (3) performing anisotropic convolution kernel estimation based on the first predicted particle field and the reconstructed particle field;

[0011] (4) performing a deconvolution operation on the reconstructed particle field using the anisotropic convolution kernel obtained by the latest estimation to obtain a second predicted particle field;

[0012] (5) performing anisotropic convolution kernel estimation based on the second predicted particle field and the reconstructed particle field;

[0013] (6) Repeat steps (4) and (5) until the iteration stop condition is met;

[0014] (7) Deconvolution operation is performed on the reconstructed particle field using the anisotropic convolution kernel estimated in step (6) to obtain a repaired particle field.

[0015] Furthermore, the step (1) includes:

[0016] (11) obtaining multiple sets of two-dimensional images and calibration data of the flow field to be measured at different angles where tracer particles are spread;

[0017] (12) Using the calibration data, multiple sets of two-dimensional images are subjected to tomographic reconstruction to obtain a reconstructed particle field.

[0018] Furthermore, in step (12), the MART algorithm is used for tomographic reconstruction, and its calculation formula is as follows:

[0019]

[0020] Where E is the intensity distribution of the reconstructed particle field; the superscript m is the number of iterations; the subscript j represents the voxel in the three-dimensional particle field; I is the intensity distribution of the two-dimensional image; the subscript i represents the pixel in the two-dimensional image; w i,j is the weight coefficient of the Gaussian distribution, which indicates the influence of the j-th voxel on the i-th pixel; N i are all pixels in the neighborhood of pixel i affected by voxel j; μ is the relaxation coefficient.

[0021] Furthermore, in the step (2), a first deconvolution operation is performed on the reconstructed particle field using a unit convolution kernel.

[0022] Furthermore, in steps (3) and (5), the anisotropic convolution kernel estimation is achieved by the following formula:

[0023]

[0024] Among them, k represents the anisotropic convolution kernel; u represents the predicted particle field; v represents the reconstructed particle field; * represents the convolution operation; To find the first-order derivative operator; Represents the square of the L2 norm; γ1, γ2, γ3 are the first, second, and third weight coefficients respectively. The first-order derivative operators are respectively used for the x, y, and z directions.

[0025] Furthermore, in step (4), the deconvolution operation on the reconstructed particle field is achieved by the following formula:

[0026]

[0027] Wherein, λ is the fourth weight coefficient.

[0028] Furthermore, in step (7), the deconvolution operation on the reconstructed particle field is achieved by the following formula:

[0029]

[0030] Among them, α and θ are the fifth and sixth weight coefficients respectively; ‖‖0 represents the L0 norm; ‖‖ TV represents the total variation norm; χ and 1-χ represent the probability of a voxel being a real particle and a ghost particle at different intensities, respectively.

[0031] Furthermore, λ is equal to 0.2, γ1 and γ2 are equal to 30, γ3 is equal to 200, α is between 2 and 3, and θ is between 10 and 20.

[0032] In a second aspect, the present invention provides a post-processing system for reconstructing a particle field by tomographic particle image velocimetry, characterized in that it comprises: a computer-readable storage medium and a processor;

[0033] The computer-readable storage medium is used to store executable instructions;

[0034] The processor is configured to read the executable instructions stored in the computer-readable storage medium and execute the post-processing method for reconstructing a particle field by tomographic particle image velocimetry as described in the first aspect.

[0035] In general, the above technical solutions conceived by the present invention can achieve the following technical effects compared with the prior art:

[0036] 1. Effectively restore the spherical shape of particles in the real particle field, which can reduce the uncertainty of velocity in the depth direction in subsequent velocity measurements, making the velocity measurement in the depth direction more accurate;

[0037] 2. It effectively suppresses ghost particles in the reconstructed particle field. Suppressing ghost particles is equivalent to denoising the particle field, which is beneficial to the peak signal-to-noise ratio of the cross-correlation algorithm in the subsequent velocity measurement algorithm, and makes the noise-sensitive variational optical flow algorithm more accurate.

[0038] 3. This method effectively suppresses ghost particles, making Tomo-PIV measurements possible at higher particle concentrations. Increasing particle concentration can reveal fine flow structures in the flow field and improve the resolution of the final velocity measurement. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Figure 1 This is a flow chart of a post-processing method for reconstructing a particle field by tomographic particle image velocimetry provided by the present invention;

[0040] Figure 2 Schematic diagram of a post-processing method for reconstructing a particle field by tomographic particle image velocimetry provided by the present invention;

[0041] Figure 3 is a schematic diagram of the arrangement of a measuring device provided in an embodiment of the present invention;

[0042] Figure 4 are images captured by four cameras provided in an embodiment of the present invention;

[0043] Figure 5 : The particle fields reconstructed by four different reconstruction methods provided in the embodiments of the present invention, wherein (a) adopts the MART algorithm, (b) adopts the SFIT-MART algorithm, (c) adopts the IntE-MART algorithm, and (d) adopts the method proposed in the present invention;

[0044] Figure 6 4 are intensity histograms of particle fields reconstructed using four different reconstruction methods provided in the embodiments of the present invention. DETAILED DESCRIPTION

[0045] In order to make the objectives, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely for the purpose of explaining the present invention and are not intended to limit the present invention. In addition, the technical features involved in the various embodiments of the present invention described below may be combined with each other as long as they do not conflict with each other.

[0046] In the present invention, the terms "first", "second", etc. (if any) in the present invention and the drawings are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence.

[0047] See Figure 1 , combined with Figure 2 The present invention provides a post-processing method for reconstructing a particle field by tomographic particle image velocimetry, comprising steps (1) to (7).

[0048] (1) Obtain the reconstructed particle field of the flow field to be measured.

[0049] In this embodiment, the reconstructed particle field of the flow field to be measured is obtained specifically by the following method:

[0050] (11) Particle image acquisition and camera calibration

[0051] In Tomo-PIV reconstruction, the camera position must first be arranged so that it is aligned at a certain angle to the flow field to be measured, where the tracer particles are scattered. The particle images are then captured at the same time by the synchronizer trigger signal. In the subsequent MART reconstruction algorithm, a weight matrix (i.e., the contribution of the physical spatial position of the voxel to the pixel position it affects) is required. At this time, each camera needs to be calibrated to obtain the camera's intrinsic and extrinsic parameter matrices and distortion parameters, which are expressed as:

[0052]

[0053] where f x ,f y ,u c ,v c are the camera internal parameters, representing the focal length in the x-direction, the focal length in the y-direction, the horizontal center position of the image, and the vertical center position of the image; γ is the image tangential distortion parameter; R 3×3 and T 3×1 are the camera external parameters, representing the camera rotation matrix and translation matrix respectively.

[0054] (12) Tomographic particle field reconstruction

[0055] After obtaining the particle images taken at multiple angles at the same time and the corresponding camera calibration parameters, the influence of the three-dimensional voxel j on each pixel i can be calculated as the weight matrix w i,j , the MART algorithm is used to iteratively calculate the voxel intensity field E(X j ,Y j ,Z j ) to update:

[0056]

[0057] After m iterations, the reconstructed particle field is obtained.

[0058] Where E is the intensity distribution of the reconstructed particle field; the superscript m is the number of iterations; the subscript j represents the voxel in the three-dimensional particle field; I is the intensity distribution of the two-dimensional image; the subscript i represents the pixel in the two-dimensional image; w i,j is the weight coefficient of the Gaussian distribution, which indicates the influence of the j-th voxel on the i-th pixel; N i are all pixels in the neighborhood of pixel i affected by voxel j; μ is the relaxation coefficient, which controls the stability of MART iteration.

[0059] (2) Performing a first deconvolution operation on the reconstructed particle field to obtain a first predicted particle field.

[0060] It is understandable that, because Tomo-PIV reconstruction is a view-limited reconstruction, the reconstructed particles are mainly elongated along the depth direction of the camera. This elongation effect can be regarded as the result of the interaction between real spherical particles and an anisotropic convolution kernel, which can be described mathematically as:

[0061] v=k*u+ε,

[0062] Among them, u is the real particle field, v is the reconstructed particle field with limited viewing angle, k is the anisotropic convolution kernel, and ε is random noise. In the above formula, the only known thing is the reconstructed particle field v. This is an underdetermined problem, and there are countless combinations that can satisfy the above formula. Therefore, in order to solve the real particle field u more accurately, we first need to estimate the anisotropic convolution kernel. Here, we can first predict a rough real particle field (that is, the first predicted particle field), and then estimate the anisotropic convolution kernel based on the difference between the rough real particle field and the reconstructed particle field. The effect of convolution often blurs edge information. Therefore, the edges of the rough real particle field that needs to be predicted first must be sharp, and a rough particle field is predicted based on the sharpened edges and the edge information blurred by the convolution kernel. This process can be described as:

[0063]

[0064] Where * represents the convolution operation, λ is the fourth weight coefficient, To find the first-order derivative operator.

[0065] It should be noted that, in step (2), the first deconvolution operation on the reconstructed particle field can be performed using an anisotropic convolution kernel or an isotropic convolution kernel (for example, a unit convolution kernel), and is not limited to the above-described process.

[0066] (3) Estimating anisotropic convolution kernel based on the first predicted particle field and the reconstructed particle field.

[0067] With a rough estimate of the real particle field, we can estimate the anisotropic convolution kernel. In order to satisfy the anisotropy of the estimated convolution kernel, we need to impose different constraints on different directions of convolution. From the reconstruction with limited viewing angle, it can be seen that the reconstructed particle field is mainly elongated along the depth direction. Therefore, the intensity information in the depth direction and perpendicular to the depth direction can be constrained. The final estimated convolution kernel has the largest intensity distribution in the depth direction, while the intensity distribution perpendicular to the depth direction is very small. Doing so can satisfy the requirement that during deconvolution, the intensity in the depth direction is suppressed to be greater than that in the direction perpendicular to the depth direction, so that the particles that are mainly elongated in the depth direction are restored to a spherical shape. This process can be described as:

[0068]

[0069] Among them, γ1, γ2, and γ3 are the first, second, and third weight coefficients respectively. represents the square of the L2 norm, The first-order derivative operators are respectively used for the x, y, and z directions.

[0070] (4) Deconvolution operation is performed on the reconstructed particle field using the anisotropic convolution kernel obtained by the latest estimation to obtain a second predicted particle field.

[0071] In step (4), the deconvolution operation of the reconstructed particle field is achieved by the following formula:

[0072]

[0073] The description of this process is similar to the specific description process given in step (2), the difference is that in step (4) it is clearly stated that the reconstructed particle field is deconvolved using the anisotropic convolution kernel obtained by the latest estimation.

[0074] (5) Estimating anisotropic convolution kernel based on the second predicted particle field and the reconstructed particle field.

[0075] In step (5), the process of anisotropic convolution kernel estimation is similar to step (3), the only difference is that the input predicted particle field is different, which will not be repeated here.

[0076] (6) Repeat steps (4) and (5) until the iteration stop condition is met.

[0077] It is understandable that the accuracy of the final estimated anisotropic convolution kernel can be improved through multiple iterations. In this embodiment, the iteration stopping condition can be reaching a preset number of iterations or reaching a preset anisotropic convolution kernel accuracy, for example, the difference between the anisotropic convolution kernel estimated in the last iteration and the anisotropic convolution kernel estimated in the previous iteration is within an allowable range.

[0078] (7) Deconvolution operation is performed on the reconstructed particle field using the anisotropic convolution kernel estimated in step (6) to obtain a repaired particle field.

[0079] In step (7), the formula given in step (4) can also be used to perform a deconvolution operation on the reconstructed particle field using the anisotropic convolution kernel obtained by the final estimation to obtain the repaired particle field.

[0080] However, considering that ghost particles are an inevitable source of error in Tomo-PIV reconstruction, they are also the main factor restricting Tomo-PIV high-precision velocity measurements. Therefore, after estimating an accurate anisotropic convolution kernel, while deconvolving the elongated particle field, constraints on ghost particles can be added to achieve the effect of restoring the particle shape while suppressing ghost particles. In the reconstructed particle field, ghost particles usually have lower intensity values. In the grayscale histogram of the particle field, ghost particles are more likely to be distributed in low-intensity areas, while real particles have higher intensity values ​​and are more likely to be distributed in high-intensity areas. By using this prior knowledge of the probability distribution of ghost particles and real particles, combined with the L0 norm to remove the number and the TV norm's edge-preserving denoising characteristics, these two norms can be combined to design a deconvolution constraint equation to effectively suppress ghost particles in the restored particle field while deconvolving the particle field. This process can be described as:

[0081]

[0082] Among them, α and θ are the fifth and sixth weight coefficients respectively; ‖‖0 represents the L0 norm; ‖‖ TV represents the total variation norm; χ and 1-χ represent the probability of a voxel being a real particle and a ghost particle at different intensities, respectively.

[0083] The following are examples:

[0084] Taking the jet measurement as an example, the jet experiment was carried out in a glass tank with a length of 2.5 meters, a width of 0.2 meters, and a height of 0.45 meters. The jet nozzle outlet diameter was 15 mm, the ambient temperature was 18.1 ° C, the pipe flow rate was 0.77 L / min, and the jet density was 998.576 kg / mm 3 , Reynolds number is 1032.83. The arrangement of the whole device is as follows Figure 3 Next, the implementation steps of the post-processing method for three-dimensional reconstruction of particle fields by tomographic particle image velocimetry of the invention are described in detail:

[0085] (1) Obtain the reconstructed particle field of the flow field to be measured.

[0086] In this embodiment, the reconstructed particle field of the flow field to be measured is obtained specifically by the following method:

[0087] (11) Particle image acquisition and camera calibration

[0088] The camera is arranged as Figure 3 As shown, four cameras are used to shoot the jet under test, where the camera resolution is 1920×1080 pixels. 2The four cameras are divided into two groups and placed on both sides of the measured flow field, with the angle between each group of cameras being 60°. In order to ensure that the four cameras are exposed at the same time, a synchronizer is used to set the trigger signal to control the four cameras to be exposed at time t to obtain particle images, such as Figure 4 As shown. Each camera takes images of the calibration plate at different positions, and then uses OpenCV to calibrate each camera, in preparation for the subsequent calculation of the weight matrix by the MART algorithm. The physical size of the flow field area measured is 120×80×24mm 3 .

[0089] (12) Tomographic particle field reconstruction

[0090] Combined with the particle images taken in (11) and the camera calibration parameters obtained, MATLAB is used to read the image pixel values, and a three-dimensional particle field with an initialization intensity of 1 is constructed for the MART iterative algorithm. After updating 6 times, the reconstructed particle field with large errors, particle elongation, and many ghost particles is obtained, as shown in the figure. Figure 2 As shown. The reconstructed particle field size is 1157×557×307 voxel 3 .

[0091] (2) Performing a first deconvolution operation on the reconstructed particle field to obtain a first predicted particle field.

[0092] For the reconstructed particle field v obtained by the MART algorithm in (1), a sharpened particle field u is predicted using the following equation:

[0093]

[0094] The first term in the formula is called the data fidelity term, which ensures the similarity between the reconstructed particle field and the potential sharp particle field convolved with the estimated kernel. The second term is called the regularization term, which penalizes the first-order gradient of the potential sharp particle field to ensure sharpness. In step (2), the unit convolution kernel is used (the central element is 1 and the rest are 0). Obviously, the above equation is a least squares problem, and the equation satisfied by the optimal solution can be obtained by solving the Euler-Lagrange equation:

[0095]

[0096] in, Denotes the conjugate of the convolution kernel k, and Δ is the Laplace operator. Assuming that all convolutions and gradients satisfy the loop boundary condition assumption, the famous convolution theorem can be used to quickly obtain the optimal solution of the above equation:

[0097]

[0098] in, represents the fast Fourier transform, represents the inverse Fourier transform, represents the complex conjugate, is the fast Fourier transform of the Laplace convolution kernel, and |·| represents the amplitude operation.

[0099] (3) Estimating anisotropic convolution kernel based on the first predicted particle field and the reconstructed particle field.

[0100] Once a rough particle field u is predicted, the anisotropic convolution kernel k can be estimated. As mentioned before, the anisotropic convolution kernel estimation can be modeled as minimizing the following variational problem:

[0101]

[0102] Assuming z represents the depth direction, γ3 must be greater than γ1 and γ2 to ensure that the gradient of the convolution kernel in the z direction is minimized when the above variational equation is minimized (that is, the convolution kernel is smoothest in the z direction). Since the particles remain close to their original size in the plane perpendicular to the depth direction, γ1 and γ2 are set to the same small value. When the total viewing angle between the cameras decreases, the effect of particle elongation becomes more significant, so γ3 should be increased. In addition, the size of the convolution kernel is a key parameter and must be set to be roughly equal to the size of the elongated particle to accurately estimate the strength of the convolution kernel. For simplicity, assume that the convolution kernel is a square with an odd size. Minimization of the above equation can be achieved by solving the Euler-Lagrange equation in the frequency domain:

[0103]

[0104]

[0105]

[0106]

[0107] Where, · represents the inner product, and the superscript T represents the matrix transpose. represents the difference operator in three directions. After solving for the estimated convolution kernel, the kernel size parameter is used to obtain only the required portion of the kernel. In addition, intensities less than 5% of the kernel's maximum intensity are removed, and the kernel is then renormalized to reduce noise. The resulting anisotropic convolution kernel is used to predict the particle field in subsequent iterations.

[0108] (4) Deconvolution operation is performed on the reconstructed particle field using the anisotropic convolution kernel obtained by the latest estimation to obtain a second predicted particle field.

[0109] In step (4), the deconvolution operation of the reconstructed particle field is achieved by the following formula:

[0110]

[0111] The description of this process is similar to the specific description process given in step (2), the difference is that in step (4) it is clearly stated that the reconstructed particle field is deconvolved using the anisotropic convolution kernel obtained by the latest estimation.

[0112] (5) Estimating anisotropic convolution kernel based on the second predicted particle field and the reconstructed particle field.

[0113] In step (5), the process of anisotropic convolution kernel estimation is similar to step (3), the only difference is that the input predicted particle field is different, which will not be repeated here.

[0114] (6) Repeat steps (4) and (5) until the iteration stop condition is met.

[0115] It is understandable that the accuracy of the final estimated anisotropic convolution kernel can be improved through multiple iterations. In this embodiment, the iteration stopping condition can be reaching a preset number of iterations or reaching a preset anisotropic convolution kernel accuracy, for example, the difference between the anisotropic convolution kernel estimated in the last iteration and the anisotropic convolution kernel estimated in the previous iteration is within an allowable range.

[0116] (7) Deconvolution operation is performed on the reconstructed particle field using the anisotropic convolution kernel estimated in step (6) to obtain a repaired particle field.

[0117] The accurate anisotropic convolution kernel estimated in step (6) can be used to deconvolve the reconstructed particle field. In addition to restoring the shape of the spherical particles, some valuable information can also be used to make the potential actual particle field more accurate. Due to the existence of ghost particles, the sparsity of the original particle field is destroyed. A noteworthy change is that the non-zero gradient of the reconstructed particle field increases significantly. Therefore, the L0 norm can be used to constrain the particle intensity gradient in the potential actual particle field to reduce the number of non-zero particles to meet the sparsity. In general, total variation (TV) regularization can effectively preserve edges. By constraining the potential actual particle field with the total variation norm, the sharp edges of the particles can be successfully preserved and weak noise can be eliminated (the total variation norm has the effect of edge preservation and denoising). Assuming that it is possible to know in advance which voxels belong to ghost particles and which voxels belong to real particles, they can be constrained by the L0 norm and the total variation norm respectively. In this case, most of the ghost particles can be removed while preserving the boundaries of the actual particles. In the proposed deconvolution model, the concept of probability is adopted to assign probabilities to voxels with different intensities based on the prior information that ghost particles are mainly distributed at lower intensities. When the intensity of voxels is low, the possibility that they are ghost particles increases; therefore, the gradients of these voxels should use the L0 norm to limit the number of ghost particles. When voxels have higher intensities, they are more likely to be actual particles. The total variation norm is then used to constrain the intensity to provide edge protection. We use the prior information of the probability density distribution curves of ghost particles and real particles to mark the probability of a voxel being a real particle at different intensities as χ, and 1-χ is the probability of the voxel being a ghost particle. The probability-based hybrid variational deconvolution model can be summarized as:

[0118]

[0119] The split Bregman algorithm is used to split the equation into three sub-problems u, g, and d, and solve them separately.

[0120] ①For the u subproblem, the corresponding constraint equation is:

[0121]

[0122] Among them, b is the split Bregman variable, which makes the solution more stable. The solution of this subproblem is:

[0123]

[0124] Where div represents the divergence.

[0125] ②For the g subproblem, the corresponding constraint equation is:

[0126]

[0127] The solution to this subproblem is:

[0128]

[0129] Among them, the subscript p represents each element.

[0130] ③For the d sub-problem, the corresponding constraint equation is:

[0131]

[0132] The solution to this subproblem is:

[0133]

[0134] Finally, the Bregman variable b is updated as follows:

[0135]

[0136] The above three sub-problems are iterated alternately to obtain the final solution, and the termination condition is:

[0137]

[0138] Where t represents the number of Bregman iterations, N represents the number of all voxels in the particle field, and TOL is the tolerance, which is set to 0.01.

[0139] After solving step (7), the ghost particles in the particle field can be effectively suppressed while deconvolving the reconstructed particle field to restore the original particle shape, significantly improving the accuracy of the reconstructed particle field. The final result is as follows: Figure 2 In order to illustrate the advantages of the present invention in particle field reconstruction, the jet experiment in the embodiment is carried out using the MART algorithm, the SFIT-MART algorithm and the IntE-MART algorithm as comparisons to compare the shape of the reconstructed particle field and the curve of the intensity histogram distribution. Figure 5 As shown in FIG, the particle shapes reconstructed by the four reconstruction algorithms are visualized. It can be seen that the particle field reconstructed by the MART algorithm is relatively fragmented and the particle shapes are irregular, while the method of the present invention ( Figure 5 The reconstructed particles in (d) are closest to spherical and each particle is more complete. Figure 6 In the figure, the particle field reconstructed by the method of the present invention (solid line) has the smallest probability density in the low-intensity area and the largest probability density in the high-intensity area, which also shows that the method of the present invention effectively suppresses ghost particles (ghost particles are mainly distributed in the low-intensity area) and enhances the real particles with high intensity. Figure 5 and Figure 6The results jointly demonstrate the effectiveness of the method of the present invention in both restoring the true particle shape and suppressing ghost particles.

[0140] It will be easily understood by those skilled in the art that the above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A post-processing method for reconstructing a particle field using tomographic particle image velocimetry, characterized in that: The following steps are involved: (1) Obtain the reconstructed particle field of the flow field to be measured; (2) performing a first deconvolution operation on the reconstructed particle field to obtain a first predicted particle field; (3) performing anisotropic convolution kernel estimation based on the first predicted particle field and the reconstructed particle field; (4) performing a deconvolution operation on the reconstructed particle field using the anisotropic convolution kernel obtained by the latest estimation to obtain a second predicted particle field; (5) performing anisotropic convolution kernel estimation based on the second predicted particle field and the reconstructed particle field; (6) Repeat steps (4) and (5) until the iteration stop condition is met; (7) performing a deconvolution operation on the reconstructed particle field using the anisotropic convolution kernel estimated in step (6) to obtain a repaired particle field; In steps (3) and (5), the anisotropic convolution kernel estimation is achieved by the following formula: Among them, k represents the anisotropic convolution kernel; u represents the predicted particle field; v represents the reconstructed particle field; Represents the convolution operation; To find the first-order derivative operator; represents the square of the L2 norm; are the first, second and third weight coefficients respectively, 、 、 The first-order derivative operators are respectively used for the x, y, and z directions; In step (4), the deconvolution operation of the reconstructed particle field is achieved by the following formula: in, is the fourth weight coefficient; In step (7), the deconvolution operation of the reconstructed particle field is achieved by the following formula: in, are the fifth and sixth weight coefficients respectively; represents the L0 norm; represents the total variation norm; and They represent the probability of a voxel being a real particle and a ghost particle at different intensities.

2. The post-processing method for reconstructing a particle field by tomographic particle image velocimetry according to claim 1, wherein: Step (1) includes: (11) Obtain multiple sets of two-dimensional images and calibration data of the flow field to be measured at different angles where tracer particles are spread; (12) Performing tomographic reconstruction on multiple sets of two-dimensional images using the calibration data to obtain a reconstructed particle field.

3. The post-processing method for reconstructing a particle field by tomographic particle image velocimetry according to claim 2, wherein in step (12), the tomographic reconstruction is performed using the MART algorithm, and the calculation formula is as follows: Where E is the intensity distribution of the reconstructed particle field; the superscript m is the number of iterations; the subscript j represents the voxel in the three-dimensional particle field; I is the intensity distribution of the two-dimensional image; the subscript i represents the pixel in the two-dimensional image; is the weight coefficient of the Gaussian distribution, which indicates the influence of the j-th voxel on the i-th pixel; are all pixels in the neighborhood of pixel i affected by voxel j; is the relaxation coefficient.

4. The post-processing method for reconstructing a particle field by tomographic particle image velocimetry according to claim 1, wherein: In step (2), a first deconvolution operation is performed on the reconstructed particle field using a unit convolution kernel.

5. The post-processing method for reconstructing a particle field by tomographic particle image velocimetry according to claim 1, wherein: is equal to 0.2, and is equal to 30, is equal to 200, Between 2 and 3, Between 10 and 20.

6. A post-processing system for particle field reconstruction using tomographic particle image velocimetry, characterized in that: include: Computer-readable storage medium and processor; The computer-readable storage medium is used to store executable instructions; The processor is configured to read the executable instructions stored in the computer-readable storage medium and execute the post-processing method for reconstructing a particle field by tomographic particle image velocimetry according to any one of claims 1 to 5.

Images