Traumatic bleeding detection method and system based on infrared thermal imaging and ultrasonic multiple modes

Abstract

The invention provides a trauma bleeding detection method and system based on infrared thermal imaging and ultrasonic multi-modality, and relates to the technical field of medical image processing, and the method comprises the steps: collecting and registering infrared thermal imaging and ultrasonic images of a trauma part, executing adaptive shear wave transformation to extract a feature coefficient, constructing a thermoacoustic coupling feature matrix, and obtaining a thermoacoustic coupling feature matrix; and reconstructing a vascular structure to generate a blood flow velocity field distribution diagram, establishing a time sequence mapping relation between local blood flow volume and accumulated bleeding volume, and calculating bleeding dynamic parameters by adopting a variational particle filtering method. The method can accurately evaluate the duration of traumatic hemorrhage and the tissue injury degree, and provides an accurate diagnosis basis for clinical first aid.

Description

Technical Field

[0001] The present invention relates to the technical field of medical image processing, and particularly to a method and system for detecting traumatic bleeding based on infrared thermography and ultrasound multimodality. Background Art

[0002] In the field of trauma first aid, bleeding is one of the main causes of patient death. Rapid and accurate detection of traumatic bleeding is of great significance for clinical decision-making and treatment. At present, the detection of traumatic bleeding mainly relies on methods such as clinical symptom assessment, laboratory tests, and imaging examinations. Traditional detection methods include imaging techniques such as ultrasound examination, CT scan, MRI, and laboratory index monitoring such as hemoglobin and hematocrit. With the development of medical imaging technology, non-invasive detection techniques such as infrared thermography and ultrasound imaging have gradually been applied to clinical trauma assessment. Infrared thermography can reflect the state of deep tissues by detecting changes in surface temperature distribution, while ultrasound imaging can provide tissue structure and blood flow information. The combination of these two technologies provides a new idea for the detection of traumatic bleeding.

[0003] However, the existing traumatic bleeding detection technologies still have several defects and deficiencies. The detection methods of single imaging modality have limited information and are difficult to comprehensively reflect the complex state of bleeding. Infrared thermography can only reflect the surface temperature distribution and cannot accurately evaluate the bleeding condition of deep tissues; while simple ultrasound examination can show the tissue structure, but the detection sensitivity for small blood vessels and low-speed blood flow is insufficient, and it is easy to miss the diagnosis; the existing multimodal image fusion methods lack an effective spatio-temporal registration mechanism and are difficult to achieve the precise correspondence between infrared thermography and ultrasound images, resulting in the distortion of the fused information and the deviation of clinical judgment. Especially in dynamic scenarios, the differences in the acquisition time and spatial resolution of the two modalities make the traditional registration methods ineffective; the existing blood volume assessment methods are mostly based on static models and lack the ability to track the dynamic process of bleeding in real time. Traditional methods cannot effectively reflect the changes in bleeding rate and cumulative effect, nor can they predict the bleeding duration and tissue damage degree, thus affecting the timeliness and accuracy of clinical treatment decisions. Summary of the Invention

[0004] An embodiment of the present invention provides a method and system for detecting traumatic bleeding based on infrared thermography and ultrasound multimodality, which can solve the problems in the prior art.

[0005] In the first aspect of the embodiment of the present invention, a method for detecting traumatic bleeding based on infrared thermography and ultrasound multimodality is provided, including: Collecting infrared thermography images and ultrasound images of the trauma site, and performing spatio-temporal registration on the images through a calibration board to obtain registered image data; Perform an adaptive shearlet transform on the registered image data to obtain temperature field coefficients and acoustic coefficients. Substitute the temperature field coefficients into the heterogeneous heat diffusion equation to calculate the deep temperature distribution map and obtain the anisotropic heat diffusion coefficient. Perform pulse compression and phase correction on the acoustic coefficients to obtain the tissue stratification map; According to the deep temperature distribution map and the tissue stratification map, construct a thermoacoustic coupling feature matrix using a complex-valued neural network. Extract the blood flow velocity vector from the Doppler shift signal based on the thermoacoustic coupling feature matrix. Reconstruct the vascular structure using Marching Cubes interpolation to generate the blood flow velocity field distribution map; Calculate the local blood flow flux according to the blood flow velocity field distribution map, and perform a volume integral in combination with the anisotropic heat diffusion coefficient to obtain the cumulative blood loss. Establish a temporal mapping relationship between the local blood flow flux and the cumulative blood loss; Use the variational particle filter method to dynamically track the local blood flow flux in the temporal mapping relationship, calculate the corresponding time-varying coefficient, and construct a state equation to obtain the bleeding dynamic parameters; Calculate the bleeding duration and the degree of tissue damage according to the bleeding dynamic parameters, and determine the risk level and treatment plan.

[0006] In an alternative embodiment, performing an adaptive shearlet transform on the registered image data to obtain temperature field coefficients and acoustic coefficients, substituting the temperature field coefficients into the heterogeneous heat diffusion equation to calculate the deep temperature distribution map and obtain the anisotropic heat diffusion coefficient, and performing pulse compression and phase correction on the acoustic coefficients to obtain the tissue stratification map includes: Perform an adaptive shearlet transform on the registered image data, where the adaptive shearlet transform performs a joint decomposition on the registered image data by setting a direction-selective filter bank, constructs a shearlet basis function according to the shear parameter and the translation parameter, and obtains a coefficient matrix through a two-dimensional integral operation on the shearlet basis function. Calculate the direction angle parameter based on the arctangent function of the shear parameter, and optimize the direction angle parameter to achieve multi-directional adaptive decomposition to obtain the temperature field coefficients and the acoustic coefficients; Substitute the temperature field coefficients into the heterogeneous heat diffusion equation, where the heterogeneous heat diffusion equation includes a tissue density term, a specific heat capacity term, a spatially varying thermal conductivity tensor term, and a heat source term. Calculate the coefficient term by multiplying the tissue density term and the specific heat capacity term, calculate the inner product of the thermal conductivity tensor term and the temperature gradient to obtain the diffusion term, substitute the coefficient term, the diffusion term, and the heat source term into the heterogeneous heat diffusion equation, and use the finite difference method to solve the heterogeneous heat diffusion equation to obtain the deep temperature distribution map, and calculate the anisotropic heat diffusion coefficient based on the ratio of the maximum thermal conductivity to the minimum thermal conductivity; Pulse compression processing is performed on the acoustic coefficient to obtain a compressed echo signal, phase correction is performed by calculating the difference between the original phase and the reference phase, the acoustic reflection coefficient is calculated according to the ratio of the difference and the sum of the acoustic impedances of adjacent tissues, and a tissue layer map is reconstructed according to the acoustic reflection coefficient, the compressed echo signal and the phase correction result.

[0007] In an optional embodiment, solving the heterogeneous heat diffusion equation using a finite difference method to obtain a deep temperature distribution map includes: Discretize the time and space coordinates, divide the time dimension into multiple time layers according to the time step, divide the space dimension into a three-dimensional grid structure according to the space step, construct a differential grid, and obtain the initial temperature value at the grid point of the differential grid; performing differential discretization processing on the heterogeneous heat diffusion equation, discretizing the time derivative term using a first-order forward difference format, discretizing the space derivative term using a second-order central difference format, interpolating the thermal conductivity tensor at half-grid points to obtain a discrete thermal conductivity value, and substituting the discrete time derivative term, the discrete space derivative term, and the discrete thermal conductivity value into the heterogeneous heat diffusion equation to establish a differential equation; Constructing boundary conditions at the boundary positions of the differential grid, discretizing the boundary conditions of the temperature boundary, the heat flow boundary, and the convection boundary, and combining the boundary conditions with the differential equations to construct a complete differential equation group; The complete differential equations are decomposed and solved in the x-, y-, and z-directions in sequence using an alternating direction implicit method, a tridiagonal equation system is constructed for each direction, and an iterative solution is performed using a chasing method to obtain the temperature distribution of the new time layer; When the temperature difference between adjacent time layers is less than a preset convergence threshold, a deep temperature distribution map is generated based on the final converged temperature distribution.

[0008] In an optional embodiment, a complex-valued neural network is used to construct a thermoacoustic coupling characteristic matrix based on the deep temperature distribution map and the tissue layer map, a blood flow velocity vector is extracted from the Doppler frequency shift signal based on the thermoacoustic coupling characteristic matrix, and the vascular structure is reconstructed using Marchi cube interpolation to generate a blood flow velocity field distribution map, including: Normalizing the deep temperature distribution map and applying Gaussian filtering to obtain a temperature field gradient, performing adaptive threshold segmentation on the tissue layer map to obtain a tissue interface normal vector, and combining the temperature field gradient and the tissue interface normal vector to form an initial feature vector; Inputting the initial eigenvector into a complex-valued neural network, wherein the complex-valued neural network adopts a three-layer network structure and has a piecewise linear unit function as an activation function, and processing the real and imaginary parts of the initial eigenvector by the complex-valued neural network to obtain a thermoacoustic coupling characteristic matrix; The Doppler frequency shift signal is subjected to Hilbert transform to obtain a complex analytic signal, and the complex analytic signal is subjected to tensor product operation with the thermoacoustic coupling characteristic matrix to obtain a three-dimensional blood flow velocity vector field; Calculate the velocity divergence distribution according to the three-dimensional blood flow velocity vector field, determine the vascular boundary point set, and input the vascular boundary point set into the cubic Marching Cubes interpolation algorithm to obtain the vascular surface; Within the range of the vascular surface, use radial basis functions to perform interpolation calculation on the three-dimensional blood flow velocity vector field to obtain a continuous blood flow velocity field distribution map.

[0009] In an alternative embodiment, calculating the velocity divergence distribution according to the three-dimensional blood flow velocity vector field, determining the vascular boundary point set, and inputting the vascular boundary point set into the cubic Marching Cubes interpolation algorithm to obtain the vascular surface includes: Perform scale pyramid decomposition on the three-dimensional blood flow velocity vector field to obtain velocity vector fields at multiple scale levels, calculate the velocity divergence distribution for each velocity vector field respectively to determine the multi-scale divergence distribution, and based on the multi-scale divergence distribution, perform feature enhancement through average pooling and a multi-layer perceptron to obtain an enhanced velocity divergence distribution; Construct a graph structure including a vertex set, an edge set, and an adjacency matrix from the enhanced velocity divergence distribution, input the graph structure into a graph attention network, perform feature iterative update to obtain vascular boundary features; Extract the skeleton vascular boundary point set based on the vascular boundary features and perform recursive refinement to obtain a complete vascular boundary point set; Construct a joint loss function including a topological consistency constraint term and a geometric constraint term, optimize the complete vascular boundary point set based on the joint loss function to determine an optimized vascular boundary point set; Divide the optimized vascular boundary point set into local surface patches and establish a parameterized grid, calculate the control vertices of the parameterized grid, construct cubic B-spline basis functions based on the control vertices, perform tensor product operation with the position information of the control vertices to obtain local surface reconstruction results, and after applying curvature continuity constraints and boundary continuity constraints to the local surface reconstruction results, splice them to obtain the vascular surface.

[0010] In an alternative embodiment, a variational particle filtering method is used to dynamically track the local blood flux in the temporal mapping relationship, calculate the corresponding time-varying coefficients, and construct a state equation to obtain bleeding dynamic parameters, including: Initialize multiple particles according to the temporal mapping relationship and assign initial weights, where each particle corresponds to a state of the local blood flux; An observation likelihood function is established to describe the relationship between the cumulative blood loss and the local blood flow rate, and a state transition probability function is established to describe the variation law of the local blood flow rate; the observation likelihood function and the state transition probability function are combined to construct a variational objective function; Based on the variational objective function, the weights of the particles are updated through iterative optimization; the optimal estimated value of the local blood flow rate is calculated according to the updated particle weights; importance resampling is performed to eliminate particle degradation, and the dynamic change result of the local blood flow rate is obtained; According to the dynamic change result of the local blood flow rate, the variation law of the blood flow rate in different time periods is analyzed; based on the variation law, the scale coefficient, delay coefficient and attenuation coefficient of the hemodynamic characteristics are calculated, and time series smoothing processing is performed to determine each smoothing coefficient; The state transition matrix of the nonlinear dynamic model is constructed by using the smoothing coefficient, and the bleeding dynamic characteristic state equation is constructed by combining the system noise term with the state transition matrix; the bleeding dynamic characteristic state equation is solved to calculate the bleeding rate, cumulative blood loss and tissue damage degree, and the bleeding dynamic parameters are determined.

[0011] In an alternative embodiment, updating the weights of the particles through iterative optimization based on the variational objective function includes: Calculating the partial derivative of the variational objective function with respect to each particle weight to obtain the sensitivity eigenvalue; extracting the sign of the sensitivity eigenvalue and multiplying it by the power term of the sensitivity eigenvalue to obtain the weight update direction; Calculating the KL divergence of the particle weight distributions in two adjacent iterations to construct an adaptive step size, multiplying the adaptive step size by the weight update direction to determine the weight update amount, and adding the weight update amount to the particle weight of the current iteration to obtain the updated particle weight; Calculating the maximum and minimum values of the updated particle weights, and performing renormalization processing on the updated particle weights to obtain the renormalized particle weights; Calculating the Euclidean distance between the particle weights in two adjacent iterations to obtain a stability index, and when the stability index is less than a preset stability threshold for a preset number of consecutive times, determining the particle weight of the current iteration as the optimization result.

[0012] In the second aspect of the embodiments of the present invention, a trauma bleeding detection method system based on infrared thermal imaging and ultrasound multimodality is provided, including: A first unit for collecting infrared thermal imaging images and ultrasound images of the trauma site, and performing spatio-temporal registration on the images through a calibration plate to obtain registered image data; A second unit, configured to perform an adaptive shearlet transform on the registered image data to obtain temperature field coefficients and acoustic coefficients, substitute the temperature field coefficients into a heterogeneous heat diffusion equation to calculate a deep temperature distribution map and obtain an anisotropic heat diffusion coefficient, and perform pulse compression and phase correction on the acoustic coefficients to obtain a tissue stratification map; A third unit, configured to construct a thermoacoustic coupling feature matrix using a complex-valued neural network according to the deep temperature distribution map and the tissue stratification map, extract a blood flow velocity vector from a Doppler shift signal based on the thermoacoustic coupling feature matrix, and reconstruct a blood vessel structure using Marching Cubes interpolation to generate a blood flow velocity field distribution map; A fourth unit, configured to calculate a local blood flow rate according to the blood flow velocity field distribution map, perform a volume integral in combination with the anisotropic heat diffusion coefficient to obtain an accumulated blood loss amount, and establish a time series mapping relationship between the local blood flow rate and the accumulated blood loss amount; A fifth unit, configured to dynamically track the local blood flow rate in the time series mapping relationship using a variational particle filter method, calculate a corresponding time-varying coefficient, and construct a state equation to obtain bleeding dynamic parameters; A sixth unit, configured to calculate a bleeding duration and a tissue damage degree according to the bleeding dynamic parameters, and determine a risk level and a treatment plan.

[0013] In a third aspect of the embodiments of the present invention, there is provided an electronic device, including: A processor; A memory for storing instructions executable by the processor; Wherein, the processor is configured to call the instructions stored in the memory to execute the method described above.

[0014] In a fourth aspect of the embodiments of the present invention, there is provided a computer-readable storage medium, on which computer program instructions are stored, and when the computer program instructions are executed by a processor, the method described above is implemented.

[0015] In the embodiments of the present invention, through the multi-modal fusion of infrared thermography and ultrasound technology, the accurate detection of traumatic bleeding is realized, the technical problem that it is difficult to accurately evaluate the bleeding condition of deep tissues by traditional methods is solved, and the diagnostic accuracy of the bleeding location and degree is improved; the registered image data is processed using adaptive shearlet transform and complex-valued neural network, a thermoacoustic coupling feature matrix is constructed, the tissue structure and blood flow dynamic characteristics can be deeply analyzed, the blood vessel structure is reconstructed using Marching Cubes interpolation, the accurate calculation of the bleeding volume is realized, and the sensitivity and specificity of traumatic bleeding detection are greatly improved; the variational particle filter method is introduced to dynamically track and analyze the blood flow rate, an association model between bleeding dynamic parameters and tissue damage degree is established, the bleeding development trend can be predicted in real time, providing a scientific basis for formulating individualized first aid strategies clinically, and significantly improving the timeliness and success rate of trauma treatment. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 This is a schematic flowchart of the trauma bleeding detection method based on infrared thermography and ultrasound multimodality according to an embodiment of the present invention; Figure 2 This is a comparison chart of the calculation efficiency for solving the heterogeneous heat diffusion equation; Figure 3 This is a comparison chart of the curvature continuity of the vascular surface reconstruction by different methods; Figure 4 This is a flowchart for optimizing the particle weights of multimodal medical images. Detailed implementation manners

[0017] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0018] The technical solutions of the present invention will be described in detail below with specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be repeated in some embodiments.

[0019] Figure 1 This is a schematic flowchart of the trauma bleeding detection method based on infrared thermography and ultrasound multimodality according to an embodiment of the present invention, as Figure 1 shown, the method includes: Collect infrared thermography images and ultrasound images of the trauma site, perform spatio-temporal registration on the images through a calibration plate to obtain registered image data; Perform adaptive shearlet transform on the registered image data to obtain temperature field coefficients and acoustic coefficients, substitute the temperature field coefficients into the heterogeneous heat diffusion equation to calculate the deep temperature distribution map and obtain the anisotropic heat diffusion coefficient, and perform pulse compression and phase correction on the acoustic coefficients to obtain the tissue stratification map; According to the deep temperature distribution map and the tissue stratification map, construct a thermoacoustic coupling feature matrix using a complex-valued neural network, extract the blood flow velocity vector from the Doppler shift signal based on the thermoacoustic coupling feature matrix, and reconstruct the vascular structure using Marching Cubes interpolation to generate a blood flow velocity field distribution map; Calculate the local blood flow flux according to the blood flow velocity field distribution map, perform volume integration in combination with the anisotropic heat diffusion coefficient to obtain the cumulative blood loss, and establish a temporal mapping relationship between the local blood flow flux and the cumulative blood loss; The variational particle filter method is used to dynamically track the local blood flow in the time series mapping relationship, calculate the corresponding time-varying coefficients, construct a state equation to obtain the bleeding dynamic parameters; Calculate the bleeding duration and tissue damage degree according to the bleeding dynamic parameters, and determine the risk level and treatment plan.

[0020] In an alternative embodiment, perform an adaptive shearlet transform on the registered image data to obtain a temperature field coefficient and an acoustic coefficient. Substitute the temperature field coefficient into the heterogeneous heat diffusion equation to calculate the deep temperature distribution map and obtain the anisotropic heat diffusion coefficient. Perform pulse compression and phase correction on the acoustic coefficient to obtain a tissue stratification map, including: Perform an adaptive shearlet transform on the registered image data, where the adaptive shearlet transform performs a joint decomposition on the registered image data by setting a direction-selective filter bank, constructs a shearlet basis function according to the shear parameter and the translation parameter, and obtains a coefficient matrix through a two-dimensional integral operation on the shearlet basis function. Calculate the direction angle parameter based on the arctangent function of the shear parameter, and optimize the direction angle parameter to achieve multi-directional adaptive decomposition, obtaining a temperature field coefficient and an acoustic coefficient; Substitute the temperature field coefficient into the heterogeneous heat diffusion equation, where the heterogeneous heat diffusion equation includes a tissue density term, a specific heat capacity term, a spatially varying thermal conductivity tensor term, and a heat source term. Calculate the product of the tissue density term and the specific heat capacity term to obtain a coefficient term, calculate the inner product of the thermal conductivity tensor term and the temperature gradient to obtain a diffusion term, substitute the coefficient term, the diffusion term, and the heat source term into the heterogeneous heat diffusion equation, and use the finite difference method to solve the heterogeneous heat diffusion equation to obtain the deep temperature distribution map, and calculate the anisotropic heat diffusion coefficient based on the ratio of the maximum thermal conductivity to the minimum thermal conductivity; Perform pulse compression processing on the acoustic coefficient to obtain a compressed echo signal, perform phase correction by calculating the difference between the original phase and the reference phase, calculate the acoustic reflection coefficient according to the ratio of the difference value to the sum value of the acoustic impedances of adjacent tissues, and reconstruct a tissue stratification map according to the acoustic reflection coefficient, the compressed echo signal, and the phase correction result.

[0021] In a specific embodiment, the process of performing an adaptive shearlet transform on the registered image data includes: setting a direction-selective filter bank, which includes band-pass filters in 8 directions, with center frequencies of 0.1 Hz, 0.2 Hz, 0.4 Hz, 0.8 Hz, 1.6 Hz, 3.2 Hz, 6.4 Hz, and 12.8 Hz respectively, and the bandwidth is ±15% of the center frequency. When performing a joint decomposition on the registered image data, use an image block of 256×256 pixels as the processing unit, and perform a convolution operation on each image block with all filters to obtain an initial decomposition result.

[0022] When constructing shearlet basis functions according to shear parameters and translation parameters, the value range of shear parameters is [-2.0, 2.0], the step size is 0.1, the translation parameters are uniformly sampled in the image space, and 32 sampling points are taken in both the x - direction and the y - direction. For each group of parameter combinations, the corresponding shearlet basis functions are constructed. Specifically, the shearlet basis functions are in the form of Gabor functions, which include a real part and an imaginary part. The real - part function is expressed as a cosine function modulated by a Gaussian function, and the imaginary - part function is expressed as a sine function modulated by a Gaussian function. The standard deviation of the Gaussian envelope function is set to 1 / 8 of the image block size.

[0023] When obtaining the coefficient matrix through two - dimensional integral operation on the shearlet basis functions, the discrete summation method is used to approximate the integral process. For each 256×256 image block, its dot - product with the shearlet basis function is calculated to obtain the corresponding coefficient value. For all parameter combinations, a coefficient matrix is formed, and the matrix dimension is (the number of shear parameters×the number of translation parameters×the number of directions).

[0024] When calculating the direction - angle parameter based on the arctangent function of the shear parameter, for the shear parameter s, the direction angle θ is calculated by θ = arctg(s). When s=-2.0, θ is approximately - 63.4 degrees; when s = 0, θ is 0 degrees; when s = 2.0, θ is approximately 63.4 degrees. To optimize the direction - angle parameter, an adaptive threshold mechanism is introduced. According to the energy distribution of the coefficient matrix, the direction angles with energy contribution exceeding 5% of the total energy are retained to achieve multi - direction adaptive decomposition.

[0025] The temperature - field coefficient and the acoustic coefficient are obtained through processing. The temperature - field coefficient corresponds to the infrared thermal imaging modality, and the acoustic coefficient corresponds to the ultrasonic modality. For typical trauma bleeding cases, the temperature - field coefficient shows higher low - frequency components and energy concentration in specific directions in the bleeding area, and the acoustic coefficient shows echo enhancement and tissue - structure change characteristics in the bleeding area.

[0026] Substitute the temperature - field coefficient into the heterogeneous heat - diffusion equation for solution. The tissue density value range is 900 - 1100 kg / m³, the density of the bleeding area is about 1060 kg / m³, and the density of normal tissue is about 1020 kg / m³. The specific - heat - capacity term value range is 3500 - 4200 J / (kg·℃), the specific heat capacity of the bleeding area is about 3900 J / (kg·℃), and the specific heat capacity of normal tissue is about 3700 J / (kg·℃).

[0027] Calculate the product of the tissue density term and the specific heat capacity term to obtain the coefficient term. Taking the bleeding area as an example, the coefficient term is 1060×3900 = 4,134,000 J / (m³·℃). The thermal conductivity tensor term is represented by a 2×2 matrix, and the range of the thermal conductivity values in the principal axis direction is 0.3 - 0.6 W / (m·℃). The principal axis thermal conductivity of the bleeding area is approximately [0.53, 0.42] W / (m·℃), and the thermal conductivity in the non-principal axis direction is approximately 0.1 W / (m·℃).

[0028] Calculate the inner product of the thermal conductivity tensor term and the temperature gradient to obtain the diffusion term. The temperature gradient is calculated by the central difference method, and the spatial step size is 1 mm. At the boundary of the bleeding area, the typical temperature gradient is approximately 0.5 - 1.5 ℃ / mm. The heat source term is set according to physiological characteristics. The heat source value of normal tissue is approximately 420 W / m³, and the heat source of the bleeding area is enhanced due to blood extravasation, with a value of approximately 650 W / m³.

[0029] When using the finite difference method to solve the heterogeneous heat diffusion equation, the time step size is set to 0.1 s, the spatial step size is 1 mm, and the explicit format is used for iterative calculation. The boundary condition is set as the third type of boundary condition, the surface heat transfer coefficient is taken as 12 W / (m²·℃), and the ambient temperature is set to 25℃. After 100 time steps of iteration, the temperature field reaches a steady state, and the deep temperature distribution map is obtained.

[0030] Calculate the anisotropic thermal diffusion coefficient based on the ratio of the maximum thermal conductivity to the minimum thermal conductivity. Taking the bleeding area as an example, the anisotropic thermal diffusion coefficient is 0.53÷0.42 = 1.26, indicating that there are directional differences in heat conduction in the bleeding area. The anisotropic thermal diffusion coefficient of normal tissue is approximately 1.05 - 1.15, and that of the bleeding area is approximately 1.20 - 1.35.

[0031] Perform pulse compression processing on the acoustic coefficient to obtain the compressed echo signal. Use a 13-bit Barker code as the modulation sequence, and the compression gain is approximately 10.8 dB. The compression processing uses the matched filtering technique, and the filter coefficient is the time-reversed conjugate form of the modulation sequence. For clinical data, the compressed axial resolution is improved to 0.3 mm.

[0032] Perform phase correction by calculating the difference between the original phase and the reference phase. The reference phase is obtained by statistical averaging of the healthy tissue area. For the bleeding area, the phase difference is usually 15 - 40 degrees, and the corrected phase error is controlled within ±5 degrees.

[0033] Calculate the acoustic reflection coefficient according to the ratio of the difference to the sum of the acoustic impedances of adjacent tissues. The acoustic impedance of normal soft tissue is approximately 1.60×10 6 kg / (m²·s), and the acoustic impedance of the bleeding area is approximately 1.66×10 6The acoustic reflection coefficient at the interface between the two is approximately 0.018. The acoustic reflection coefficient typically varies between 0.001 and 0.05. A bleeding area can be identified if the difference between the acoustic reflection coefficient of the bleeding area and the surrounding tissue is greater than 0.01.

[0034] A tissue layer map was reconstructed based on the acoustic reflection coefficient, compressed echo signals, and phase correction results. A layer-stripping algorithm was used in the reconstruction process, calculating the reflected signal intensity layer by layer from superficial to deep. The synthetic aperture was 32 elements, with a focal depth from the surface to 60 mm and a step size of 5 mm. The reconstructed tissue layer map had a resolution of 0.3 mm × 0.5 mm, clearly demonstrating the boundaries and internal structures of the hemorrhage area.

[0035] In this embodiment, a direction-selective filter group and optimized directional angle parameters are used to accurately extract features in different directions in the image, thereby improving the spatial resolution and directional discrimination of shear wave decomposition; the temperature field coefficient is substituted into the heterogeneous heat diffusion equation and solved by the finite difference method, so that a high-precision deep temperature distribution map can be obtained inside the tissue, which is beneficial for temperature control monitoring in medical applications such as thermal therapy and ablation; based on the ratio of maximum thermal conductivity to minimum thermal conductivity, the anisotropic thermal diffusion coefficient is calculated, which can quantitatively describe the directional dependence of tissue heat conduction and provide a basis for the design of personalized thermal therapy plans; the acoustic coefficients are pulse compressed, phase corrected, and the reflection coefficient is calculated, which greatly improves the signal-to-noise ratio and phase accuracy of the echo signal, which contributes to the clarity and reliability of subsequent tissue imaging; the compressed echo and phase correction results are combined, and tissue tomography reconstruction is performed based on the acoustic impedance difference, which can accurately depict the tissue distribution at different acoustic levels and achieve fine visualization of subtle structures.

[0036] In an optional embodiment, solving the heterogeneous heat diffusion equation using a finite difference method to obtain a deep temperature distribution map includes: Discretize the time and space coordinates, divide the time dimension into multiple time layers according to the time step, divide the space dimension into a three-dimensional grid structure according to the space step, construct a differential grid, and obtain the initial temperature value at the grid point of the differential grid; performing differential discretization processing on the heterogeneous heat diffusion equation, discretizing the time derivative term using a first-order forward difference format, discretizing the space derivative term using a second-order central difference format, interpolating the thermal conductivity tensor at half-grid points to obtain a discrete thermal conductivity value, and substituting the discrete time derivative term, the discrete space derivative term, and the discrete thermal conductivity value into the heterogeneous heat diffusion equation to establish a differential equation; Constructing boundary conditions at the boundary positions of the differential grid, discretizing the boundary conditions of the temperature boundary, the heat flow boundary, and the convection boundary, and combining the boundary conditions with the differential equations to construct a complete differential equation group; The alternating direction implicit method is used to decompose and solve the complete difference equations in the x-direction, y-direction, and z-direction in sequence. A tridiagonal system of equations is constructed for each direction and the iterative solution is performed using the chasing method to obtain the temperature distribution at the new time level. When the temperature difference between adjacent time levels is less than the preset convergence threshold, a deep temperature distribution map is generated based on the finally converged temperature distribution.

[0037] In a specific embodiment, when discretizing the spatio-temporal coordinates, the discretization of the time dimension is carried out in a fixed time step manner. In this embodiment, the time step is set to 0.1 s, the total simulation duration is 10 s, and a total of 100 time levels are planned to be divided. The discretization of the spatial dimension is carried out in a structured grid manner, and the spatial step sizes in the x-direction, y-direction, and z-direction are set to 0.5 mm, 0.5 mm, and 0.5 mm respectively. For a typical traumatic bleeding detection scenario, the coverage range of the constructed three-dimensional difference grid is 100 mm × 100 mm × 60 mm, corresponding to 201×201×121 grid points.

[0038] When obtaining the initial temperature value at the grid points of the difference grid, for the grid points on the skin surface, the temperature value is directly obtained through an infrared thermal imager, and the typical temperature range is 32°C to 36°C. For the grid points of the subcutaneous tissue, the initial temperature is estimated by an empirical formula. The temperature at 1 cm below the surface is about 36.5°C, the temperature at 2 cm is about 37.0°C, and the temperature at 3 cm and deeper is close to 37.2°C. For the bleeding area, the initial temperature is set to 37.3°C to 37.5°C, slightly higher than the temperature of the surrounding normal tissue.

[0039] When performing difference discretization on the heterogeneous heat diffusion equation, the time derivative term is discretized using the first-order forward difference format. Specifically, the derivative of temperature with respect to time is approximated as the difference between the temperature at the current time level and the next time level divided by the time step. In this embodiment, for the grid point (i, j, k), the discrete term of the time derivative is expressed as the temperature at the next time level minus the temperature at the current time level, and then divided by 0.1 s.

[0040] The spatial derivative term is discretized using the second-order central difference format. For the second-order derivative term, the temperature values of three adjacent grid points are used for calculation. In the x-direction, the discrete term of the second-order derivative is expressed as the temperature value of the grid point (i + 1, j, k) plus the temperature value of the grid point (i - 1, j, k), minus 2 times the temperature value of the grid point (i, j, k), and then divided by the square of the x-direction step size. The second-order derivative discretization methods in the y-direction and z-direction are similar. For the mixed derivative term, the temperature values of four adjacent diagonal grid points are used for calculation.

[0041] The discrete value of thermal conductivity is obtained by interpolating the thermal conductivity tensor at the semi-grid points. The thermal conductivity tensor contains nine components, corresponding to three principal directions and their cross terms. In this embodiment, the values of the main diagonal elements of the thermal conductivity tensor in the bleeding area are 0.52 W / (m·°C), 0.48 W / (m·°C), and 0.45 W / (m·°C), and the values of the off-diagonal elements are 0.08 W / (m·°C), 0.06 W / (m·°C), and 0.07 W / (m·°C). The values of the main diagonal elements of the thermal conductivity tensor in normal tissue are 0.45 W / (m·°C), 0.43 W / (m·°C), and 0.42 W / (m·°C), and the values of the off-diagonal elements are approximately 0.02 W / (m·°C).

[0042] For the thermal conductivity at the grid point (i + 1 / 2, j, k), the arithmetic mean of the thermal conductivities at the grid points (i, j, k) and (i + 1, j, k) is used. If the grid point (i, j, k) is located at the boundary of the bleeding area, the thermal conductivity value is taken as the weighted average of the thermal conductivities of adjacent grid points, and the weight coefficient is inversely proportional to the distance of each grid point to the boundary. Through this interpolation method, the smooth transition of thermal conductivity at the tissue boundary is ensured, and the instability of numerical calculation is avoided.

[0043] The time derivative discrete term, the spatial derivative discrete term, and the discrete value of thermal conductivity are substituted into the heterogeneous heat diffusion equation to establish a difference equation. In this embodiment, the left side of the difference equation is the product of the time derivative discrete term and the product of density and specific heat capacity. The density value is 1060 kg / m 3 , and the specific heat capacity value is 3900 J / (kg·°C). The right side is the discrete term of the inner product of the discrete value of thermal conductivity and the temperature gradient, plus the heat source term. The heat source term takes the value of 650 W / m 3 [[ID=~10]] in the bleeding area and 420 W / m 3 in the normal tissue area.

[0044] Boundary conditions are constructed at the boundary positions of the difference grid. For the temperature boundary, the temperature values of the boundary grid points are directly specified. In this embodiment, the temperature boundary on the skin surface is set according to the infrared thermal imaging measurement results, and the temperature value range is 32°C to 36°C. For the heat flux boundary, the heat flux density on the boundary surface is specified. On the deep tissue boundary surface, the heat flux density is set to 20 W / m 2 , and the direction points to the inside of the body. For the convective boundary, the heat flux is calculated according to the ambient temperature and the heat transfer coefficient. On the skin surface, the ambient temperature is set to 25°C, and the heat transfer coefficient is set to 12 W / (m 2 ·°C).

[0045] The boundary conditions are combined with the difference equations to construct a complete system of difference equations. For interior grid points, standard difference equations are applied. For boundary grid points, the coefficients and constants of the difference equations are modified based on the various boundary conditions. For example, for grid points on convective boundaries, an additional term is added to the difference equations representing the product of the heat transfer coefficient and the difference between the ambient and surface temperatures. The complete system of difference equations forms a large, sparse linear system.

[0046] The complete system of difference equations is solved using the alternating direction implicit method. Within each time step, the solution is performed in three substeps. In the x-direction substep, a tridiagonal system of equations is constructed, with coefficients consisting of implicit terms in the x direction and explicit terms in the y and z directions. In the y-direction substep, a tridiagonal system of equations is constructed, with coefficients consisting of implicit terms in the y direction and explicit terms in the x and z directions. In the z-direction substep, a tridiagonal system of equations is constructed, with coefficients consisting of implicit terms in the z direction and explicit terms in the x and y directions.

[0047] The chase method is used to iteratively solve the tridiagonal system of equations in each direction. This method solves the tridiagonal system through two processes: forward and backward iteration. During the forward iteration, intermediate variables are calculated and the lower diagonal elements are eliminated. During the backward iteration, the temperature values of each grid point are solved sequentially from back to front. After completing the iterations in all three directions, the temperature distribution of the new time layer is obtained.

[0048] When the temperature difference between adjacent time layers is less than the preset convergence threshold, the temperature field is considered to have reached a steady state. In this embodiment, the preset convergence threshold is set to 0.01°C, that is, the iteration is stopped when the maximum temperature difference between the corresponding grid points in adjacent time layers is less than 0.01°C. Based on the final converged temperature distribution, a deep temperature distribution map is generated. The temperature distribution map uses a pseudo-color display method, with a temperature range from 32°C to 40°C and a total of 16 color levels, each corresponding to a temperature interval of 0.5°C.

[0049] Traditional methods typically use explicit difference formats, which require very small time steps to ensure numerical stability and have low computational efficiency. Some methods use implicit formats, but their processing of the thermal conductivity tensor is not precise enough, and numerical oscillations are easily generated at the tissue interface. The method of this embodiment uses the alternating direction implicit method, which not only ensures numerical stability but also improves computational efficiency. At the same time, by accurately interpolating the thermal conductivity tensor at half-grid points, numerical oscillations at the tissue interface are effectively eliminated. These improvements increase the temperature field reconstruction accuracy from the original ±0.5°C to ±0.1°C, and shorten the calculation time from several minutes to more than ten seconds, providing a reliable guarantee for the rapid and accurate detection of traumatic bleeding.

[0050] like Figure 2As shown, the computational times of three different methods at various grid sizes are presented on a logarithmic scale. The alternating direction implicit method of the present technical solution shows significant performance advantages at all grid sizes. At the largest 250 3 grid, the explicit difference method requires 1734.5 seconds, the standard implicit method requires 912.8 seconds, while the present technical solution only requires 229.4 seconds, which is 7.6 times and 4.0 times faster than the former two respectively. Notably, as the grid size increases, the performance advantage of the present technical solution becomes more obvious, indicating that this method has significant computational efficiency advantages in dealing with high-resolution three-dimensional temperature field reconstruction. This high efficiency enables near-real-time temperature field analysis in a clinical environment.

[0051] In an optional implementation manner, according to the deep temperature distribution map and the tissue stratification map, a complex-valued neural network is used to construct a thermoacoustic coupling feature matrix. Based on the thermoacoustic coupling feature matrix, the blood flow velocity vector is extracted from the Doppler shift signal, and the Marching Cubes interpolation is used to reconstruct the vascular structure. Generating the blood flow velocity field distribution map includes: Normalize the deep temperature distribution map and use Gaussian filtering to obtain the temperature field gradient. Perform adaptive threshold segmentation on the tissue stratification map to obtain the tissue interface normal vector. Combine the temperature field gradient and the tissue interface normal vector to form an initial feature vector; Input the initial feature vector into the complex-valued neural network. The complex-valued neural network adopts a three-layer network structure and has a piecewise linear unit function as the activation function. Process the real part and the imaginary part of the initial feature vector through the complex-valued neural network to obtain the thermoacoustic coupling feature matrix; Obtain the complex analytic signal by performing Hilbert transform on the Doppler shift signal. Perform tensor product operation on the complex analytic signal and the thermoacoustic coupling feature matrix to obtain a three-dimensional blood flow velocity vector field; Calculate the velocity divergence distribution according to the three-dimensional blood flow velocity vector field, determine the vascular boundary point set, and input the vascular boundary point set into the cubic Marching Cubes interpolation algorithm to obtain the vascular surface; Within the range of the vascular surface, use the radial basis function to perform interpolation calculation on the three-dimensional blood flow velocity vector field to obtain a continuous blood flow velocity field distribution map.

[0052] In a specific implementation manner, obtain the patient's deep temperature distribution map and tissue stratification map. Among them, the deep temperature distribution map can be obtained by a thermoacoustic imaging device and is represented as a two-dimensional image of 256×256 pixels, and the pixel value range is 20°C to 42°C; the tissue stratification map is obtained by a medical tomography device and is represented as a two-dimensional image of the same size, where different tissue types are represented by different gray values.

[0053] Normalize the deep temperature distribution map and map the temperature values to the range of 0 to 1. The specific implementation method is that for the original temperature value T(x,y), calculate the normalized temperature value through (T(x,y)-Tmin) / (Tmax-Tmin), where Tmin = 20°C and Tmax = 42°C. Apply Gaussian filtering to the normalized temperature distribution map and perform a convolution operation using a Gaussian kernel with a standard deviation of 1.5 to obtain a smoothed temperature field. Then calculate the temperature field gradient, calculate the gradient components in the x and y directions respectively, and form a gradient vector field.

[0054] Perform adaptive threshold segmentation on the tissue layer map and automatically determine the optimal threshold using the OTSU method. In this embodiment, for different tissue types such as muscle, fat, and blood vessels, the thresholds are determined to be 85, 120, and 165 (gray value range 0 - 255) respectively. After segmentation, obtain a binary tissue interface map, then use the Sobel operator to detect the edges of the tissue interface, calculate the normal vectors at the edges, and obtain the tissue interface normal vector field.

[0055] Combine the temperature field gradient and the tissue interface normal vector to form an initial feature vector. For each pixel point (x,y) in the image, extract the temperature gradient vector (Gx,Gy) and the tissue interface normal vector (Nx,Ny), and combine them to form a four-dimensional feature vector (Gx,Gy,Nx,Ny). To convert it into a complex-valued representation, use the temperature gradient as the real part and the tissue interface normal vector as the imaginary part to form a complex-valued feature vector (Gx + jNx, Gy + jNy).

[0056] Construct a complex-valued neural network to process the above initial feature vector. The network adopts a three-layer structure, including an input layer, a hidden layer, and an output layer. The number of nodes in the input layer is 2, corresponding to the two components of the complex-valued feature vector; the number of nodes in the hidden layer is 16; the number of nodes in the output layer is 4, corresponding to the elements of the output feature matrix. The network uses a piecewise linear unit function as the activation function, which is defined as: when the real part or imaginary part of the input is less than 0, the output is 0; when the input is greater than or equal to 0 and less than 1, the output is equal to the input value; when the input is greater than or equal to 1, the output is 1.

[0057] The network training uses the stochastic gradient descent method with a batch size of 64, the initial learning rate is set to 0.01, and the number of training epochs is 200. Use the mean squared error as the loss function and stop training when the error on the validation set reaches below 0.05. After training is completed, input the initial feature vector into the network to obtain a thermoacoustic coupling feature matrix, which is a 2×2 complex-valued matrix representing the thermoacoustic coupling coefficient.

[0058] Process the Doppler frequency shift signal. The collected Doppler frequency shift signal (sampling rate: 10 kHz, number of sampling points: 1024) is transformed into a complex analytic signal through Hilbert transform. In specific implementation, first perform Fourier transform on the original frequency shift signal s(t) to obtain the frequency spectrum, then set the negative frequency components to zero, and finally perform inverse Fourier transform to obtain the complex analytic signal s_a(t) = s(t) + j·H[s(t)], where H[s(t)] is the Hilbert transform of s(t).

[0059] Perform a tensor product operation on the complex analytic signal and the thermoacoustic coupling characteristic matrix. For each time point t and spatial position (x, y), calculate the product of the complex analytic signal s_a(t) and the thermoacoustic coupling characteristic matrix at this position to obtain the three-dimensional blood flow velocity vector field v(x, y, t).

[0060] Calculate the velocity divergence distribution based on the three-dimensional blood flow velocity vector field, that is, calculate the divergence value of the velocity vector for each point (x, y). Regions with higher velocity divergence values usually correspond to blood vessel boundaries. Set the divergence threshold to 0.15, and extract the set of points with divergence values greater than the threshold as the blood vessel boundary point set. In this embodiment, approximately 2500 boundary points are obtained.

[0061] Input the blood vessel boundary point set into the cubic Marching Cubes interpolation algorithm to construct a smooth blood vessel surface. The algorithm inserts control points between adjacent boundary points to ensure that the generated surface has C2 continuity. Specifically, set the tension parameter to 0.5 and the deviation parameter to 0.0 to obtain a smooth blood vessel surface model, which contains approximately 15000 surface meshes.

[0062] Within the range of the blood vessel surface, use the radial basis function to perform interpolation calculation on the three-dimensional blood flow velocity vector field. Select the Gaussian radial basis function, whose form is exp(-r 2 / σ 2 ), where r is the spatial distance and σ is the shape parameter, with a value of 3.0. For any point within the blood vessel, calculate the velocity value of this point based on the weighted average of the surrounding known velocity vectors, and the weights are determined by the radial basis function. After interpolation, a continuous blood flow velocity field distribution map with a spatial resolution of 0.1 mm is generated, clearly showing the direction and magnitude distribution of the blood flow within the blood vessel.

[0063] In this embodiment, the whole process of generating a high-precision blood flow velocity field distribution map is successfully realized starting from the deep temperature distribution map and the tissue stratification map, via constructing a thermoacoustic coupling feature matrix by a complex-valued neural network. By fusing the normalized temperature gradient with the tissue interface normal vector and processing it by a three-layer complex-valued neural network, the weak signal features of thermoacoustic coupling can be effectively captured, enhancing the integrity and distinctiveness of feature expression; performing a tensor product operation on the complex-analytic Doppler signal and the thermoacoustic feature matrix can directly output a three-dimensional blood flow velocity vector field, realizing a high-precision description of complex hemodynamics; determining the vascular boundary point set through the velocity divergence distribution and introducing the cubic Marching Cubes interpolation algorithm can reconstruct a smooth and continuous three-dimensional vascular surface model without manual intervention; using a radial basis function to interpolate the velocity vector field within the vascular surface range to generate a continuous and smooth blood flow velocity field distribution map is helpful for subsequent hemodynamic analysis and simulation; the full-link algorithm design from the temperature field and tissue stratification to the blood flow field interpolation realizes seamless data transfer and processing, improving the automation degree of the system and the robustness to noise / variations.

[0064] In an alternative embodiment, calculating the velocity divergence distribution according to the three-dimensional blood flow velocity vector field, determining the vascular boundary point set, and inputting the vascular boundary point set into the cubic Marching Cubes interpolation algorithm to obtain the vascular surface includes: Performing scale pyramid decomposition on the three-dimensional blood flow velocity vector field to obtain velocity vector fields at multiple scale levels, respectively calculating the velocity divergence distribution of the velocity vector fields to determine the multi-scale divergence distribution, and based on the multi-scale divergence distribution, performing feature enhancement through average pooling and a multi-layer perceptron to obtain an enhanced velocity divergence distribution; Constructing the enhanced velocity divergence distribution into a graph structure including a vertex set, an edge set, and an adjacency matrix, inputting the graph structure into a graph attention network for iterative feature update to obtain vascular boundary features; Extracting the skeleton vascular boundary point set based on the vascular boundary features and performing recursive refinement to obtain the complete vascular boundary point set; Constructing a joint loss function including a topological consistency constraint term and a geometric constraint term, optimizing the complete vascular boundary point set based on the joint loss function to determine the optimized vascular boundary point set; Dividing the optimized vascular boundary point set into local surface patches and establishing a parameterized grid, calculating the control vertices of the parameterized grid, constructing cubic B-spline basis functions based on the control vertices, performing a tensor product operation with the position information of the control vertices to obtain the local surface reconstruction result, and after imposing curvature continuity constraints and boundary continuity constraints on the local surface reconstruction result, splicing to obtain the vascular surface.

[0065] In a specific embodiment, when performing scale pyramid decomposition on the three-dimensional blood flow velocity vector field, a Gaussian filter and downsampling operations are used to construct a multi-scale representation. In this embodiment, 5 scale levels are constructed, and the scale parameters are 1.0, 2.0, 4.0, 8.0, and 16.0 in sequence. The corresponding spatial resolutions are 1 times, 1 / 2 times, 1 / 4 times, 1 / 8 times, and 1 / 16 times of the original resolution respectively. For a typical traumatic bleeding detection scenario, the spatial resolution of the original blood flow velocity vector field is 0.1 mm, and the spatial resolutions of each layer of the scale pyramid are 0.1 mm, 0.2 mm, 0.4 mm, 0.8 mm, and 1.6 mm in sequence. Through multi-scale decomposition, large-scale vascular structures and small-scale vascular branch information can be captured simultaneously.

[0066] When calculating the velocity divergence distribution for the velocity vector field respectively, the divergence value represents the amount of fluid flowing out of the unit volume net. A positive value indicates outflow (source), and a negative value indicates inflow (sink). In the traumatic bleeding area, since blood flows from blood vessels to tissues, it usually shows a relatively large positive divergence value. For the velocity vector field of each scale level, calculate the velocity differences of adjacent grid points in three directions and divide by the grid spacing in the corresponding direction to obtain the velocity divergence value. In this embodiment, the divergence value range in the normal vascular area is -0.05 to 0.05 s -1 , while the divergence value in the bleeding area is usually greater than 0.2 s -1 , and can reach up to 2.0 s -1 .

[0067] After determining the multi-scale divergence distribution, feature enhancement is performed through average pooling and a multi-layer perceptron. The average pooling operation uses a 3×3×3 pooling window with a stride of 1 to smooth the divergence distribution and reduce the influence of noise. The multi-layer perceptron contains three layers of neural networks. The number of nodes in the input layer is the same as the number of pixels in the divergence image, the number of nodes in the hidden layer is 64, and the number of nodes in the output layer is the same as that in the input layer. The activation function uses the ReLU function to improve the non-linear expression ability of the network. During the training process, the learning rate is set to 0.001, the Adam optimizer is used, and the number of training epochs is 100. Through feature enhancement, the contrast of the divergence distribution is improved, and the distinguishability between the bleeding area and the normal area is enhanced, obtaining an enhanced velocity divergence distribution.

[0068] When constructing the enhanced velocity divergence distribution into a graph structure, each voxel point serves as a vertex of the graph. The vertex feature vector contains the divergence value of this point and its divergence values at each scale, with a total of 6 components. The edges of the graph connect spatially adjacent voxel points, and the neighborhood radius is set to 2 voxel distances. On average, each vertex is connected to about 26 neighbor vertices. The adjacency matrix is a sparse matrix, and the non-zero element value is the similarity between vertices. The calculation method is the exponential function of the cosine similarity of the vertex feature vectors, and the similarity range is from 0 to 1.

[0069] The graph structure is input into the graph attention network for feature iterative update. The graph attention network contains 3 layers of graph convolutional layers. Each layer of graph convolution uses the multi-head attention mechanism, with the number of heads set to 8, the input feature dimension to 6, the hidden layer feature dimension to 32, and the output feature dimension to 16. The attention weights are normalized by the Softmax function to ensure that the sum of the weights of each neighbor vertex is 1. After each layer of graph convolution, a residual connection and batch normalization are added to improve the network stability. Through the feature propagation of 3 layers of graph convolution, each vertex synthesizes the information within its local neighborhood, more accurately representing the vascular boundary features. For the bleeding area, due to the obvious differences in its features from the normal vascular area, a unique representation pattern will be formed during the feature update process, which helps with subsequent boundary extraction.

[0070] When extracting the skeleton vascular boundary point set based on the vascular boundary features, a threshold segmentation and thinning algorithm is adopted. First, the feature map is binarized, and the threshold is set to the feature mean plus 0.8 times the standard deviation to obtain the initial vascular area mask. The distance transform is applied to the mask to obtain the shortest distance from each point to the background. The local maximum points of the distance transform are marked as skeleton points. In this embodiment, a typical skeleton vascular network contains 500 - 1000 skeleton points, which varies according to the complexity of the vascular network.

[0071] The complete vascular boundary point set is obtained by recursively refining the skeleton vascular boundary point set. The recursive refinement adopts an adaptive grid encryption strategy. The distance between the initial boundary points is about 1mm, and encryption is performed in areas with larger curvature, with the minimum inter-point distance reaching 0.1mm. During the refinement process, the positions of the newly added points are determined by cubic spline interpolation and adjusted considering the local curvature information. The recursive refinement is executed 5 times, increasing the boundary point density by about 10 times in the high-curvature area, and finally obtaining a complete vascular boundary point set containing 5000 - 10000 points.

[0072] A joint loss function containing a topological consistency constraint term and a geometric constraint term is constructed to optimize the complete vascular boundary point set. The topological consistency constraint term ensures the connectivity and branch structure of the vascular network remain unchanged, with the weight coefficient set to 0.6. The geometric constraint term includes a curvature smoothing constraint and a volume preservation constraint, with weight coefficients of 0.3 and 0.1 respectively. The curvature smoothing constraint calculates the discrete curvature using the Laplace operator to limit the gradient magnitude of the curvature change. The volume preservation constraint ensures that the overall volume change of the blood vessels before and after optimization does not exceed 3%. The optimization algorithm uses the L - BFGS method, with a maximum number of iterations of 200 times and a convergence threshold of 0.0001. Through optimization, while maintaining the topological structure, the boundary point set achieves smoothness and accuracy in geometric shape.

[0073] The optimized vascular boundary point set is divided into local surface patches and a parametric grid is established. The division of local surface patches is based on the region growing method. Bounded by vascular branch points, the vascular network is decomposed into multiple tubular segments. In this embodiment, a typical vascular network is divided into 15 - 25 local surface patches. For each local surface patch, a parametric grid is established with a grid resolution of 16 dissection points circumferentially, and the number of sampling points axially is adaptively set according to the vascular length, approximately 8 sampling points per millimeter. Parametrization uses cylindrical coordinate mapping to map the three - dimensional space surface to a two - dimensional parametric space.

[0074] The control vertices of the parametric grid are calculated. Each local surface patch is set with a control point grid of 4×4×4, and the initial positions of the control points are obtained by least - square fitting of the point set on the parametric grid. Based on the control vertices, cubic B - spline basis functions are constructed with a function order of 3, and the knot vector uses an open uniform knot vector. The cubic B - spline basis functions are tensor - product - operated with the position information of the control vertices to obtain the parametric representation of the local surface. In this embodiment, each local surface patch is determined by 64 control points, providing sufficient degrees of freedom to represent complex vascular shapes.

[0075] Curvature continuity constraints and boundary continuity constraints are imposed on the local surface reconstruction results. The curvature continuity constraint ensures continuous curvature changes at the connection of adjacent surface patches, technically achieved by matching the first two - order derivatives of the boundaries of adjacent surface patches. The boundary continuity constraint ensures that adjacent surface patches completely coincide at the shared boundary, technically achieved by sharing boundary control points. The two constraints are combined into the objective function by the Lagrange multiplier method, and the weight coefficients are both 0.5. After constraint optimization, all local surface patches are stitched together to obtain a complete and smooth vascular surface reconstruction result.

[0076] In the prior art, Marching Cubes or direct triangular mesh methods are mainly used for vascular surface reconstruction, which cannot effectively handle complex topological changes in bleeding areas. The method of this embodiment starts from hemodynamic characteristics, extracts vascular boundaries using the divergence field, enhances boundary features through a graph attention network, and combines geometric constraints and topological constraints to reconstruct a smooth and continuous vascular surface. This method improves the boundary accuracy by 30% compared with traditional techniques, increases the bleeding area detection rate by 25%, and reduces the average error of the reconstructed surface from 0.5 mm to 0.15 mm. Especially for the reconstruction of complex topological structures around small vascular branches and bleeding points, the method of this embodiment shows significant advantages, providing reliable technical support for precise clinical positioning of traumatic bleeding.

[0077] As Figure 3As shown, the surface curvature continuity of blood vessels reconstructed by three different methods was evaluated. The horizontal axis represents the normalized position (0 - 1) of points on the blood vessel surface, and the vertical axis represents the rate of change of mean curvature. For the technical solution (blue line), the rate of change of curvature is basically stable between 0.02 - 0.08 across the entire surface, with an average of 0.043. The curve is continuous and smooth, and the fluctuation is only 0.032 in the key area (the branch intersection corresponding to the normalized position of 0.4 - 0.6), indicating that the reconstructed surface has a high degree of geometric continuity. The rate of change of curvature of the Marching Cubes method (red line) has an average of 0.176, with a spike reaching 0.412 near the branch intersection (around 0.5), indicating obvious geometric discontinuity. The average rate of change of curvature of the direct triangular mesh method (green line) is 0.231, with irregular fluctuations at multiple positions, and the maximum fluctuation reaches 0.538 (at position 0.53), and there are obvious step-like changes at the surface connection points (at positions 0.25 and 0.75). The mean square deviations of the rate of change of curvature of the three methods are 0.0019 (for the technical solution), 0.0187 (for the Marching Cubes method), and 0.0226 (for the direct triangular mesh method), respectively, indicating that the technical solution effectively ensures the smooth transition of the surface, reduces geometric artifacts, and is particularly suitable for medical application scenarios that require high-precision hydrodynamic analysis by using cubic Marching basis functions and imposing curvature continuity constraints.

[0078] In an alternative embodiment, a variational particle filtering method is used to dynamically track the local blood flow rate in the time-sequence mapping relationship, calculate the corresponding time-varying coefficients, and construct a state equation to obtain bleeding dynamic parameters including: Initialize multiple particles according to the time-sequence mapping relationship and assign initial weights, where each particle corresponds to a state of the local blood flow rate; Establish an observation likelihood function to describe the relationship between the cumulative blood loss and the local blood flow rate, and establish a state transition probability function to describe the change law of the local blood flow rate; combine the observation likelihood function and the state transition probability function to construct a variational objective function; Update the weights of the particles by iterative optimization based on the variational objective function; calculate the optimal estimated value of the local blood flow rate according to the updated particle weights; perform importance resampling to eliminate particle degeneracy and obtain the dynamic change result of the local blood flow rate; Analyze the change law of the blood flow rate in different time periods according to the dynamic change result of the local blood flow rate; calculate the scale coefficient, delay coefficient, and attenuation coefficient of the hemodynamic characteristics based on the change law, and perform time-sequence smoothing processing to determine each smoothing coefficient; Construct a state transition matrix of a non-linear dynamics model using a smoothing coefficient, and combine the system noise term with the state transition matrix to construct a state equation for the bleeding dynamic characteristics; solve the state equation for the bleeding dynamic characteristics, calculate the bleeding rate, the cumulative blood loss, and the degree of tissue damage, and determine the bleeding dynamic parameters.

[0079] In a specific embodiment, a plurality of particles are initialized based on a time series mapping relationship and initial weights are assigned. Assuming that the local blood flow state is X, 200 particles are initialized, and each particle represents a possible blood flow state. For each particle P(i), where i = 1, 2,..., 200, the same initial weight w(i) = 1 / 200 is assigned. For example, if the initially estimated average blood flow rate is 5 ml / min, particles obeying a normal distribution can be generated around this value, with a standard deviation set to 1 ml / min.

[0080] An observation likelihood function and a state transition probability function are established. The observation likelihood function is used to describe the relationship between the cumulative blood loss Z and the local blood flow X. In practical applications, a Gaussian distribution can be used for modeling, where the error between the observed value and the predicted value follows a normal distribution with a mean of 0 and a standard deviation of 2 ml. For the state transition probability function, a first-order Markov process is used to describe the variation law of the blood flow rate, that is, the blood flow state at the current moment is only related to the state at the previous moment. For example, if the blood flow rate at the previous moment is 4.5 ml / min, the blood flow rate at the current moment may fluctuate within the range of 4.3 - 4.7 ml / min.

[0081] Based on the above two functions, a variational objective function is constructed. This function combines the observation likelihood probability and the state transition probability, and is used to evaluate the credibility of each particle. The design of the variational objective function takes into account the balance between data fitting and model prior knowledge, and optimizes the particle distribution by minimizing the information entropy difference.

[0082] The particle weights are updated through iterative optimization. At each time step t, the new weight of each particle is calculated based on the variational objective function. In a specific implementation, 20 iterations are used, and in each iteration, the particle weights are adjusted by the gradient descent method, with the learning rate set to 0.05. For example, if the blood flow rate corresponding to a certain particle at t = 10 min is 6.2 ml / min, and the observed data shows that this value is more likely to be 6.0 ml / min, then the weight of this particle will be correspondingly reduced.

[0083] Calculate the optimal estimated value based on the updated particle weights. The weighted average method is used to calculate the estimated value of the local blood flow rate, that is, the sum of all particle values multiplied by their corresponding weights. If there are 5 particle values of 5.8, 6.0, 6.1, 5.9, 6.2 ml / min, and the weights are 0.15, 0.25, 0.3, 0.2, 0.1 respectively, then the optimal estimated value is 6.01 ml / min.

[0084] Perform importance resampling to eliminate particle degeneracy. Resampling is triggered when the effective number of particles is below the threshold of 100. The resampling process generates a new set of particles based on the particle weights, and particles with larger weights have a higher probability of being selected and replicated. For example, a particle with a weight of 0.3 may be replicated into 3 identical particles with a weight of 0.1. The particle set obtained after resampling better represents the actual distribution of local blood flow.

[0085] Based on the dynamic change results of local blood flow, analyze the change patterns in different time periods. The entire observation period is divided into multiple time windows, each window being 15 minutes. For example, within 0 - 15 minutes after trauma, the blood flow may rapidly increase from 2 ml / min to 8 ml / min; it may stabilize at 7 - 8 ml / min within 15 - 30 minutes; and it may gradually decrease to 5 ml / min within 30 - 60 minutes.

[0086] Calculate the hemodynamic characteristic parameters. The scale coefficient represents the blood flow intensity. For example, the maximum blood flow can reach 9 ml / min, corresponding to a scale coefficient of 1.8; the delay coefficient reflects the blood response time. For example, the blood flow starts to change significantly 2 minutes after trauma, corresponding to a delay coefficient of 0.033; the attenuation coefficient represents the blood flow weakening rate. For example, during the period of 40 - 60 minutes, the blood flow decreases from 7 ml / min to 5 ml / min, corresponding to an attenuation coefficient of 0.017.

[0087] Perform time - series smoothing on the above coefficients. The sliding window averaging method is used, with a window size of 5 time points. For example, for the scale coefficient of 1.75 at a certain time point, averaging is performed by combining the values of the two adjacent time points 1.70, 1.72, 1.75, 1.78, 1.76, and the smoothed value is 1.742.

[0088] Construct the state transition matrix of the nonlinear dynamics model using the smoothed coefficients. This matrix describes how the system state transfers from one time point to the next. For example, the state transition matrix may represent the relationship between the current bleeding rate and the bleeding rate at the previous moment as: current bleeding rate = 0.95×previous bleeding rate + system noise. Here, 0.95 is the state transition coefficient, reflecting the memory effect of the system.

[0089] Combine the system noise term with the state transition matrix to construct the state equation of the bleeding dynamic characteristics. The system noise term is set as Gaussian noise with a mean of 0 and a standard deviation of 0.5 ml / min, representing the model uncertainty.

[0090] Solve the state equation of bleeding dynamic characteristics, and calculate the bleeding rate, cumulative blood loss, and degree of tissue damage. For example, 30 minutes after trauma, the bleeding rate is 6.5 ml / min, the cumulative blood loss is 180 ml, and the degree of tissue damage is moderate (corresponding injury index 0.65). By analyzing these parameters at multiple time points, the bleeding situation can be comprehensively evaluated and guide clinical intervention decisions.

[0091] In this embodiment, through particle filtering and importance resampling based on the variational objective function, high-precision and dynamic tracking and estimation of local blood flow are achieved; the mapping relationship between cumulative blood loss and local blood flow is characterized by the observation likelihood function, making the bleeding monitoring results more physiologically meaningful and accurate; the time-series smoothing process is performed on the particle weight update results, and the scale, delay, and attenuation coefficients are extracted to effectively suppress noise interference and ensure the stability and reliability of parameter estimation; based on the smoothing coefficient, a state transition matrix is constructed and combined with system noise to form a complete state equation of bleeding dynamic characteristics, supporting the modeling and prediction of complex pathological changes; by solving the state equation, key indicators such as bleeding rate, cumulative blood loss, and degree of tissue damage are directly obtained, providing a quantitative basis for clinical diagnosis and intervention.

[0092] In an alternative embodiment, updating the weights of the particles by iterative optimization based on the variational objective function includes: Calculate the partial derivative of the variational objective function with respect to the weights of each particle to obtain the sensitivity eigenvalue; extract the sign of the sensitivity eigenvalue and multiply it by the power term of the sensitivity eigenvalue to obtain the weight update direction; Calculate the KL divergence of the particle weight distributions in two adjacent iterations, construct an adaptive step size, multiply the adaptive step size by the weight update direction to determine the weight update amount, and add the weight update amount to the particle weight of the current iteration to obtain the updated particle weight; Calculate the maximum and minimum values of the updated particle weights, and perform renormalization processing on the updated particle weights to obtain the renormalized particle weights; Calculate the Euclidean distance between the particle weights in two adjacent iterations to obtain the stability index. When the stability index is less than the preset stability threshold for a preset number of consecutive times, determine the particle weight of the current iteration as the optimization result.

[0093] In a specific embodiment, when calculating the sensitivity eigenvalue by taking the partial derivative of the variational objective function with respect to the weights of each particle, the variational objective function consists of two parts: an image similarity metric term and a spatial smoothing regularization term. The image similarity metric term uses normalized mutual information, and its numerical range is usually from 0.5 to 1.5. The mutual information value in the bleeding area is about 0.8, and the mutual information value in the normal tissue area is about 1.2. The spatial smoothing regularization term uses a second-order gradient penalty term, and the weight coefficient is set to 0.3. For a typical traumatic bleeding detection scenario, the number of particles is set to 1024, distributed in a detection area of 100mm×100mm, and the particle spacing is about 3mm. The partial derivative of the variational objective function with respect to the weights of each particle is calculated by the finite difference method, and the central difference step size is set to 0.01. For each particle, by increasing and decreasing its weight by 0.01 respectively, the objective function values are calculated, and the difference between the two divided by 0.02 is the sensitivity eigenvalue of the particle weight. The numerical range of the sensitivity eigenvalue is usually from -2.0 to 2.0. A positive value indicates that increasing the weight of the particle will increase the objective function value, and a negative value indicates that decreasing the weight of the particle will increase the objective function value.

[0094] Extract the sign of the sensitivity eigenvalue and multiply it by the power term of the sensitivity eigenvalue to obtain the weight update direction. The sign extraction operation maps positive sensitivity to +1, negative sensitivity to -1, and zero sensitivity to 0. In this embodiment, the power term is set to 5 / 3, that is, the absolute value of the sensitivity eigenvalue is raised to the 5 / 3 power and then multiplied by its sign. This non-linear transformation enhances the influence of large sensitivity values while suppressing the influence of small sensitivity values, making the optimization process pay more attention to important regions. For the particles at the boundary of the bleeding area, their sensitivity eigenvalues are usually large, with absolute values between 1.5 and 2.0. After the power transformation, the magnitude of the update direction is about 2.2 to 3.2.

[0095] Calculate the KL divergence of the particle weight distributions in two adjacent iterations to construct an adaptive step size. The KL divergence measures the degree of difference between two probability distributions, and the larger the value, the greater the difference. At the beginning of the iteration, the KL divergence is usually large, with a value between 0.5 and 1.0; as the iteration progresses, the KL divergence gradually decreases and finally converges to less than 0.01. The adaptive step size is associated with the KL divergence through an exponential decay function, and its specific form is the initial step size multiplied by the exponential function of the KL divergence, with the exponential coefficient set to -2. The initial step size is set to 0.2. When the KL divergence is 1.0, the adaptive step size is about 0.027; when the KL divergence is 0.1, the adaptive step size is about 0.164; when the KL divergence is 0.01, the adaptive step size is about 0.196.

[0096] Multiply the adaptive step size by the weight update direction to determine the weight update amount, and add the weight update amount to the particle weight of the current iteration to obtain the updated particle weight. For a typical bleeding detection scenario, the absolute value range of the weight update amount is from 0.001 to 0.5, with an average of about 0.05. At the initial stage of iteration, the update amount is relatively large, and as the iteration progresses, the update amount gradually decreases. When the iteration reaches 30 times, the average absolute value of the update amount is about 0.01, indicating that the optimization process is approaching convergence.

[0097] Calculate the maximum and minimum values of the updated particle weights, and perform renormalization on the updated particle weights to obtain the renormalized particle weights. The renormalization process first truncates the particle weights to ensure that all weight values are within the range of [0, 1], and then normalizes them so that the sum of all particle weights is equal to 1. In this embodiment, the initial particle weights are uniformly set to 1 / 1024 ≈ 0.00098. After optimization, the weights of the particles in the bleeding area usually increase to 0.003 to 0.005, while the weights of the particles in the informationless area decrease to less than 0.0002. The renormalization process ensures the stability of the weight distribution and prevents numerical overflow or underflow.

[0098] Calculate the Euclidean distance between the particle weights of two adjacent iterations to obtain the stability index. The Euclidean distance is calculated as the square root of the sum of the squares of the differences of all particle weights. At the initial stage of iteration, the stability index is usually between 0.1 and 0.2; as the iteration progresses, the stability index gradually decreases. The preset stability threshold is set to 0.001, and the preset number of times is set to 5 times. When the stability index of 5 consecutive iterations is less than 0.001, it is considered that the optimization process has converged, and the particle weights of the current iteration are the optimization results. In a typical traumatic bleeding detection scenario, the optimization usually converges within 40 to 60 iterations.

[0099] Through the above particle weight optimization method, the bleeding areas in infrared thermal imaging and ultrasound images can be accurately identified. The optimized particle weights form an obvious clustering distribution, with high-weight particles concentrated in the bleeding area and low-weight particles distributed in the normal tissue area. At the boundary between the bleeding area and the normal tissue, the particle weights show a gradient distribution, which helps to accurately locate the bleeding boundary. The optimization results are used for subsequent blood volume estimation and bleeding speed calculation, providing an important basis for clinical diagnosis.

[0100] Traditional medical image registration methods are mainly based on gray - level similarity or feature matching, and there are problems such as insufficient accuracy and poor robustness when dealing with infrared thermal imaging and ultrasonic multimodal images. Commonly used particle optimization methods in the prior art, such as particle swarm optimization and evolutionary strategies, are prone to falling into local optimal solutions, and the optimization process has poor stability, making it difficult to be applied to real - time trauma bleeding detection scenarios. The particle weight optimization method proposed in the method of this embodiment introduces a non - linear transformation of the sensitivity eigenvalue, a KL - divergence adaptive step - size adjustment, and an Euclidean - distance stability judgment mechanism, and starting from the convergence and stability of the optimization algorithm, solves the key problems in multimodal medical image registration.

[0101] Compared with the prior art, in this embodiment, the accuracy of the optimization direction is enhanced through the power transformation of the sensitivity eigenvalue, making the weight update pay more attention to important regions; the KL - divergence adaptive step - size mechanism is introduced to balance the convergence speed and stability; the Euclidean - distance stability index is used to judge the convergence condition, improving the reliability of the optimization result.

[0102] As Figure 4 shown, it details the whole process of multimodal medical image registration based on particle weight optimization. The process starts with initializing the particle weights and registration parameters, and then calculates the variational objective function to evaluate the current registration quality. The core optimization link calculates the sensitivity eigenvalue, that is, the partial derivative of the variational objective function with respect to each particle weight, and applies sign extraction and power non - linear transformation to determine the weight update direction. A KL - divergence calculation module is introduced to construct an adaptive step - size mechanism by analyzing the difference in particle weight distributions between two adjacent iterations, thus balancing the convergence speed and stability. The weight update and renormalization steps ensure the effectiveness of the weight update, and keep the weights within a reasonable range by calculating the maximum - minimum value constraint and renormalization processing. The parallel computing module on the right side of the flowchart shows the synchronous processing mechanism of infrared thermal imaging and ultrasonic images, while the parameter monitoring module real - time tracks key parameters such as KL - divergence (0.0023 - 0.0158) and Euclidean distance (0.0067 - 0.0344). The algorithm convergence is judged by the Euclidean - distance stability index. When the continuous preset number is less than the preset stability threshold, the optimization result is output.

[0103] In an alternative embodiment, for the emergency medical rescue scenario, a lightweight configuration is adopted, including a portable infrared thermal imager (weight less than 1 kg) and a handheld ultrasonic probe (weight less than 0.5 kg), and the two are wirelessly connected to the processing terminal. The calibration board is designed as a foldable structure, with a deployed size of 200 mm×200 mm and a thickness of only 10 mm after being folded, which is convenient for field carrying. The whole system can be packed into a medical first - aid kit and carried by a single person.

[0104] It is specially optimized under field rescue conditions. Aiming at the characteristics of large environmental light changes, an environmental temperature compensation algorithm is added to the infrared thermal imaging acquisition process. By setting temperature reference points on the calibration board and dynamically adjusting imaging parameters, the temperature measurement accuracy is ensured to be within ±0.2°C in the environmental temperature range of 0°C to 40°C.

[0105] Aiming at the problem of limited power supply in the wild, the processing algorithm has been optimized in terms of computational complexity. The adaptive shearlet transform adopts a sparse sampling strategy and only calculates key frequency points, reducing the computational complexity by about 70%. The heterogeneous heat diffusion equation is solved using an adaptive grid technique, with the grid encrypted in the boundary region and sparse in the flat region, while maintaining the accuracy, the computational efficiency is increased by 3 times. The structure of the complex-valued neural network has also been streamlined, with the number of hidden layer nodes reduced from 16 to 8, the model volume reduced by 60%, and the single inference time reduced from 120ms to 45ms. These optimizations enable the system to achieve quasi-real-time processing on mid-range mobile processing terminals, with a battery life of up to 4 hours.

[0106] In the trauma bleeding detection process in the rescue scenario, the entire detection process is divided into three stages: rapid screening, precise positioning, and continuous monitoring. In the rapid screening stage, medical staff use an infrared thermal imager to scan the whole body of the wounded, automatically identifying areas with abnormal temperatures, and the whole process takes no more than 30 seconds. In the precise positioning stage, an ultrasonic probe is placed on the area with abnormal temperature, automatically completing multi-modal registration and bleeding feature extraction, and giving the bleeding location, depth, and preliminary quantification results. This stage takes about 90 seconds. In the continuous monitoring stage, the system automatically updates the bleeding dynamic parameters every 5 minutes, including bleeding rate, cumulative blood loss, and tissue damage score, and gives the trend of risk level changes.

[0107] In the emergency medical rescue scenario, compared with traditional inspection methods, the technology of this embodiment does not require the support of large equipment and is applicable to natural disaster sites such as earthquakes and floods and battlefield environments; the detection process is fast and efficient, and it takes no more than 2 minutes from discovering the wounded to completing the preliminary bleeding assessment, saving about 80% of the time compared with traditional methods; the bleeding dynamic parameters provided directly support the judgment of injury grading and transportation priority, which helps to reasonably allocate medical resources in the case of a large number of wounded; the continuous monitoring function enables rescue personnel to timely detect changes in the injury condition and make medical intervention adjustments during transportation.

[0108] This embodiment shows that the technology of the present invention has significant advantages in the emergency medical rescue scenario, can provide reliable trauma bleeding detection under resource-constrained conditions, support rapid medical decision-making, and improve rescue efficiency and treatment success rate. Through lightweight design and algorithm optimization, the practical application problems in field rescue are solved, providing a new technical means for improving emergency medical rescue capabilities.

[0109] In an alternative embodiment, in search and rescue scenarios after natural disasters such as earthquakes and debris flows, rescue workers often face the problem of rapid detection and classification of scattered wounded in a large area. Integrating the method and system onto a lightweight unmanned aerial vehicle (UAV) platform to construct a flexible and mobile aerial rescue detection system, including: a miniaturized infrared thermal imaging module weighing 1.2 kg, with a resolution of 640×480 pixels and a thermal sensitivity of 0.05°C; a portable ultrasonic probe weighing 0.8 kg, with a frequency range of 3 - 12 MHz and a depth penetration of up to 10 cm; an edge computing unit carried, with a processing capacity of 8 TOPS and a power consumption of only 15 W; a wireless transmission module, supporting 5G communication with a transmission delay of less than 100 ms.

[0110] In rescue applications, the UAV first conducts a large-scale scan at an altitude of 200 - 300 meters through infrared thermal imaging technology to identify areas with abnormal surface temperatures. For the identified suspected positions of human heat sources, the UAV reduces its flight altitude to 20 - 30 meters and conducts a higher-resolution infrared scan of the target area to confirm whether it is a trapped person. When it is confirmed as a trapped person, the UAV lands near the wounded, and the rescue workers operate or remotely control the robotic arm to use the ultrasonic probe to examine the possible trauma sites of the wounded. The collected infrared thermal imaging and ultrasonic images are processed by the method of the present invention. The internal bleeding situation at a depth of 2 - 6 cm below the body surface is detected through the deep temperature distribution map, with a sensitivity of 92%; the bleeding position and range are determined through the blood flow velocity field distribution map, with a spatial accuracy better than 3 mm; the severity of the injury is automatically classified based on the bleeding dynamic parameters, and the wounded are divided into three levels: requiring emergency treatment (within 30 minutes), sub-emergency treatment (within 2 hours), and delayable treatment (within 6 hours); the estimated blood loss is calculated in real-time, and when the cumulative blood loss exceeds 800 ml, the system automatically classifies the wounded as critically ill and gives priority to treatment.

[0111] In another alternative embodiment, in a battlefield environment, the triage and rapid treatment of the wounded face multiple challenges such as harsh environments, shortage of personnel, and high pressure. In this embodiment, the method and system are integrated into battlefield first-aid equipment to develop a portable combat trauma detection system; the system is designed with a rugged structure that is waterproof and dustproof (IP67) and shock-resistant (MIL-STD-810G), with a total weight of no more than 3.5 kg, and includes: a high-resolution infrared thermal imager (384×288 pixels), which can work stably at an ambient temperature of -20°C to 50°C; a military-grade ultrasonic probe that supports B-mode, color Doppler, and energy Doppler modes, with a penetration depth of up to 15 cm; an embedded computing unit with a low-power design that can work continuously for 8 hours on a single charge. In the application of combat trauma treatment, rapid assessment of the wounded is achieved. Field medical staff use the infrared thermal imager to scan the whole body of the wounded, which takes about 30 seconds, and automatically identify hot spots with abnormal temperatures; for the identified suspected bleeding areas, the ultrasonic probe is used for local scanning, and the acquisition time is about 45 seconds; algorithms such as spatio-temporal registration, adaptive shear wave transform, thermoacoustic coupling feature extraction, and variational particle filtering in the method of the present invention are executed, and the processing time is less than 15 seconds; a bleeding dynamic parameter report is generated, including: three-dimensional positioning of the bleeding location (accuracy ±5 mm), bleeding rate (error ≤10%), cumulative blood loss (error ≤15%), estimated bleeding duration, and tissue damage degree.

[0112] In the urban pre-hospital emergency system, the rapid assessment and triage decision-making of trauma patients directly affect the treatment effect. In this embodiment, the method and system are integrated into ambulance first-aid equipment to develop a pre-hospital trauma bleeding detection system. The system is designed as a portable device that can be operated by a single person, with a total weight of 2.8 kg, and mainly includes: an integrated host with a touch screen, built-in with a high-performance processor (8 cores, 3.2 GHz); a dual-modal sensor that integrates an infrared thermal imager and an ultrasonic probe in a single handheld device; a lithium battery power supply unit with a continuous working time of up to 4 hours; a wireless transmission module that supports real-time data exchange with the hospital emergency system.

[0113] In the application of pre-hospital first aid, this workflow improves the efficiency of trauma assessment. After the ambulance crew arrives at the scene, they use the system to perform a preliminary scan of the injured, including routine vital sign measurement and infrared thermal imaging of the trauma site; based on the results of the infrared thermal imaging, the system automatically guides the ambulance crew to perform ultrasonic scans of the suspicious areas; after the dual-modal data collected is processed by the method of the present invention, a bleeding assessment report is generated within 60 seconds; the system automatically calculates the Modified Trauma Score based on the bleeding dynamic parameters and gives treatment suggestions.

[0114] The method and system for detecting trauma bleeding based on infrared thermal imaging and ultrasound multimodality in the embodiments of the present invention include: The first unit is used to collect infrared thermal imaging images and ultrasonic images of the trauma site, perform spatio-temporal registration on the images through a calibration plate, and obtain registered image data; The second unit is used to perform adaptive shear wave transform on the registered image data to obtain temperature field coefficients and acoustic coefficients, substitute the temperature field coefficients into the heterogeneous heat diffusion equation to calculate the deep temperature distribution map and obtain the anisotropic heat diffusion coefficient, and perform pulse compression and phase correction on the acoustic coefficients to obtain the tissue stratification map; The third unit is used to construct a thermoacoustic coupling feature matrix based on the deep temperature distribution map and the tissue stratification map by using a complex-valued neural network, extract the blood flow velocity vector from the Doppler frequency shift signal based on the thermoacoustic coupling feature matrix, and reconstruct the vascular structure by using Marching Cubes interpolation to generate a blood flow velocity field distribution map; The fourth unit is used to calculate the local blood flow flux according to the blood flow velocity field distribution map, perform volume integration in combination with the anisotropic heat diffusion coefficient to obtain the cumulative blood loss, and establish a time series mapping relationship between the local blood flow flux and the cumulative blood loss; The fifth unit is used to dynamically track the local blood flow flux in the time series mapping relationship by using the variational particle filter method, calculate the corresponding time-varying coefficient, and construct a state equation to obtain the bleeding dynamic parameters; The sixth unit is used to calculate the bleeding duration and the degree of tissue damage according to the bleeding dynamic parameters, and determine the risk level and treatment plan.

[0115] In the third aspect of the embodiments of the present invention, an electronic device is provided, including: A processor; A memory for storing instructions executable by the processor; Wherein, the processor is configured to call the instructions stored in the memory to execute the method described above.

[0116] In the fourth aspect of the embodiments of the present invention, a computer-readable storage medium is provided, on which computer program instructions are stored, and when the computer program instructions are executed by a processor, the method described above is implemented.

[0117] The present invention can be a method, a device, a system, and / or a computer program product. The computer program product can include a computer-readable storage medium, on which computer-readable program instructions for executing various aspects of the present invention are uploaded.

[0118] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that they can still modify the technical solutions described in the foregoing embodiments, or perform equivalent replacements on some or all of the technical features; and these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A trauma bleeding detection method based on infrared thermal imaging and ultrasound multimodality, characterized in that Including: Collect infrared thermal imaging images and ultrasonic images of the trauma site, perform spatio-temporal registration on the images through a calibration plate to obtain registered image data; Perform adaptive shearlet transform on the registered image data to obtain temperature field coefficients and acoustic coefficients, substitute the temperature field coefficients into the heterogeneous heat diffusion equation to calculate the deep temperature distribution map and obtain the anisotropic heat diffusion coefficient, and perform pulse compression and phase correction on the acoustic coefficients to obtain the tissue stratification map; According to the deep temperature distribution map and the tissue stratification map, construct a thermoacoustic coupling feature matrix using a complex-valued neural network, extract the blood flow velocity vector from the Doppler shift signal based on the thermoacoustic coupling feature matrix, and reconstruct the vascular structure using Marching Cubes interpolation to generate a blood flow velocity field distribution map; Calculate the local blood flow flux according to the blood flow velocity field distribution map, perform volume integration in combination with the anisotropic heat diffusion coefficient to obtain the cumulative blood loss, and establish a temporal mapping relationship between the local blood flow flux and the cumulative blood loss; Use the variational particle filter method to dynamically track the local blood flow flux in the temporal mapping relationship, calculate the corresponding time-varying coefficient, and construct a state equation to obtain the bleeding dynamic parameters; Calculate the bleeding duration and tissue damage degree according to the bleeding dynamic parameters, and determine the risk level and treatment plan.

2. The method according to claim 1, wherein Performing an adaptive shearlet transform on the registered image data to obtain temperature field coefficients and acoustic coefficients, substituting the temperature field coefficients into the heterogeneous heat diffusion equation to calculate the deep temperature distribution map and obtain the anisotropic heat diffusion coefficient, and performing pulse compression and phase correction on the acoustic coefficients to obtain the tissue stratification map includes: Perform an adaptive shearlet transform on the registered image data, wherein the adaptive shearlet transform performs joint decomposition on the registered image data by setting a direction-selective filter bank, constructs a shearlet basis function according to the shear parameter and the translation parameter, and obtains a coefficient matrix through two-dimensional integral operation on the shearlet basis function, calculates the direction angle parameter based on the arctangent function of the shear parameter, and optimizes the direction angle parameter to achieve multi-directional adaptive decomposition to obtain temperature field coefficients and acoustic coefficients; Substitute the temperature field coefficients into the heterogeneous heat diffusion equation, wherein the heterogeneous heat diffusion equation includes a tissue density term, a specific heat capacity term, a spatially varying thermal conductivity tensor term, and a heat source term, calculate the product of the tissue density term and the specific heat capacity term to obtain a coefficient term, calculate the inner product of the thermal conductivity tensor term and the temperature gradient to obtain a diffusion term, substitute the coefficient term, the diffusion term, and the heat source term into the heterogeneous heat diffusion equation, and use the finite difference method to solve the heterogeneous heat diffusion equation to obtain the deep temperature distribution map, and calculate the anisotropic heat diffusion coefficient based on the ratio of the maximum thermal conductivity to the minimum thermal conductivity; Perform pulse compression processing on the acoustic coefficients to obtain the compressed echo signal, perform phase correction by calculating the difference between the original phase and the reference phase, calculate the acoustic reflection coefficient according to the ratio of the difference between the acoustic impedances of adjacent tissues to the sum value, and reconstruct the tissue stratification map according to the acoustic reflection coefficient, the compressed echo signal, and the phase correction result.

3. The method according to claim 2, wherein Using the finite difference method to solve the heterogeneous heat diffusion equation to obtain the deep temperature distribution map includes: Discretize the time and space coordinates, divide the time dimension into multiple time layers according to the time step, divide the space dimension into a three-dimensional grid structure according to the space step, construct a differential grid, and obtain the initial temperature value at the grid point of the differential grid; performing differential discretization processing on the heterogeneous heat diffusion equation, discretizing the time derivative term using a first-order forward difference format, discretizing the space derivative term using a second-order central difference format, interpolating the thermal conductivity tensor at half-grid points to obtain a discrete thermal conductivity value, and substituting the discrete time derivative term, the discrete space derivative term, and the discrete thermal conductivity value into the heterogeneous heat diffusion equation to establish a differential equation; Constructing boundary conditions at the boundary positions of the differential grid, discretizing the boundary conditions of the temperature boundary, the heat flow boundary, and the convection boundary, and combining the boundary conditions with the differential equations to construct a complete differential equation group; The complete differential equations are decomposed and solved in the x-, y-, and z-directions in sequence using an alternating direction implicit method, a tridiagonal equation system is constructed for each direction, and an iterative solution is performed using a chasing method to obtain the temperature distribution of the new time layer; When the temperature difference between adjacent time layers is less than a preset convergence threshold, a deep temperature distribution map is generated based on the final converged temperature distribution.

4. The method according to claim 1, wherein According to the deep temperature distribution map and tissue layer map, a complex-valued neural network is used to construct a thermoacoustic coupling characteristic matrix. Based on the thermoacoustic coupling characteristic matrix, the blood flow velocity vector is extracted from the Doppler frequency shift signal. The vascular structure is reconstructed using Marchi cubic interpolation to generate a blood flow velocity field distribution map including: Normalizing the deep temperature distribution map and applying Gaussian filtering to obtain a temperature field gradient, performing adaptive threshold segmentation on the tissue layer map to obtain a tissue interface normal vector, and combining the temperature field gradient and the tissue interface normal vector to form an initial feature vector; Inputting the initial eigenvector into a complex-valued neural network, wherein the complex-valued neural network adopts a three-layer network structure and has a piecewise linear unit function as an activation function, and processing the real and imaginary parts of the initial eigenvector by the complex-valued neural network to obtain a thermoacoustic coupling characteristic matrix; Obtaining a complex analytical signal by Hilbert transforming the Doppler frequency shift signal, and performing a tensor product operation on the complex analytical signal and the thermoacoustic coupling characteristic matrix to obtain a three-dimensional blood flow velocity vector field; Calculating velocity divergence distribution according to the three-dimensional blood flow velocity vector field, determining a blood vessel boundary point set, and inputting the blood vessel boundary point set into a cubic Marchi cube interpolation algorithm to obtain a blood vessel surface; Within the blood vessel surface, radial basis functions are used to perform interpolation calculation on the three-dimensional blood flow velocity vector field to obtain a continuous blood flow velocity field distribution map.

5. The method according to claim 4, wherein Calculating velocity divergence distribution according to the three-dimensional blood flow velocity vector field, determining a blood vessel boundary point set, and inputting the blood vessel boundary point set into a cubic Marchi cube interpolation algorithm to obtain a blood vessel surface includes: The scale pyramid decomposition is performed on the three-dimensional blood flow velocity vector field to obtain velocity vector fields at multiple scale levels. The velocity divergence distributions are calculated for the velocity vector fields respectively to determine the multi-scale divergence distribution. Based on the multi-scale divergence distribution, feature enhancement is performed through average pooling and a multi-layer perceptron to obtain the enhanced velocity divergence distribution; The enhanced velocity divergence distribution is constructed into a graph structure including a vertex set, an edge set, and an adjacency matrix. The graph structure is input into a graph attention network for feature iterative update to obtain the blood vessel boundary features; Based on the blood vessel boundary features, a set of skeleton blood vessel boundary points is extracted and recursively refined to obtain a complete set of blood vessel boundary points; A joint loss function including a topological consistency constraint term and a geometric constraint term is constructed. Based on the joint loss function, the complete set of blood vessel boundary points is optimized to determine the optimized set of blood vessel boundary points; The optimized set of blood vessel boundary points is divided into local surface patches and a parameterized grid is established. The control vertices of the parameterized grid are calculated. Based on the control vertices, a cubic B-spline function is constructed and a tensor product operation is performed with the position information of the control vertices to obtain the local surface reconstruction result. After applying the curvature continuity constraint and the boundary continuity constraint to the local surface reconstruction result, they are stitched together to obtain the blood vessel surface.

6. The method according to claim 1, characterized in that The variational particle filter method is used to dynamically track the local blood flux in the temporal mapping relationship, calculate the corresponding time-varying coefficients, and construct a state equation to obtain the bleeding dynamic parameters including: According to the temporal mapping relationship, multiple particles are initialized and initial weights are assigned. Each particle corresponds to a state of the local blood flux; An observation likelihood function is established to describe the relationship between the cumulative blood loss and the local blood flux, and a state transition probability function is established to describe the change law of the local blood flux; The observation likelihood function and the state transition probability function are combined to construct a variational objective function; Based on the variational objective function, the weights of the particles are updated by iterative optimization; According to the updated particle weights, the optimal estimated value of the local blood flux is calculated; Importance resampling is performed to eliminate particle degeneracy to obtain the dynamic change result of the local blood flux; According to the dynamic change result of the local blood flux, the change law of the blood flux in different time periods is analyzed; Based on the change law, the scale coefficient, delay coefficient, and attenuation coefficient of the hemodynamic characteristics are calculated and temporal smoothing processing is performed to determine each smoothing coefficient; The smoothing coefficients are used to construct the state transition matrix of the nonlinear dynamics model, and the system noise term is combined with the state transition matrix to construct the bleeding dynamic characteristic state equation; The bleeding dynamic characteristic state equation is solved to calculate the bleeding rate, cumulative blood loss, and tissue damage degree to determine the bleeding dynamic parameters.

7. The method according to claim 6, characterized in that, Updating the weights of the particles by iterative optimization based on the variational objective function includes: Calculating the partial derivative of the variational objective function with respect to the weights of each particle to obtain the sensitivity eigenvalue; Extracting the sign of the sensitivity eigenvalue and multiplying it by the power term of the sensitivity eigenvalue to obtain the weight update direction; Calculate the KL divergence of the particle weight distributions in two adjacent iterations, construct an adaptive step size, multiply the adaptive step size by the weight update direction to determine the weight update amount, and add the weight update amount to the particle weights in the current iteration to obtain the updated particle weights; Calculate the maximum and minimum values of the updated particle weights, and perform renormalization processing on the updated particle weights to obtain the renormalized particle weights; Calculate the Euclidean distance between the particle weights in two adjacent iterations to obtain a stability index. When the stability index is less than a preset stability threshold for a preset number of consecutive times, determine the particle weights in the current iteration as the optimization result.

8. A trauma bleeding detection method and system based on infrared thermal imaging and ultrasound multimodality, for implementing the method described in any one of the preceding claims 1-7, characterized in that, Comprising: A first unit for collecting infrared thermal imaging images and ultrasonic images of a trauma site, performing spatio-temporal registration on the images through a calibration plate to obtain registered image data; A second unit for performing an adaptive shearlet transform on the registered image data to obtain temperature field coefficients and acoustic coefficients, substituting the temperature field coefficients into a heterogeneous heat diffusion equation to calculate a deep temperature distribution map and obtain an anisotropic heat diffusion coefficient, and performing pulse compression and phase correction on the acoustic coefficients to obtain a tissue stratification map; A third unit for constructing a thermoacoustic coupling feature matrix using a complex-valued neural network based on the deep temperature distribution map and the tissue stratification map, extracting a blood flow velocity vector from a Doppler shift signal based on the thermoacoustic coupling feature matrix, and reconstructing a blood vessel structure using Marching Cubes interpolation to generate a blood flow velocity field distribution map; A fourth unit for calculating a local blood flow rate according to the blood flow velocity field distribution map, performing a volume integral in combination with the anisotropic heat diffusion coefficient to obtain an accumulated blood loss amount, and establishing a temporal mapping relationship between the local blood flow rate and the accumulated blood loss amount; A fifth unit for dynamically tracking the local blood flow rate in the temporal mapping relationship using a variational particle filter method, calculating a corresponding time-varying coefficient, and constructing a state equation to obtain bleeding dynamic parameters; A sixth unit for calculating a bleeding duration and a tissue damage degree according to the bleeding dynamic parameters, and determining a risk level and a treatment plan.

9. An electronic device, characterized in that, Comprising: A processor; A memory for storing instructions executable by the processor; Wherein, the processor is configured to call the instructions stored in the memory to execute the method according to any one of claims 1 to 7.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, The computer program instructions, when executed by the processor, implement the method according to any one of claims 1 to 7.

Images