Tunnel blasting lumpiness statistical method based on infrared imaging

Abstract

The invention discloses a tunnel blasting fragmentation statistical method based on infrared imaging, and relates to the technical field of tunnel engineering blasting, and the method comprises the following steps: S1, dual-band infrared imaging: employing long-wave infrared and medium-wave infrared dual-band cooperative imaging, penetrating smoke and enhancing fragment edge recognition; s2, fragment motion trajectory modeling: continuously capturing a fragment flight trajectory through high-speed infrared camera shooting, and calculating an initial velocity, an acceleration and parabola characteristics of the fragments in combination with an inertial motion equation; according to the invention, through a dual-band cooperative imaging technology, the LWIR band has excellent capability of penetrating through smoke and dust, and fragment contours can be obtained in a severe environment; the MWIR band is sensitive to a high-temperature target, fragment edge recognition is enhanced, the image signal-to-noise ratio and the target segmentation precision are remarkably improved through dual-band data fusion, and the problem that the imaging quality of traditional visible light camera shooting is poor in a complex blasting environment is effectively solved.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of tunnel engineering blasting technology, and particularly relates to a tunnel blasting block size statistical method based on infrared imaging. BACKGROUND

[0002] In the tunnel engineering blasting operation, the size and distribution of the blasting block size have an important influence on the subsequent construction process and engineering cost. Suitable blasting block size can improve the efficiency of slag removal and reduce the wear and tear of mechanical equipment, while too large or too small block size can lead to construction obstruction and cost increase. The traditional blasting block size analysis relies on visible light photography or manual sampling. The imaging quality of visible light photography is greatly limited in the night, underground tunnel or blasting dust scene. When the light is insufficient, the image clarity decreases, and it is difficult to accurately identify the fragments; a large amount of dust generated by blasting can block the line of sight, so that the visible light cannot effectively penetrate, resulting in that the complete fragment information cannot be obtained. Manual sampling is not only low in efficiency, but also has certain danger. In the complex environment after blasting, it is difficult for manual sampling to comprehensively and accurately collect representative samples, thereby affecting the accuracy of the block size analysis. The existing infrared technology is mostly used for static temperature monitoring, and has not been effectively combined with the fragment kinematics model and block size correlation analysis in the tunnel blasting scene. This makes it difficult to estimate the blasting block size in real time and accurately in the complex tunnel blasting environment, and it is difficult to meet the actual engineering needs. SUMMARY

[0003] The present application aims to overcome the deficiencies of the prior art and proposes a block size statistical method based on infrared imaging and fragment motion feature fusion. The method comprehensively utilizes key technologies such as dual-band infrared imaging, fragment motion trajectory modeling, thermodynamics-morphology correlation analysis and block size grading statistics, to realize accurate and real-time statistics of the tunnel blasting block size.

[0004] The present application is realized by the following technical solutions: A tunnel blasting block size statistical method based on infrared imaging, comprising the following steps: S1, dual-band infrared imaging: long-wave infrared and middle-wave infrared dual-band cooperative imaging is adopted to penetrate smoke and enhance fragment edge recognition; S2, fragment motion trajectory modeling: the flight trajectory of the fragment is continuously captured by a high-speed infrared camera, and the initial velocity, acceleration and parabolic characteristics of the fragment are calculated by combining the inertial motion equation; S3, thermodynamics-morphology correlation: the temperature image of the scene containing the blasting fragments is collected by utilizing the difference in the temperature field distribution of the fragment surface; S4, segmentation of the target region; a fragment mass-resistance model is established, and the equivalent volume is inversely deduced through the trajectory offset; S5, block size classification statistics: combined with the trajectory variance of the stability of the fragment movement and the thermal radiation area parameter, the blasting block size is divided into large / medium / small particle size grades, and the particle size distribution curve is output.

[0005] Further, the wavelength range of the long-wave infrared is 8-14 μm, and the wavelength range of the medium-wave infrared is 3-5 μm.

[0006] Further, the fragment motion trajectory modeling specifically includes the following steps: S21, identification point arrangement and coordinate system construction: a plurality of infrared reflective identification points are arranged around the shooting area, the identification points should be uniformly distributed, and the positional relationship between the identification points is clear, forming a stable reference coordinate system, the material of the identification point needs to have high infrared reflectivity, so that it can be clearly presented in the infrared image, and it is easy to identify and track; S22, high-speed infrared shooting and data acquisition: using a high-speed infrared camera (frame rate ≥ 200 fps) to continuously shoot the fragments after blasting; S23, establish motion equation: under the condition that the fragments are only subjected to gravity in the movement process, according to the kinematics principle, the inertial motion equation is established, in the two-dimensional plane, in the horizontal direction, the fragments do uniform linear motion, the equation is ; in the vertical direction, the fragments do uniform linear motion, the equation is ; In the formula, (x0, y0) is the initial position of the fragments, (v 0x ,v 0y ) is the initial velocity of the fragments in the horizontal and vertical directions, t is the motion time, and g is the acceleration of gravity; S24, trajectory calculation: from the constructed fragment motion trajectory, the position information (x i ,y i ) of the feature points at different times is obtained and the position sequence is generated, the position information is substituted into the inertial motion equation, and the parameters in the equation are fitted through numerical calculation methods such as least squares method, so as to calculate the initial velocity (v 0x ,v 0y ) of the fragments, the acceleration and the parabolic feature; S25, through experiments, the motion of fragments of different shapes at different speeds is simulated, the air resistance of the fragments is measured, and the experimental data is analyzed.

[0007] Further, the motion equation correction specifically includes the following steps: S231, considering the air resistance, the inertial motion equation is corrected, the direction of the air resistance is opposite to the direction of the fragment movement, and the size is proportional to the square of the speed, that is ; S232, in the horizontal direction, the motion equation becomes ; S233, in the vertical direction, the motion equation becomes ; wherein ，F d is the air resistance force generated by the air to the fragment motion, C d is the air resistance coefficient, v is the fragment motion speed, p is the air density, A is the projection area of the fragment in the motion direction, m is the mass of the fragment, is the acceleration of the fragment in the horizontal direction, indicating the rate of change of the horizontal direction speed with time, vx is the component of the fragment speed in the horizontal direction, vy is the component of the fragment speed in the vertical direction, is the acceleration of the fragment in the vertical direction, that is, the rate of change of the vertical direction speed with time.

[0008] Further, the thermodynamic-morphological correlation specifically comprises the following steps: S31, infrared thermal imaging data acquisition: at the moment of blasting and in a short time thereafter, an infrared thermal imaging device is used to quickly acquire a scene temperature image containing blasting fragments; S32, temperature threshold range segmentation: analyze the collected temperature image, and count the temperature distribution range of the background area in the image, while observing the temperature characteristics of the fragment area.

[0009] Further, the segmented target area specifically comprises the following steps: S41, binary processing: based on the set temperature threshold range, the temperature image is processed by using an image binary method, the pixel points with temperature higher than the threshold range are set as foreground and valued as 1; the pixel points with temperature lower than the threshold range are set as background and valued as 0; S42, trajectory offset measurement and calculation: the target area is segmented by using the difference in the surface temperature field distribution of the fragment; In a short time after blasting, the surface temperature of the fragment is higher than the surrounding environment, and in the temperature image obtained by infrared thermal imaging, the fragment area is segmented according to the temperature threshold range; By analyzing the difference between the actual motion trajectory of the fragment in the image and the ideal trajectory calculated theoretically, the trajectory offset is obtained; S43, equivalent volume backstepping model: the trajectory offset is substituted into the inertial motion equation, and the equivalent volume of the fragment is backstepped by iteration calculation or numerical solution method; The inertial motion equation involves multiple parameters, mass m, air resistance coefficient Cd , projection area A, etc., to establish the relationship between A and V.

[0010] Further, the trajectory variance is obtained by the following way: Based on the inertial motion equation, the air resistance is corrected to fit the theoretical trajectory, and the ideal position sequence is generated ; Calculate the deviation between the actual position and the theoretical position: , ; Horizontal direction variance calculation: ; Vertical direction variance calculation: ; Trajectory variance ; wherein, is the weight coefficient, ; Set the trajectory variance threshold range, if the trajectory variance is less than the threshold range, it is determined as large block degree, if the trajectory variance is greater than the threshold range, it is determined as small block degree, and if the trajectory variance is within the threshold range, it is determined as medium block degree.

[0011] Further, the thermal radiation area is obtained by the following steps: According to the temperature distribution of the explosion scene, manually set the threshold range to separate the high-temperature debris from the low-temperature background, and use the Otsu algorithm or the maximum inter-class variance method to automatically calculate the optimal threshold range:

[0012] In the formula, T represents the temperature, is the proportion of pixels below the threshold range, is the proportion of pixels above the threshold range, respectively represent the average temperature of the pixels below the threshold range and the average temperature of the pixels above the threshold range; Identify the independent debris area by 8-neighbor connectivity analysis, and assign a unique label to each connected domain. Count the number of pixels N pixel in each connected domain, and calculate the actual thermal radiation area according to the spatial resolution R of the camera: .

[0013] Set a thermal radiation actual area threshold range, when the thermal radiation area is greater than its threshold range and the trajectory variance is less than its threshold range, it is large particle size debris, when the thermal radiation area and the trajectory variance are within its threshold range, it is medium particle size debris; when the thermal radiation area is less than its threshold range and the trajectory variance is greater than its threshold range, it is small particle size debris.

[0014] Through the trajectory variance and the thermal radiation area grading result, the number of fragments of different particle size grades is counted, and finally a particle size distribution curve is output, so that the distribution of the blockiness after tunnel blasting is intuitively displayed.

[0015] Compared with the prior art, the present application has the following beneficial effects: 1. Dual-band cooperative imaging technology: the LWIR band has excellent penetration ability for smoke and dust, and can obtain fragment contours in harsh environments; the MWIR band is sensitive to high-temperature targets, enhances fragment edge recognition, and the fusion of dual-band data significantly improves the image signal-to-noise ratio and target segmentation accuracy, effectively solving the problem of poor imaging quality of traditional visible light cameras in complex blasting environments.

[0016] 2. Fine motion modeling and air resistance correction: a stable coordinate system is constructed by arranging infrared reflective markers, and the position information of the fragments is extracted from the high-speed image sequence; an air resistance model is introduced to correct the inertial motion equation, and the motion parameters of the fragments are accurately calculated by numerical integration iteration optimization.

[0017] 3. Temperature-morphology dynamic correlation analysis: based on the high-temperature characteristics of the fragments, the temperature threshold is adaptively set, combined with binary segmentation and trajectory offset, a mass-resistance-volume correlation model is established to realize non-contact size measurement, avoiding the safety risks and errors of manual sampling.

[0018] 4. Multi-dimensional grading statistics and visualization: fusion of motion stability and thermal radiation area dual parameters, dynamic division of particle size grades, output of cumulative distribution curve, intuitive feedback of blasting effect, providing scientific basis for adjusting charge quantity, hole spacing and other parameters. BRIEF DESCRIPTION OF DRAWINGS

[0019] The drawings described herein are used to provide further understanding of the present application, and form a part of the present application. The illustrative embodiments of the present application and their descriptions are used to explain the present application, and do not constitute an improper limitation on the present application. In the drawings: Figure 1 It is a schematic diagram of dual-band infrared camera and marker arrangement; Figure 2 It is a schematic diagram of fragment motion trajectory modeling; Figure 3 It is a temperature threshold segmentation and equivalent volume backstepping flowchart; Figure 4 It is a blasting blockiness distribution curve. DETAILED DESCRIPTION

[0020] The embodiments of the present application will be described in detail below with reference to the drawings.

[0021] The following specific examples illustrate the implementation of the present invention. Those skilled in the art can easily understand other advantages and effects of the present invention from the content disclosed in this specification. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of them. The present invention can also be implemented or applied through other different specific embodiments, and the details in this specification can also be modified or changed based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that, in the absence of conflict, the following embodiments and features in the embodiments can be combined with each other. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0022] Figure 4 The number of fragments in each grade is counted according to grading standards (such as large / medium / small particle size), and denoted as N. large N medium N small Cumulative percentage d i The horizontal axis represents the particle size threshold. The horizontal axis represents the particle size level (e.g., small → medium → large). The vertical axis represents the cumulative percentage (0%~100%). Use tools such as Matplotlib or Origin to generate step-like or smooth curves.

[0023] See Figures 1-4 One method includes the following steps: S1. Dual-band infrared imaging: Employs dual-band collaborative imaging of long-wave infrared and mid-wave infrared to penetrate smoke and enhance debris edge recognition. S2. Debris trajectory modeling: By continuously capturing the flight trajectory of debris using high-speed infrared cameras and combining it with the inertial motion equations, the initial velocity, acceleration, and parabolic characteristics of the debris are calculated. S3, Thermodynamic-morphological correlation: Utilizing the differences in temperature field distribution on the surface of fragments, acquire scene temperature images containing explosive fragments; S4. Segment the target region; establish a fragment mass-resistance model and infer the equivalent volume from the trajectory offset. S5. Particle Size Classification Statistics: Combining the trajectory variance of fragment motion stability and thermal radiation area parameters, the blasted particle size is classified into large / medium / small particle size levels, and the particle size distribution curve is output.

[0024] The wavelength range of long-wave infrared is 8-14 μm, and the wavelength range of mid-wave infrared is 3-5 μm, such as the FLIR X8580 series high-speed infrared camera, which supports synchronous dual-band acquisition. The camera frame rate is ≥200 fps, and the resolution is not less than 640*512 pixels. The long-wave infrared has strong ability to penetrate smoke and dust, and can obtain the basic object contour information in the harsh environment generated by blasting. The mid-wave infrared is more sensitive to high-temperature objects, and in the instant of blasting, the fragments generate high temperature due to explosion, and the mid-wave infrared can highlight these high-temperature fragments and enhance the edge contrast thereof; The fragment motion trajectory modeling specifically includes the following steps: S21, marker arrangement and coordinate system construction: a plurality of infrared reflective markers are arranged around the shooting area. These markers should be uniformly distributed and have a clear positional relationship with each other, forming a stable reference coordinate system. The material of the marker should have high infrared reflectivity so as to be clearly presented in the infrared image and easily identified and tracked; S22, high-speed infrared shooting and data acquisition: a high-speed infrared camera (frame rate ≥200 fps) is used to continuously shoot the fragments after blasting; S23, establishment of motion equation: it is assumed that the fragments are only subjected to gravity action (other factors such as air resistance are initially ignored) during the motion process, and an inertial motion equation is established according to the kinematics principle. In a two-dimensional plane, in the horizontal direction, the fragments do uniform linear motion, and the equation is ; in the vertical direction, the fragments do uniform accelerated linear motion, and the equation is , wherein (x0, y0) is the initial position of the fragments, (v 0x ,v 0y ) is the initial velocity of the fragments in the horizontal and vertical directions, t is the motion time, and g is the gravity acceleration; S24, trajectory calculation: the position information (x i ,y i ) of the feature points at different times is obtained from the constructed fragment motion trajectory, and a position sequence is generated. The position information is substituted into the inertial motion equation. The parameters in the equation are fitted by a numerical calculation method such as least square method, so as to calculate the initial velocity (v 0x ,v 0y ) of the fragments, the acceleration (the gravity acceleration g in the vertical direction) and the parabolic feature (the shape and trajectory of the parabola are determined by the initial velocity and the acceleration); S25, through multiple experiments, the motion of fragments of different shapes at different speeds is simulated, and the air resistance of the fragments is measured. The experimental data are analyzed, and a regression analysis method is used to determine the relationship between the air resistance coefficient and the shape and speed of the fragments.

[0025] The motion equation correction specifically includes the following steps: S231, considering the air resistance, the inertial motion equation is corrected, the direction of the air resistance is opposite to the direction of the fragment motion, and the size is proportional to the square of the speed, that is, S232, in the horizontal direction, the motion equation becomes S233, in the vertical direction, the motion equation becomes wherein ，F d is the air resistance force generated by the air to the fragment motion, C d is the air resistance coefficient, v is the fragment motion speed, p is the air density, A is the projection area of the fragment in the motion direction, m is the mass of the fragment, is the acceleration of the fragment in the horizontal direction, which represents the rate of change of the horizontal direction speed with time, vx is the component of the fragment speed in the horizontal direction, vy is the component of the fragment speed in the vertical direction, is the acceleration of the fragment in the vertical direction, that is, the rate of change of the vertical direction speed with time, the corrected equation is combined with the previously obtained characteristic point position information, and the numerical integration method is used to re-calculate the initial speed, acceleration and parabolic characteristic of the fragment, so as to improve the accuracy of the motion parameter calculation; The thermodynamic-morphological correlation specifically includes the following steps: S31, infrared thermal imaging data acquisition: at the moment of blasting and in a short time thereafter, an infrared thermal imaging device is used to quickly acquire a scene temperature image containing the blasting fragments. Since blasting will make the fragments obtain high temperature, the surface temperature is significantly higher than the surrounding environment, and there will be obvious temperature difference in the infrared thermal imaging image; S32, temperature threshold range segmentation: analyze the collected temperature image, count the temperature distribution range of the background area (non-fragment area) in the image, and observe the temperature characteristics of the fragment area at the same time. Through multiple experiments or based on experience, a suitable temperature threshold range is set. For example, if the background temperature range is 20-30℃, and the blasting fragment temperature is 100-500℃, the temperature threshold range can be initially set to 50℃. This threshold range should be able to effectively distinguish the fragment area and the background area, and avoid missegmentation; The segmented target area specifically includes the following steps: ​​​S41, binary processing: based on the set temperature threshold range, the temperature image is processed by image binarization method. The pixel points with temperature higher than the threshold range are set as foreground (i.e. the fragment area), and the value is 1; the pixel points with temperature lower than the threshold range are set as background, and the value is 0. In this way, the target area where the fragments are located can be segmented from the original temperature image, and a binary image containing only fragments is obtained; S42, trajectory offset measurement and calculation: the target area is segmented by the difference of the surface temperature field distribution of the fragments (high temperature residue at the moment of blasting). In a short time after blasting, the surface temperature of the fragments is higher than the surrounding environment, and the fragment area is segmented according to the temperature threshold range in the temperature image obtained by infrared thermal imaging. By analyzing the difference between the actual motion trajectory of the fragments in the image and the ideal trajectory calculated theoretically, the trajectory offset is obtained. For example, in the vertical direction, the vertical distance between the actual landing point of the fragments and the landing point of the ideal parabolic trajectory is measured as the trajectory offset in the vertical direction; in the horizontal direction, the difference between the actual horizontal displacement and the ideal horizontal displacement is measured as the trajectory offset in the horizontal direction; S43, equivalent volume backstepping model: the trajectory offset is substituted into the inertial motion equation, and the equivalent volume of the fragments is backstepped by iterative calculation or numerical solution method. Since the inertial motion equation involves multiple parameters, such as mass m, air resistance coefficient C d , projection area A, etc., and the projection area A is related to the equivalent volume V (for simple shapes such as spherical , , the relationship between A and V can be established; for irregular shapes, the relationship can be established by shape factor, etc.). Under the condition that other parameters (such as air density p, acceleration of gravity g, etc.) are known, the value of equivalent volume V is adjusted so that the trajectory calculated by the model is consistent with the measured trajectory offset, thereby determining the equivalent volume of the fragments; The contour of the fragments is obtained by infrared image segmentation, and the following feature parameters are extracted to establish the shape factor. The major axis a and the minor axis b are obtained by ellipse fitting. The aspect ratio (a / b) is calculated. The shape of the fragments is preliminarily judged as circular: ; flat or elongated shape: . The basic projection area is calculated. Based on the perimeter C and the area S of the fragment contour, the is calculated, which quantifies the irregularity of the edge and corrects the basic projection area. The corrected projection area is . Assuming that the fragments are ellipsoids, the volume formula is , where the third axis c can be determined as follows: if the fragments are close to circular (a / b ), take (combine irregularity correction). If the fragments are significantly elongated or flat (a / b ), take (Prefer the impact of aspect ratio).

[0026] Trajectory variance is obtained by: Considering air resistance correction fitting theoretical trajectory based on inertial motion equation, generating ideal position sequence ; Calculate the deviation between actual position and theoretical position: , ; Horizontal direction variance calculation: ; Vertical direction variance calculation: ; Trajectory variance ; Where, is the weight coefficient, ; For motion stability, calculate the variance of the trajectory of each fragment during flight, and the fragment with smaller variance moves relatively stably, usually indicating that its size is larger; The motion of the fragment with larger variance is less stable, and the size may be smaller.

[0027] Thermal radiation area is obtained by the following steps: The measurement of thermal radiation area needs to be realized by the process of threshold range segmentation, binarization, connected component analysis and pixel statistics, and is not the direct output result of thermal imaging equipment, but needs to be realized by combining data processing algorithm.

[0028] The original data collected by the infrared camera is a temperature matrix, and each pixel corresponds to a temperature value (unit: ℃ or K). Usually it is 16-bit RAW data or standard image format (such as JPEG, PNG), which contains temperature information and spatial coordinates. The spatial resolution of the camera needs to be calibrated in advance (such as 1 pixel = 0.1mm) to ensure accurate conversion of physical size.

[0029] According to the typical temperature distribution of the explosion scene, such as the background temperature (20-30℃) and the fragment temperature (100-500℃), the threshold range (such as 50℃) is manually set to separate the high-temperature fragments from the low-temperature background. Otsu algorithm or maximum inter-class variance method is used to automatically calculate the optimal threshold range:

[0030] In the formula, T represents the temperature, is the proportion of pixels below the threshold range, is the proportion of pixels above the threshold range, respectively represent the average temperature of pixels below the threshold range and the average temperature of pixels above the threshold range; Convert the original heat map to a binary image, where the debris area is marked as foreground (assigned value 1) and the background area is marked as background (assigned value 0). Then, eliminate small noise by morphological opening operation, and fill the internal holes of the debris by closing operation to ensure the integrity of the target area. Identify independent debris areas by 8-neighbor connectivity analysis, and assign a unique label to each connected domain. Count the number of pixels N pixel of each connected domain. Calculate the actual area A according to the spatial resolution R of the camera (unit: mm / pixel): .

[0031] Combine trajectory variance and thermal radiation area to improve the robustness and accuracy of block size classification: trajectory variance reflects the dynamic characteristics of debris (such as mass, air resistance), which is indirectly related to size. Large debris has high mass and inertia, so the motion trajectory is more stable (small variance). Thermal radiation area directly reflects the static geometric size of debris (projected area), but is affected by shape complexity. The combined criterion avoids misjudgment of a single parameter (such as small debris with regular shape may have small variance). At the same time, dynamic and static information are combined to reduce environmental noise interference (such as local occlusion of thermal imaging by dust).

[0032] Calculate the variance of the trajectory of each debris during flight to assess motion stability. Measure the thermal radiation area of the debris through infrared thermal imaging images. Combine the equivalent volume information of the debris to set the classification criteria. For example, if the thermal radiation area is greater than S 1 and the trajectory variance is less than σ 2 1, it is a large-diameter debris; if the thermal radiation area is between S 2- S 1 and the trajectory variance is between σ 2 1- σ 2 2, it is a medium-diameter debris; if the thermal radiation area is less than S 2 and the trajectory variance is greater than σ 2 2, it is a small-diameter debris ( S 1, S 2, σ 2 1, σ 2 2 can be determined according to experiments and field experience). According to the classification criteria, count the number of debris in different particle size grades. Take the particle size grade as the horizontal coordinate and the percentage of particles smaller than a certain particle size as the vertical coordinate to draw the particle size distribution curve, which intuitively shows the distribution of block size after tunnel blasting. According to the classification results of trajectory variance and thermal radiation area, count the number of debris in different particle size grades, and finally output the particle size distribution curve to intuitively show the distribution of block size after tunnel blasting.

[0033] The above described embodiments are only to illustrate the preferred embodiments of the present application, and are not intended to limit the scope of the present application. Any modification and improvement of the technical solutions of the present application made by those skilled in the art without departing from the design spirit of the present application shall fall within the protection scope of the present application.

Claims

1. A tunnel blasting fragmentation statistical method based on infrared imaging, characterized in that, The method comprises the following steps: S1, dual-band infrared imaging: long-wave infrared and mid-wave infrared dual-band cooperative imaging is adopted to penetrate smoke and enhance the edge recognition of fragments; S2, fragment trajectory modeling: the flight trajectory of the fragments is continuously captured by a high-speed infrared camera, and the initial speed, acceleration and parabolic characteristics of the fragments are calculated by combining the inertial motion equation; S3, thermodynamic-morphological correlation: the temperature image of the scene containing the blasting fragments is collected by utilizing the difference in the surface temperature field distribution of the fragments; S4, target region segmentation; A fragment mass-resistance model is established, and the equivalent volume is inversely deduced through the trajectory deviation; S5, block size classification statistics: the blasting block size is divided into large, medium and small particle size grades by combining the trajectory variance of the fragment motion stability and the thermal radiation area parameter, and a particle size distribution curve is output.

2. The method for tunnel blasting fragmentation statistical analysis based on infrared imaging according to claim 1, characterized in that, The wavelength range of the long-wave infrared is 8-14 μm, and the wavelength range of the mid-wave infrared is 3-5 μm.

3. The infrared imaging based tunnel blast fragmentation size statistics method according to claim 1, characterized in that, The fragment trajectory modeling specifically comprises the following steps: S21, identification point arrangement and coordinate system construction: a plurality of infrared reflective identification points are arranged around the shooting area, the identification points are uniformly distributed and have a clear positional relationship with each other, a stable reference coordinate system is formed, and the material of the identification points has high infrared reflectivity; S22, high-speed infrared shooting and data acquisition: a high-speed infrared camera (frame rate ≥ 200 fps) is used to continuously shoot the fragments after blasting; S23, establish motion equation: in the case of only gravity acting on the debris during the movement, according to the kinematic principle, the inertial motion equation is established, in the two-dimensional plane, in the horizontal direction, the debris does uniform linear motion, the equation is ; In the vertical direction, the fragments do uniform linear motion, the equation is ; where (x0, y0) is the initial position of the fragment, (v 0x ,v 0y ) is the initial velocity of the fragment in the horizontal and vertical direction components, t is the time of movement, and g is the acceleration of gravity; S24, trajectory estimation: from the constructed trajectory of the debris, the position information (x i ,y i ) of the feature points at different time is obtained and a position sequence is generated, the position information is substituted into the inertial motion equation, the parameters in the equation are fitted through the least square method, so as to estimate the initial velocity (v 0x ,v 0y ), acceleration and parabolic characteristics of the debris; S25, through experiments, the motion of fragments of different shapes at different speeds is simulated, the air resistance of the fragments is measured, and experimental data is analyzed.

4. The method for tunnel blasting fragmentation statistical analysis based on infrared imaging according to claim 3, characterized in that, The motion equation correction specifically comprises the following steps: S231、Consider the air resistance, the inertial motion equation is modified, the direction of air resistance is opposite to the direction of debris movement, its size is proportional to the square of speed, namely ; S232、In the horizontal direction, the equation of motion becomes ; S233、In the vertical direction, the equation of motion becomes ; wherein ，F d Fdragis the drag force generated by air on the debris motion, C d Cdragis the air drag coefficient, v V is the debris motion velocity, ρ P is the air density, A A is the debris projected area in the motion direction, m is the debris mass, a is the debris acceleration in the horizontal direction, representing the rate of change of the horizontal velocity with time, vx Vx is the debris velocity component in the horizontal direction, vy Vy is the debris velocity component in the vertical direction, a is the debris acceleration in the vertical direction, representing the rate of change of the vertical velocity with time.

5. The infrared imaging based tunnel blast fragmentation size statistics method of claim 1, wherein, The thermodynamic-morphological correlation specifically comprises the following steps; S31, infrared thermal imaging data acquisition: the temperature image of the scene containing the blasting fragments is quickly acquired by using an infrared thermal imaging device at the moment of blasting and in a short time thereafter; S32, temperature threshold range segmentation: the acquired temperature image is analyzed, the temperature distribution range of the background region in the image is counted, and the temperature characteristics of the fragment region are observed.

6. The method for tunnel blasting fragmentation statistical analysis based on infrared imaging according to claim 5, characterized in that, The target region segmentation specifically comprises the following steps: S41, binary processing: based on the set temperature threshold range, the temperature image is processed by using an image binary method, the pixel points with a temperature higher than the threshold range are set as the foreground and assigned a value of 1, and the pixel points with a temperature lower than the threshold range are set as the background and assigned a value of 0; S42, trajectory deviation measurement and calculation: the target region is segmented by utilizing the difference in the surface temperature field distribution of the fragments; In a short time after blasting, the surface temperature of the fragments is higher than that of the surrounding environment, the fragment region is segmented from the temperature image acquired by the infrared thermal imaging according to the temperature threshold range; The trajectory deviation is obtained by analyzing the difference between the actual motion trajectory of the fragments in the image and the ideal trajectory calculated theoretically; S43, equivalent volume inverse deduction model: the trajectory deviation is substituted into the inertial motion equation, and the equivalent volume of the fragments is inversely deduced by iterative calculation or numerical solution method. The trajectory variance is obtained in the following manner:

7. The infrared imaging based tunnel blast fragmentation size statistics method of claim 3, wherein, The trajectory variance threshold range is set, the trajectory variance calculation value is compared with the trajectory variance threshold range, and the size of the blasting block is determined. Based on the inertial motion equation, the ideal position sequence is generated by considering the air resistance correction fitting theoretical trajectory ; calculating a deviation of the actual position from the theoretical position: , ; Horizontal direction variance calculation: ; Vertical direction variance calculation: ; Trajectory variance ; wherein is a weight coefficient, ; ​ 8. The infrared imaging based tunnel blast fragmentation size statistics method of claim 1, wherein, The thermal radiation area is obtained by the following steps: According to the temperature distribution of the blasting scene, the threshold range is manually set to separate the high-temperature debris from the low-temperature background, and the optimal threshold range is automatically calculated by using the Otsu algorithm or the maximum inter-class variance method: wherein T represents temperature, is a proportion of pixels below the threshold range, is a proportion of pixels above the threshold range, respectively represent the average temperature of pixels below the threshold range and the average temperature of pixels above the threshold range; The independent fragment regions are identified by 8-neighbor connectivity analysis, and a unique label is assigned to each connected domain. The number of pixels Npixel in each connected domain is counted, and the actual area of thermal radiation is calculated according to the spatial resolution R of the camera: ; A threshold range of the actual thermal radiation area is set, the thermal radiation area is compared with the threshold range of the actual thermal radiation area, and the blasting size is determined by further combining the comparison result of the trajectory variance.

9. The infrared imaging based tunnel blast fragmentation size statistics method of claim 8, wherein, Through the trajectory variance and the classification result of the thermal radiation area, the number of debris of different particle size grades is counted, and finally a particle size distribution curve is output to intuitively show the distribution of the blasting size after the tunnel blasting.

Images