A method for testing characteristics of flying fragments based on image processing technology

Abstract

This invention proposes a method and system for testing fragment dispersion characteristics based on image processing technology. The method involves calculating the fragment distribution radius based on the warhead yield, deploying binocular cameras within a safe area around the fragment distribution radius, calibrating and setting the camera parameters, and acquiring a sequence of fragment flight trajectory images. The acquired trajectory images are then used for automatic fragment identification and coordinate extraction. The extracted fragments are matched against each other in the fields of view of the two cameras, and the three-dimensional coordinates of the fragments in space are calculated using the two-dimensional coordinates of the extracted fragments in the images. The fragment motion equations are fitted in three-dimensional space, and the fragment motion parameters are solved based on the fitted equations and the known three-dimensional coordinates of the fragments. This invention features convenient system layout, real-time recording, non-contact measurement, high measurement accuracy, and fast analysis speed, providing an effective theoretical basis for the calculation and analysis of damage to new targets.

Description

Technical Field

[0001] This invention belongs to the field of image processing technology, specifically relating to a method and system for testing the fragment scattering characteristics based on image processing technology. Background Technology

[0002] In an explosion field, parameters such as the velocity, shape, and distribution of weapon warhead fragments are crucial for evaluating their destructive effectiveness. Because the types of warhead fragments, interfering factors, testing environments, and testing parameters vary, higher demands are placed on the testing principles and methods for fragments.

[0003] Currently, the main methods for fragment velocity measurement at test ranges include cutoff velocity measurement, radar velocity measurement, and high-speed photography. Cutoff velocity measurement is a commonly used method in domestic test ranges, primarily employing a fixed-distance, time-based measurement approach. However, due to the time-consuming and labor-intensive deployment and data acquisition of this method, domestic fragment testing is currently developing towards non-contact measurement methods such as photoelectric and electromagnetic methods. Radar velocity measurement, due to the complexity of radar echo analysis, is only suitable for relatively dispersed, low-speed flying fragments. Utilizing a high-speed camera as a photoelectric acquisition device allows for the measurement of the spatial position of fragment targets based on their motion images. This method offers advantages such as convenient system layout, real-time recording, non-contact measurement, high measurement accuracy, and fast analysis speed. Therefore, this application proposes a fragment dispersion characteristic testing method based on image processing. Summary of the Invention

[0004] This invention provides a method and system for testing the fragment scattering characteristics based on image processing technology.

[0005] The technical solution for achieving the present invention is: a method for testing the fragment scattering characteristics based on image processing technology, comprising the following steps:

[0006] Step 1: Calculate the fragmentation radius based on the warhead yield, deploy binocular cameras within the safe zone of the fragmentation radius, calibrate and set the parameters of the cameras, and collect a sequence of images of the fragmentation flight trajectory;

[0007] Step 2: Automatically identify and extract the coordinates of fragments from the acquired trajectory sequence images;

[0008] Step 3: Match the extracted fragments in the fields of view of the two cameras, and calculate the three-dimensional coordinates of the fragments in space using the two-dimensional coordinates of the extracted fragments in the image.

[0009] Step 4: Fit the fragment motion equation in three-dimensional space, and solve for the fragment motion parameters based on the fitted equation and the known three-dimensional spatial coordinates of the fragment.

[0010] Furthermore, binocular cameras are deployed within a safe area according to the fragment distribution radius, and the cameras are calibrated and their parameters are set, including the following steps:

[0011] Step S11: Use two high-speed cameras of the same model, both of which meet the high frame rate requirements of the test. Arrange the binocular cameras in a safe area according to the fragment distribution radius, and make the fields of view of the two cameras intersect in the detection area.

[0012] Step S12: In the experiment, a checkerboard image is used for calibration. Two cameras simultaneously capture images of the same checkerboard and solve the transformation relationship between the coordinate systems of the two cameras. The three-dimensional spatial position parameters of any point on the checkerboard are obtained by GPS and total station, and the position of the arbitrary point in the two-dimensional image is obtained by feature point extraction algorithm.

[0013] Step S13: Connect the two cameras to the GPS timing device respectively, and synchronize the two devices to GPS time using IRIG-B code. At the same time, set the cameras to B code synchronization mode. After the acquisition is completed, only the time of the first image of the two image sequences needs to be aligned to ensure time synchronization. If there is a frame loss phenomenon, and the video information acquired by one camera cannot find the corresponding first image of the other, then take frames according to the time sequence and take the time when the corresponding images at the same time are stored as the start to ensure time synchronization.

[0014] Step S14: Obtain the trajectory sequence image of the fragment flight.

[0015] Furthermore, the acquired image sequence is subjected to automatic fragment identification and coordinate extraction, including the following steps:

[0016] Step S21: Preprocess and denoise the fragment flight trajectory sequence images acquired by the binocular camera. Based on the distribution of wavelet coefficients, different frequency bands are processed separately. A soft threshold is set. When the absolute value of a wavelet coefficient is less than the given threshold, it is set to 0; when it is greater than the threshold, it is subtracted from the threshold. The basic principle is as follows:

[0017]

[0018] In equation (1), sign(w) is the sign function, w is the wavelet coefficient of the original image, and w λ These are the wavelet coefficients regenerated after wavelet transform, and λ is the threshold.

[0019] Step S22, let f k f k+1 f k+2 Given three adjacent images, perform difference operations between the three images to obtain:

[0020]

[0021] In formula (2), S1, S2, and S3 are the target foreground gray matrices obtained by difference;

[0022] In step S23, by superimposing the foreground after difference, the gray value at the position of the target can be increased, and the gray difference between the target and the noise can be improved. By superimposing the foreground gray matrices after difference, we get:

[0023] D k = S1 + S2 + S3 (3)

[0024] In formula (3), D k is the superimposed matrix of the foreground gray matrices after difference;

[0025] In step S24, in order to further increase the contrast between the target foreground and the noise, square the gray value of each pixel point in D k to obtain:

[0026]

[0027] In formula (4), g k is the matrix obtained by squaring the gray values of the superimposed matrix;

[0028] In step S25, adopt an improved local threshold method to separate the target foreground from the noise. First, divide the image obtained by squaring the gray values into sub-images, and then perform the following operations on each sub-image:

[0029] Obtain the histogram of the sub-image f(i, j);

[0030] Traverse the histogram from the smallest to the largest gray level, find the first peak and the last peak, and record the corresponding gray values as T0 and T1 respectively. Denote Tmin and Tmax as the minimum and maximum gray levels of the sub-image respectively;

[0031] Let f(i, j) be the gray value of the pixel point in the sub-image. If f(i, j) > T0 and f(i, j) < T1, then use

[0032] calculate the new gray value to replace the gray value of the sub-image, where T is the threshold determined by the maximum between-class variance method for this sub-image;

[0033] Stitch all the sub-images processed as above into the original size to obtain an image with the background and the target separated;

[0034] In step S26, perform morphological processing to remove the isolated noise points and repair the holes that appear in the difference and segmentation processes;

[0035] Step S27: Count the number of fragments in the image and infer the shape of the fragments based on the shape of the fragments appearing in the image when flying in each trajectory.

[0036] Furthermore, the extracted fragments are matched against each other in the fields of view of the two cameras, and the three-dimensional coordinates of the fragments are calculated using the extracted two-dimensional coordinates, including the following steps:

[0037] Step S31: Extract corner points from the obtained image where the background and target are separated using the Harris corner extraction algorithm and perform initial matching of the corner points using the ZNCC strategy;

[0038] Step S32: Perform SIFT feature point detection and matching. Use the basic matrix calculated from SIFT feature points to remove mismatches from the SIFT matching points and initial matching corner points in the result of step S31.

[0039] Step S33: For the matching image after removing mismatches, calculate the fundamental matrix from the matching pairs using the eight-point method, and obtain the essential matrix from the fundamental matrix;

[0040] Step S34: Decompose the essential matrix to obtain the camera extrinsic parameters;

[0041] Step S35: Perform 3D spatial point triangulation on the optimized SIFT feature points and corner points;

[0042] Step S36: Obtain the intrinsic parameters A of the left and right cameras through calibration. l With A r The rotation matrix R between the two cameras, based on the principles of monocular imaging and coordinate transformation, has the following equation:

[0043]

[0044] In equation (5), x w y w , z w Let u1, v1 and u2, v2 be the world coordinates of the fragment, u1, v1 and u2, v2 be the pixel coordinates of the fragment in the left and right cameras respectively, s1, s2 be the working distance parameters of the left and right cameras respectively, and t be the pixel coordinates of the fragment in the left and right cameras respectively. l , t r R represents the offset vectors for the left and right cameras, respectively, transforming them from the world coordinate system to the camera coordinate system. l ,R r These are the rotation matrices for transforming the left and right cameras from the world coordinate system to the camera coordinate system, respectively.

[0045] Assuming the left camera coordinate system is the measurement coordinate system, i.e., the world coordinate system, the rotation matrix from the left camera coordinate system to the world coordinate system is R. l =E, translation vector is t l=0, meaning the left camera coordinate system coincides with the world coordinate system, and the translation vector is 0. Simplifying the camera intrinsic and extrinsic parameter matrices yields the following equation:

[0046]

[0047] The least squares method is used to calculate the coordinates of a 3D point, based on the following assumptions: the image points p1 and p2 of any point P in space on the two cameras C1 and C2 have been detected from the two images respectively, that is, p1 and p2 are known to be the corresponding points of the same point P in space; cameras C1 and C2 have been calibrated, and their corresponding projection matrices are M1 and M2 respectively, thus:

[0048]

[0049] Where (u1,v1,1) and (u2,v2,1) are the homogeneous coordinates of points P1 and P2 in their respective images; (X,Y,Z,1) are the homogeneous coordinates of point P in the world coordinate system; m ij k For M k The element in the i-th row and j-th column, eliminating s1 and s2, yields information about x. w y w , z w The coordinates of a three-dimensional point P can be obtained by using the four linear equations.

[0050] Furthermore, the fragment motion parameters are solved based on the fitted equation and the known three-dimensional spatial coordinates of the fragment, including the following steps:

[0051] Step S41, construct the differential equation of motion of the fragment in air:

[0052]

[0053] In equation (8), m is the fragment mass, and C x ρ is the frontal drag coefficient of the fragment, s is the local air density, T is the frontal area of ​​the fragment perpendicular to the flight direction, v is the flight time of the fragment, and v is the instantaneous velocity of the fragment.

[0054] Step S42: Obtain the coordinates of the projectile before static explosion and the time of explosion, and combine them with the obtained fragment trajectory to solve the problem using an indirect calculation method;

[0055] First, taking the first moment of the trajectory as the relative moment and its position as the relative starting position, we have:

[0056]

[0057] In the formula: T1 is the time corresponding to the first data point of the fragment trajectory, X′ is the distance between the fragment and the explosion point at time T1, v1 corresponds to the velocity of the fragment at time T1, v0 corresponds to the velocity of the fragment at time T0, x is the distance of the fragment from the explosion point at a certain time, and α is the velocity attenuation coefficient. Then we have:

[0058]

[0059] In equation (10), x(T) is the flight distance of the fragment at any time T, which is equivalent to:

[0060]

[0061] In equation (11): s tr Let t be the distance of the fragment relative to the coordinates at time T1. r The velocity v1 and velocity decay coefficient α can be obtained by fitting the obtained fragment trajectory data relative to the time difference at time T1.

[0062] Calculate the perpendicular coordinates from the warhead's center of mass to the fitted straight line of the fragment flight path. Since the fitted straight line of the fragment may not pass through the warhead's center of mass, the perpendicular coordinates are used as the starting point of the fragment's motion. Based on this, the distance between the first data point of the trajectory and the perpendicular coordinate is calculated as X′.

[0063] Substituting the velocity v1, velocity attenuation coefficient α, and distance X′ into equation (9), the initial velocity v0 is obtained.

[0064] A fragment scattering characteristic testing system based on image processing technology is provided, which realizes fragment scattering characteristic testing based on image processing technology through the fragment scattering characteristic testing method described above.

[0065] A computer device includes a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it performs a fragment scattering characteristic test based on image processing technology using the fragment scattering characteristic test method described above.

[0066] A computer-readable storage medium having a computer program stored thereon, wherein when the computer program is executed by a processor, the fragment scattering characteristic test method based on the image processing technology is implemented to achieve fragment scattering characteristic test.

[0067] Compared with the prior art, the significant advantages of this invention are: 1) It can solve the fragment motion and scattering characteristics in three-dimensional space more accurately; 2) The non-contact measurement system based on the principle of binocular vision is simpler to set up than the traditional target plate method and the data acquisition is faster; 3) The local maximum inter-class variance method is used when extracting fragment targets, which can effectively avoid the loss of target extraction caused by uneven lighting at the work site. Attached Figure Description

[0068] Figure 1 This is a system flowchart of the present invention.

[0069] Figure 2 This is a layout diagram of the binocular camera of the present invention.

[0070] Figure 3 The flowchart for the three-frame difference method with squared enhancement. Detailed Implementation

[0071] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific implementation processes described herein are merely illustrative and not intended to limit the scope of this application.

[0072] like Figure 1 As shown, a method for determining fragment scattering parameters based on image processing technology includes the following steps:

[0073] Step 1: Calculate the fragmentation radius based on the warhead yield, and deploy the binocular cameras within the safe zone according to the fragmentation radius. Figure 2 As shown, the camera is calibrated and its parameters are set to acquire images of the fragments passing through the intersecting fields of view of the two cameras during flight. The specific steps are as follows:

[0074] Based on the calculated fragmentation radius using the warhead yield, binocular cameras were positioned within a safe zone corresponding to this radius. Considering the initial fragment velocity of 800m / s-1500m / s, a high-speed camera with a shooting frequency of 1000s⁻¹ was selected. The lens focal length was approximately 165mm, the image resolution was 1280×800 pixels, and the pixel size was 20μm. The combined lateral and longitudinal field of view of the two cameras was approximately 10m, and the combined longitudinal field of view was approximately 5m. The spatial resolution per pixel was approximately 7mm. Both cameras were positioned 70m from the target, with a meeting angle of 60°.

[0075] Two high-speed cameras are connected to GPS timing devices, and the two devices are synchronized to GPS time using IRIG-B code. Simultaneously, the cameras are set to B-code synchronization mode. After data acquisition, only the time of the first image in the two image sequences needs to be aligned. If frame loss occurs, and the video information acquired by one camera cannot find the corresponding first image of the other, frames are retrieved sequentially from the time sequence, taking the moment when corresponding images at the same time are all stored as the starting point to ensure time synchronization.

[0076] The binocular camera can be triggered manually or by sound / vibration. Manual triggering requires connecting the high-speed camera to the control computer via a network cable to establish a local area network, and then manually controlling the triggering process through the computer.

[0077] Step 2: Automatically identify and extract coordinates of fragments from the collected sequence images;

[0078] Preprocessing and denoising of fragmented images acquired by a stereo camera can be performed by designing filters to process different frequency bands based on the distribution of wavelet coefficients. A soft threshold is set; when the absolute value of a wavelet coefficient is less than the given threshold, it is set to 0; when it is greater than the threshold, it is subtracted from the threshold. The basic principle is as follows:

[0079]

[0080] In equation (1), sign(w) is the sign function, w is the wavelet coefficient of the original image, and w λ These are the wavelet coefficients regenerated after wavelet transform, and λ is the threshold.

[0081] Let f k f k+1 f k+2 Given three adjacent images, we can obtain the following by performing difference operations on each of the three images:

[0082]

[0083] In equation (2), S1, S2, and S3 are the target foreground grayscale matrices obtained by difference.

[0084] Since noise is random, its magnitude varies at different locations and even at different times. However, the position of the fragment target within the three adjacent frames is relatively fixed. Superimposing the foreground after difference can increase the grayscale value of the target's location, thereby increasing the grayscale difference between the target and the noise. Superimposing the foreground grayscale matrices after difference yields:

[0085] D k =S1+S2+S3 (3)

[0086] In equation (3), D kis the superimposed matrix of the foreground gray level matrix after difference;

[0087] Step S24. To further increase the contrast between the target foreground and the noise, square the gray level value of each pixel in D k to obtain:

[0088]

[0089] In Equation (4), g k is the matrix obtained by squaring the gray level value of the superimposed matrix;

[0090] The contrast between the target foreground and the noise is enhanced. A threshold needs to be set to separate the target foreground from the noise. Since simple binary thresholds cannot meet the requirements of separating fragments from the background when the image illumination is uneven, there are sudden noises, or the background gray level changes greatly, an improved local threshold method is adopted. The local threshold process is as follows:

[0091] The images captured by two CCD cameras are of the same size. The image obtained by squaring the gray level value and the original image are both 1280*800. According to the size, the image obtained by squaring the gray level value is divided into sub-images of 128*80, and each gray square image is divided into 100 sub-images; the following operations are performed on each sub-image:

[0092] Obtain the histogram of the sub-image f(i, j);

[0093] Traverse the histogram from the smallest to the largest gray level, find the first peak and the last peak, and record the corresponding gray level values as T0 and T1 respectively. Denote Tmin and Tmax as the minimum and maximum gray levels of the sub-image respectively.

[0094] Let f(i, j) be the gray level value of the pixel in the sub-image. If f(i, j)>T0 and f(i, j)<T1, then use

[0095] to calculate a new gray level value to replace the gray level value of the sub-image, where T is the threshold determined by the maximum between-class variance method for this sub-image.

[0096] Stitch all the sub-images processed as above into the original size to obtain an image with the background and the target separated.

[0097] Furthermore, the between-class variance is defined as:

[0098] σ 2 =θ(μ - μ1) 2 +(1 - θ)(μ2 - μ) 2 =θ(1 - θ)(μ1 - μ2) 2

[0099] In the formula, For class C1, the gray values ​​i ∈ [0, m], and for class C2, the gray values ​​i ∈ [m+1, L-1]. μ represents the overall mean of the image, μ1 represents the class mean of C1, μ2 represents the class mean of C2, and θ is the probability of a gray value appearing in C1. Let σ 2 The largest threshold m is the optimal threshold.

[0100] Morphological processing removes isolated noise and repairs holes that occur during the difference and segmentation processes.

[0101] The above steps can extract the fragments from the background of the acquired image, thereby obtaining the two-dimensional coordinate information of the fragments in each frame of the image, counting the number of fragments in the image, and inferring and describing the shape of the fragments based on the shape of the fragments appearing in the image when flying in each trajectory.

[0102] Step 3: Match the extracted fragments in the fields of view of the two cameras, and calculate the three-dimensional coordinates of the fragments using the extracted two-dimensional coordinates;

[0103] The Harris corner extraction algorithm is used to extract corners and the ZNCC strategy is used to perform initial matching of corners. The basis matrix calculated from SIFT feature points is used to remove mismatches of SIFT matching points and initial matching corners.

[0104] Perform SIFT feature point detection and matching;

[0105] The fundamental matrix is ​​calculated from the matching pairs using the eight-point method, and the essential matrix can be obtained from the fundamental matrix.

[0106] The camera's extrinsic parameters can be obtained by decomposing the essential matrix;

[0107] Perform 3D spatial point triangulation on the optimized SIFT feature points and corner points;

[0108] Internal parameter A was obtained through calibration. l With A r Let R be the rotation matrix between the two cameras, and t be the translation vector. Based on the principles of monocular imaging and coordinate transformation, we have the following equation:

[0109]

[0110] In equation (5), x w y w , z w Let u1, v1 and u2, v2 be the world coordinates of the fragment, u1, v1 and u2, v2 be the pixel coordinates of the fragment in the left and right cameras respectively, s1, s2 be the working distance parameters of the left and right cameras respectively, and t be the pixel coordinates of the fragment in the left and right cameras respectively. l , t r R represents the offset vectors for the left and right cameras, respectively, transforming them from the world coordinate system to the camera coordinate system.l ,R r These are the rotation matrices for transforming the left and right cameras from the world coordinate system to the camera coordinate system, respectively.

[0111] Assuming the left camera coordinate system is the measurement coordinate system, i.e., the world coordinate system, the rotation matrix from the left camera coordinate system to the world coordinate system is R. l =E, translation vector is t l =0, meaning the left camera coordinate system coincides with the world coordinate system, and the translation vector is 0. Simplifying the camera intrinsic and extrinsic parameter matrices yields the following equation:

[0112]

[0113] To calculate the coordinates of a 3D point using the least squares method, we need to clarify the following basic assumptions: The image points p1 and p2 of any point P in space on the two cameras C1 and C2 have already been detected from the two images. That is, p1 and p2 are known to be corresponding points of the same point P in space. How these corresponding points are found is a matter of image registration. Another assumption is that cameras C1 and C2 are calibrated, and their projection matrices are M1 and M2, respectively. Therefore:

[0114]

[0115] Where (u1,v1,1) and (u2,v2,1) are the homogeneous coordinates of points P1 and P2 in their respective images; (X,Y,Z,1) are the homogeneous coordinates of point P in the world coordinate system; m ij k For M k The element in the i-th row and j-th column. Eliminating s1 and s2 from Equation 7 yields the information about x. w y w , z w The four linear equations are then used to determine the coordinates of the three-dimensional point P.

[0116] Step four: Simulate the trajectory of the fragments based on the image sequence and frame number, and count the number and location of the fragments.

[0117] Under certain assumptions, construct the differential equation of motion for the fragment in air. The assumptions are as follows:

[0118] (1) Ignoring the effect of gravity, the fragments only decelerate under the action of air resistance;

[0119] (2) Ignore the effect of fragment rotation, assume that the air resistance coefficient does not change, and that the air resistance acts on the centroid of the fragment;

[0120] (3) Ignore the effect of air temperature change on air density during static explosion.

[0121] The fragment decelerates under the influence of air resistance, and its differential equation of motion is:

[0122]

[0123] In equation (8), m is the fragment mass, and C x ρ is the frontal drag coefficient of the fragment, s is the local air density, T is the frontal area of ​​the fragment perpendicular to the flight direction, v is the flight time of the fragment, and v is the instantaneous velocity of the fragment.

[0124] The coordinates of the projectile before static explosion and the moment of explosion are obtained. Combined with the obtained fragment trajectory, an indirect calculation method is used to solve the problem.

[0125] First, taking the first moment of the trajectory as the relative moment and its position as the relative starting position, we have:

[0126]

[0127] In the formula: T1 is the time corresponding to the first data point of the fragment trajectory, X′ is the distance between the fragment and the explosion point at time T1, v1 corresponds to the velocity of the fragment at time T1, v0 corresponds to the velocity of the fragment at time T0, x is the distance of the fragment from the explosion point at a certain time, and α is the velocity attenuation coefficient. Then we have:

[0128]

[0129] In equation (10), x(T) is the flight distance of the fragment at any time T, which is equivalent to:

[0130]

[0131] In equation (11): s tr Let t be the distance of the fragment relative to the coordinates at time T1. r The velocity v1 and velocity decay coefficient α can be obtained by fitting the obtained fragment trajectory data relative to the time difference at time T1.

[0132] Calculate the perpendicular coordinates from the warhead's center of mass to the fitted straight line of the fragment flight path. Since the fitted straight line of the fragment may not pass through the warhead's center of mass, the perpendicular coordinates are used as the starting point of the fragment's motion. Based on this, the distance between the first data point of the trajectory and the perpendicular coordinate is calculated as X′.

[0133] Substituting the velocity v1, velocity attenuation coefficient α, and distance X′ into equation (9), the initial velocity v0 is obtained.

[0134] In summary, this invention can measure the fragment dispersion characteristics with high precision, and features convenient and simplified layout, high measurement accuracy, and fast analysis speed, providing an effective theoretical basis for the calculation and analysis of damage to new targets.

[0135] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0136] The embodiments described above are merely illustrative of several implementation methods of this application, and while the descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the invention patent. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of this application, and these all fall within the protection scope of this application. Therefore, the protection scope of this patent application should be determined by the appended claims.

Claims

1. A method for testing the flying characteristics of fragments based on image processing technology, characterized in that, The method comprises the following steps: Step 1: calculating the radius of fragment distribution according to the equivalent of the warhead, arranging binocular cameras in the safe area according to the radius of fragment distribution, calibrating and setting parameters of the cameras, and collecting trajectory sequence images of the fragments flying; Step 2: automatically identifying and extracting coordinates of the fragments from the collected trajectory sequence images, comprising the following steps: Step S21, pre-processing and denoising the trajectory sequence images of the fragments flying collected by the binocular cameras, processing different frequency bands respectively according to the distribution of wavelet coefficients, setting a threshold as a soft threshold, and setting the wavelet coefficient as 0 when the absolute value of the wavelet coefficient is less than the given threshold, and setting the wavelet coefficient as the threshold minus the threshold when the absolute value of the wavelet coefficient is greater than the threshold, and the basic principle is as follows: In formula (1), sign (w) is a sign function, w is a wavelet coefficient of an original image, w λ is a wavelet coefficient regenerated after wavelet transform, and λ is a threshold value. Step S22, set f k , f k+1 , f k+2 For three adjacent images, difference between three images is obtained respectively: In formula (2), S1, S2 and S3 are target foreground gray matrices obtained through difference; Step S23, superimposing the foreground after difference, so as to increase the gray value of the position where the target is located, improve the gray difference between the target and noise, superimposing the foreground gray matrix after difference, and obtaining D k = S1 + S2 + S3 (3) In formula (3), D k is a superposition matrix of the difference post-foreground grayscale matrix; Step S24, in order to further increase the contrast of target foreground and noise, square the gray value of D k Each pixel point is subjected to gray value square processing to obtain: In formula (4), g k is a gray value square processing matrix of the superposition matrix; Step S25, using an improved local threshold method to separate the target foreground from the noise, dividing the image obtained by squaring the gray value into sub-images, and then performing the following operations on each sub-image: Obtaining a histogram of the sub-image f(i,j); Traversing the histogram from small to large according to the gray level, finding the first peak and the last peak, and recording the corresponding gray values as T0 and T1 respectively, and recording Tmin and Tmax as the minimum and maximum gray of the sub-image respectively; If f(i,j)>T0 and f(i,j)<T1, then the new gray value is calculated to replace the gray value of the sub-image, wherein T is a threshold value determined by the maximum inter-class variance method for the sub-image; Splicing all the sub-images processed according to the above method into the original size to obtain an image in which the background and the target are separated; Step S26, performing morphological processing to remove isolated noise points and repair holes appearing in the difference and segmentation processing; Step S27, counting the number of fragments in the image, and making an inference description of the shape of the fragments according to the shape of the fragments appearing in the image when the fragments fly in each trajectory; Step 3: matching the extracted fragments in the fields of view of the two cameras with each other, and calculating the three-dimensional coordinates of the fragments in space by using the two-dimensional coordinates of the extracted fragments in the image; Step 4: fitting the motion equation of the fragments in the three-dimensional space, and solving the motion parameters of the fragments according to the fitted equation and the known three-dimensional coordinates of the fragments. Arranging binocular cameras in the safe area according to the radius of fragment distribution, calibrating and setting parameters of the cameras, comprising the following steps:

2. The method of claim 1, wherein, Step S11, using two high-speed cameras, and arranging binocular cameras in the safe area according to the radius of fragment distribution, and making the fields of view of the two cameras intersect in the detection area. ​ Step S12, in the test, a chessboard image is used for calibration, two cameras simultaneously capture the same chessboard, and the conversion relationship between the two camera coordinate systems is solved, wherein the three-dimensional space position parameters of any point on the chessboard are obtained through GPS and total station, and the position of the any point in the two-dimensional image is obtained through a feature point extraction algorithm; Step S13, a GPS time device is connected to each of the two cameras, the two devices are synchronized to GPS time through IRIG-B code, at the same time, the camera is set to B code synchronization mode, after the collection is completed, only the time of the first image of the two image sequences needs to be aligned to ensure the synchronization of the time, if there is a frame loss phenomenon, the video information collected by one camera cannot find the corresponding first image of the other camera, then the corresponding image stored at the same time is taken as the start time according to the time sequence backwardly, and the time synchronization is ensured; Step S14, a trajectory sequence image of the flying fragments is obtained.

3. The method of claim 1, wherein the method is characterized by, The extracted fragments are matched with each other in the fields of view of the two cameras, and the three-dimensional coordinates of the fragments are calculated based on the extracted two-dimensional coordinates, including the following steps: Step S31, the obtained image in which the background and the target are separated is used to extract corner points through a Harris corner point extraction algorithm, and the corner points are initially matched through a ZNCC strategy; Step S32, SIFT feature point detection and matching are performed, the fundamental matrix calculated based on the SIFT feature points is used to remove the mismatched points from the results of the initial matching of the corner points in step S31; Step S33, after the mismatched points are removed, the fundamental matrix is calculated based on the matching pairs through an eight-point method, and the essential matrix is obtained based on the fundamental matrix; Step S34, the essential matrix is decomposed to obtain the camera external parameters; Step S35, the optimized SIFT feature points and the corner points are triangulated in the three-dimensional space; Step S36, the internal parameters A of the left and right cameras are obtained by calibration l With A r The rotation matrix R between the two cameras has the following equation according to the monocular imaging principle and the coordinate transformation principle: In formula (5), x w， y w， z w are coordinates of the fragment in the world coordinate system, u1, v1 and u2, v2 are respectively coordinates of the fragment in the pixel coordinate system of the left and right cameras, s1, s2 are respectively working distance parameters of the left and right cameras, t l , t r are respectively offset vectors of the left and right cameras from the world coordinate system to the camera coordinate system, R l , R r are respectively rotation matrices of the left and right cameras from the world coordinate system to the camera coordinate system. Assuming that the left camera coordinate system is the measurement coordinate system, i.e. the position of the world coordinate system, the rotation matrix of the left camera coordinate system to the world coordinate system is R l = E, and the translation vector is t l = 0, i.e. the left camera coordinate system and the world coordinate system coincide, and the translation vector is 0. The camera intrinsic matrix and the extrinsic matrix are simplified to obtain the following formula: The three-dimensional point coordinates are calculated through the least square method based on the following assumptions: the image points p1 and p2 of any point P in the space on the two cameras C1 and C2 have been detected from the two images, that is, p1 and p2 are the corresponding points of the same point P in the space; the cameras C1 and C2 have been calibrated, and the corresponding projection matrices are M1 and M2, and thus: where (u1, v1, 1) and (u2, v2, 1) are the image homogeneous coordinates of P1 and P2 in their respective images; (X, Y, Z, 1) is the homogeneous coordinates of P in the world coordinate system; m ij k is the (i, j)th element of M k , and eliminating s1 and s2 gives four linear equations in x w， y w， z w , from which the coordinates of the three-dimensional point P can be solved.

4. The method of claim 1, wherein, The motion parameters of the fragments are solved based on the fitting equation and the known three-dimensional space coordinates of the fragments, including the following steps: Step S41, a differential equation of the motion of the fragments in the air is constructed: In formula (8), m is the mass of the fragment, C x is the frontal resistance coefficient of the fragment, p is the local air density, s is the windward area of the fragment perpendicular to the flight direction, T is the flight time of the fragment, and v is the instantaneous speed of the fragment. Step S42, the coordinates of the projectile before the static explosion and the coordinates at the explosion moment are obtained, the obtained trajectory of the fragments is combined, and an indirect solving method is adopted for solving; Firstly, the first time point of the trajectory is taken as a relative time point, and the position thereof is taken as a relative initial position, and thus: In the formula, T1 is the time corresponding to the first data point of the trajectory of the fragment, X' is the distance between the fragment and the explosion point at T1, v1 corresponds to the speed of the fragment at T1, v0 corresponds to the speed of the fragment at T0, x is the distance between the fragment and the explosion point at a certain time, and a is a speed attenuation coefficient, and thus: In formula (10), x(T) is the corresponding flight distance of the fragment at any T time, which is equivalent to: In formula (11): s tr is the distance of the fragment relative to the coordinates at T1 moment, t r is the time difference relative to T1 moment, using the obtained fragment trajectory data, the velocity v1 and the velocity attenuation coefficient α can be fitted. The footpoint coordinate of the fitted straight line of the above-mentioned flying of the fragments to the center of the warhead is calculated, and since the fitted straight line of the fragments does not necessarily pass through the center of the warhead, the footpoint coordinate is taken as the starting point of the fragment movement, and the distance between the first data point of the trajectory and the footpoint is X'; The velocity v1, the velocity attenuation coefficient a and the distance X' are substituted into formula (9) to obtain the initial velocity v0.

5. A system for testing the flying characteristics of fragments based on image processing technology, characterized in that, The fragment scattering characteristic test method based on the image processing technology is realized by the fragment scattering characteristic test method based on the image processing technology of any one of claims 1-4. 6.A computer device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein when the processor executes the computer program, the fragment scattering characteristic test method based on the image processing technology of any one of claims 1-4 is realized. 7.A computer readable storage medium, having a computer program stored thereon, wherein when the computer program is executed by a processor, the fragment scattering characteristic test method based on the image processing technology of any one of claims 1-4 is realized.

Images