Traumatic bleeding detection method and system based on infrared thermal imaging and ultrasound multimodality

Abstract

This invention provides a method and system for detecting traumatic bleeding based on infrared thermal imaging and ultrasound multimodality. This method relates to the field of medical image processing technology. The method involves acquiring and registering infrared thermal imaging and ultrasound images of the wound site, performing adaptive shear wave transform to extract characteristic coefficients, constructing a thermoacoustic coupling characteristic matrix, reconstructing vascular structure to generate a blood flow velocity field distribution map, establishing a temporal mapping relationship between local blood flow and cumulative bleeding volume, and using a variational particle filter to calculate dynamic bleeding parameters. This method can accurately assess the duration of traumatic bleeding and the extent of tissue damage, providing a precise diagnostic basis for clinical emergency treatment.

Description

Technical Field

[0001] The present invention relates to the technical field of medical image processing, and in particular to a trauma bleeding detection method and system based on infrared thermal imaging and ultrasonic multimodality. Background Art

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

[0003] However, existing trauma bleeding detection technologies still have several defects and shortcomings. The detection methods of a single imaging modality have limited information and cannot fully reflect the complex state of bleeding. Infrared thermal imaging can only reflect the surface temperature distribution and cannot accurately assess the bleeding situation in deep tissues. Although simple ultrasound examination can display tissue structure, it lacks sensitivity in detecting tiny blood vessels and low-speed blood flow, making it easy to miss diagnoses. Existing multimodal image fusion methods lack an effective spatiotemporal registration mechanism, making it difficult to achieve accurate correspondence between infrared thermal imaging and ultrasound images, resulting in distortion of fused information and deviations in clinical judgment. Especially in dynamic scenarios, the difference in acquisition time and spatial resolution between the two modalities makes traditional registration methods ineffective. Existing bleeding 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 changes in blood flow rate and cumulative effects, nor can they predict the duration of bleeding and the degree of tissue damage, thus affecting the timeliness and accuracy of clinical treatment decisions. Summary of the Invention

[0004] The embodiments of the present invention provide a trauma bleeding detection method and system based on infrared thermal imaging and ultrasonic multimodality, which can solve the problems in the prior art.

[0005] A first aspect of an embodiment of the present invention provides a method for detecting traumatic bleeding based on infrared thermal imaging and ultrasound multimodality, comprising:

[0006] Collect infrared thermal imaging images and ultrasound images of the wound site, perform spatiotemporal registration of the images using a calibration plate, and obtain registered image data;

[0007] Adaptive shear wave transform is performed on the registered image data to obtain temperature field coefficients and acoustic coefficients. The temperature field coefficients are substituted into the heterogeneous heat diffusion equation to calculate the deep temperature distribution map and obtain the anisotropic thermal diffusion coefficient. Pulse compression and phase correction are performed on the acoustic coefficients to obtain the tissue layer map.

[0008] Based on 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 cube interpolation to generate a blood flow velocity field distribution map.

[0009] The local blood flow flux is calculated based on the blood flow velocity field distribution map, and the cumulative bleeding volume is obtained by volume integration combined with the anisotropic thermal diffusion coefficient. The time series mapping relationship between the local blood flow flux and the cumulative bleeding volume is established.

[0010] 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, and construct the state equation to obtain the dynamic bleeding parameters.

[0011] The duration of bleeding and the degree of tissue damage are calculated based on the dynamic parameters of bleeding to determine the risk level and treatment plan.

[0012] In an optional embodiment, performing adaptive shear wave 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 anisotropic thermal diffusion coefficients, and performing pulse compression and phase correction on the acoustic coefficients to obtain a tissue layer map includes:

[0013] Performing an adaptive shearlet transform on the registered image data, wherein the adaptive shearlet transform jointly decomposes the registered image data by setting a direction-selective filter group, constructs a shearlet basis function according to a shearing parameter and a translation parameter, obtains a coefficient matrix by performing a two-dimensional integral operation on the shearlet basis function, calculates an angular direction parameter based on an inverse tangent function of the shearing parameter, and optimizes the angular direction parameter to achieve multi-directional adaptive decomposition and obtain a temperature field coefficient and an acoustic coefficient;

[0014] Substituting the temperature field coefficient into a 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; calculating the coefficient term by multiplying the tissue density term by the specific heat capacity term; calculating the diffusion term by the inner product of the thermal conductivity tensor term and the temperature gradient; substituting the coefficient term, the diffusion term, and the heat source term into the heterogeneous heat diffusion equation; solving the heterogeneous heat diffusion equation using a finite difference method to obtain a deep temperature distribution map; and calculating the anisotropic thermal diffusion coefficient based on the ratio of maximum thermal conductivity to minimum thermal conductivity;

[0015] 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.

[0016] In an optional embodiment, solving the heterogeneous heat diffusion equation using a finite difference method to obtain a deep temperature distribution map includes:

[0017] 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;

[0018] 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, 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;

[0019] 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;

[0020] 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;

[0021] 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.

[0022] 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:

[0023] 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;

[0024] 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;

[0025] 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;

[0026] 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;

[0027] 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.

[0028] In an optional embodiment, calculating the 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 the blood vessel surface includes:

[0029] The three-dimensional blood flow velocity vector field is decomposed into a scale pyramid to obtain velocity vector fields at multiple scale levels. The velocity divergence distribution of the velocity vector field is calculated separately to determine the multi-scale divergence distribution. Based on the multi-scale divergence distribution, feature enhancement is performed through average pooling and multi-layer perceptron to obtain the enhanced velocity divergence distribution.

[0030] The enhanced velocity divergence distribution is constructed as a graph structure including a vertex set, an edge set, and an adjacency matrix, and the graph structure is input into a graph attention network to perform feature iterative update to obtain blood vessel boundary features;

[0031] Extracting a skeleton blood vessel boundary point set based on the blood vessel boundary features, and recursively refining it to obtain a complete blood vessel boundary point set;

[0032] constructing a joint loss function including a topological consistency constraint term and a geometric constraint term, optimizing the complete blood vessel boundary point set based on the joint loss function, and determining an optimized blood vessel boundary point set;

[0033] The optimized vascular boundary point set is divided into local surface patches and a parameterized grid is established. The control vertices of the parameterized grid are calculated. A cubic Marchi basis function is constructed based on the control vertices. A tensor product operation is performed with the position information of the control vertices to obtain a local surface reconstruction result. After applying curvature continuity constraints and boundary continuity constraints to the local surface reconstruction result, the vascular surface is spliced.

[0034] In an optional embodiment, a variational particle filter method is used to dynamically track the local blood flow in the time series mapping relationship, calculate the corresponding time-varying coefficient, and construct a state equation to obtain the dynamic bleeding parameters, including:

[0035] Initializing multiple particles according to the time series mapping relationship and assigning initial weights, each particle corresponding to a state of local blood flow;

[0036] Establishing an observation likelihood function to describe the relationship between the cumulative bleeding volume and the local blood flow, and establishing a state transition probability function to describe the change pattern of the local blood flow; combining the observation likelihood function and the state transition probability function to construct a variational objective function;

[0037] updating the particle weights through iterative optimization based on the variational objective function; calculating the optimal estimate of the local blood flow according to the updated particle weights; performing importance resampling to eliminate particle degradation and obtain a dynamic change result of the local blood flow;

[0038] Analyzing the blood flow variation patterns in different time periods based on the dynamic variation results of the local blood flow; calculating the scale coefficient, delay coefficient, and attenuation coefficient of the hemodynamic characteristics based on the variation patterns, and performing time series smoothing to determine the smoothing coefficients;

[0039] A state transfer matrix of a nonlinear dynamic model is constructed using a smoothing coefficient, and the system noise term is combined with the state transfer matrix to construct a state equation of the dynamic characteristics of bleeding. The state equation of the dynamic characteristics of bleeding is solved to calculate the bleeding rate, cumulative bleeding volume and degree of tissue damage, and determine the dynamic parameters of bleeding.

[0040] In an optional embodiment, updating the weights of particles through iterative optimization based on the variational objective function includes:

[0041] Calculating the partial derivative of the variational objective function with respect to each particle weight to obtain a sensitivity eigenvalue; performing symbol extraction on the sensitivity eigenvalue and multiplying it with a power term of the sensitivity eigenvalue to obtain a weight update direction;

[0042] Calculate the KL divergence of the particle weight distribution of 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;

[0043] Calculating the maximum and minimum values ​​of the updated particle weights, and performing renormalization processing on the updated particle weights to obtain renormalized particle weights;

[0044] The Euclidean distance between the particle weights of two adjacent iterations is calculated to obtain a stability index. When the stability index is less than a preset stability threshold for a consecutive preset number of times, the particle weight of the current iteration is determined to be the optimization result.

[0045] A second aspect of an embodiment of the present invention provides a trauma bleeding detection method system based on infrared thermal imaging and ultrasound multimodality, comprising:

[0046] The first unit is used to collect infrared thermal imaging images and ultrasound images of the wound site, and perform spatiotemporal registration of the images through a calibration plate to obtain registered image data;

[0047] 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 thermal diffusion coefficient, and perform pulse compression and phase correction on the acoustic coefficients to obtain the tissue layer map;

[0048] The third unit is used to construct a thermoacoustic coupling characteristic matrix based on the deep temperature distribution map and tissue layer map using a complex-valued neural network, extract the blood flow velocity vector from the Doppler frequency shift signal based on the thermoacoustic coupling characteristic matrix, reconstruct the vascular structure using Marchi cube interpolation, and generate a blood flow velocity field distribution map;

[0049] The fourth unit is used to calculate the local blood flow flux based on the blood flow velocity field distribution map, perform volume integration based on the anisotropic thermal diffusion coefficient to obtain the cumulative bleeding volume, and establish a time series mapping relationship between the local blood flow flux and the cumulative bleeding volume;

[0050] The fifth unit is used to dynamically track the local blood flow in the time series mapping relationship using the variational particle filter method, calculate the corresponding time-varying coefficients, and construct the state equation to obtain the dynamic parameters of bleeding;

[0051] Unit 6 is used to calculate the duration of bleeding and the degree of tissue damage based on the dynamic parameters of bleeding, and to determine the risk level and treatment plan.

[0052] According to a third aspect of an embodiment of the present invention, an electronic device is provided, including:

[0053] processor;

[0054] a memory for storing processor-executable instructions;

[0055] The processor is configured to call the instructions stored in the memory to execute the aforementioned method.

[0056] According to a fourth aspect of an embodiment of the present invention, a computer-readable storage medium is provided, on which computer program instructions are stored. When the computer program instructions are executed by a processor, the method described above is implemented.

[0057] In an embodiment of the present invention, through the multimodal fusion of infrared thermal imaging and ultrasound technology, accurate detection of traumatic bleeding is achieved, which solves the technical problem that traditional methods are difficult to accurately evaluate deep tissue bleeding and improves the diagnostic accuracy of bleeding location and degree; adaptive shear wave transform and complex-valued neural network are used to process and align image data, and a thermoacoustic coupling feature matrix is ​​constructed, which can deeply analyze tissue structure and blood flow dynamic characteristics, and reconstruct vascular structure through Marchi cube interpolation to achieve accurate calculation of bleeding volume, greatly improving the sensitivity and specificity of traumatic bleeding detection; the variational particle filter method is introduced to dynamically track and analyze blood flow volume, and a correlation model between bleeding dynamic parameters and tissue damage degree is established, which can predict the development trend of bleeding in real time, provide a scientific basis for clinical formulation of individualized emergency strategies, and significantly improve the timeliness and success rate of trauma treatment. BRIEF DESCRIPTION OF THE DRAWINGS

[0058] Figure 1 Schematic diagram of the process of a trauma bleeding detection method based on infrared thermal imaging and ultrasound multimodality according to an embodiment of the present invention;

[0059] Figure 2 Comparison chart of computational efficiency for solving heterogeneous heat diffusion equations;

[0060] Figure 3 Comparison of curvature continuity of vascular surface reconstruction using different methods;

[0061] Figure 4 Flowchart of particle weight optimization for multimodal medical images. DETAILED DESCRIPTION

[0062] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts shall fall within the scope of protection of the present invention.

[0063] The following specific embodiments are used to describe the technical solution of the present invention in detail. The following specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described in detail in some embodiments.

[0064] Figure 1 FIG. 1 is a flow chart of a method for detecting traumatic bleeding based on infrared thermal imaging and ultrasound multimodality according to an embodiment of the present invention. Figure 1 As shown, the method includes:

[0065] Collect infrared thermal imaging images and ultrasound images of the wound site, perform spatiotemporal registration of the images using a calibration plate, and obtain registered image data;

[0066] Adaptive shear wave transform is performed on the registered image data to obtain temperature field coefficients and acoustic coefficients. The temperature field coefficients are substituted into the heterogeneous heat diffusion equation to calculate the deep temperature distribution map and obtain the anisotropic thermal diffusion coefficient. Pulse compression and phase correction are performed on the acoustic coefficients to obtain the tissue layer map.

[0067] Based on 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 cube interpolation to generate a blood flow velocity field distribution map.

[0068] The local blood flow flux is calculated based on the blood flow velocity field distribution map, and the cumulative bleeding volume is obtained by volume integration combined with the anisotropic thermal diffusion coefficient. The time series mapping relationship between the local blood flow flux and the cumulative bleeding volume is established.

[0069] 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, and construct the state equation to obtain the dynamic bleeding parameters.

[0070] The duration of bleeding and the degree of tissue damage are calculated based on the dynamic parameters of bleeding to determine the risk level and treatment plan.

[0071] In an optional embodiment, performing adaptive shear wave 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 thermal diffusion coefficient, and performing pulse compression and phase correction on the acoustic coefficients to obtain the tissue layer map includes:

[0072] Performing an adaptive shearlet transform on the registered image data, wherein the adaptive shearlet transform jointly decomposes the registered image data by setting a direction-selective filter group, constructs a shearlet basis function according to a shearing parameter and a translation parameter, obtains a coefficient matrix by performing a two-dimensional integral operation on the shearlet basis function, calculates an angular direction parameter based on an inverse tangent function of the shearing parameter, and optimizes the angular direction parameter to achieve multi-directional adaptive decomposition and obtain a temperature field coefficient and an acoustic coefficient;

[0073] Substituting the temperature field coefficient into a 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; calculating the coefficient term by multiplying the tissue density term by the specific heat capacity term; calculating the diffusion term by the inner product of the thermal conductivity tensor term and the temperature gradient; substituting the coefficient term, the diffusion term, and the heat source term into the heterogeneous heat diffusion equation; solving the heterogeneous heat diffusion equation using a finite difference method to obtain a deep temperature distribution map; and calculating the anisotropic thermal diffusion coefficient based on the ratio of maximum thermal conductivity to minimum thermal conductivity;

[0074] 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.

[0075] In one embodiment, the process of performing an adaptive shearlet transform on the registered image data includes setting a directionally selective filter bank comprising eight directional bandpass filters 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, and a bandwidth of ±15% of the center frequency. Joint decomposition of the registered image data is performed using 256×256 pixel image blocks as processing units, with each image block convolved with all filters to obtain an initial decomposition result.

[0076] When constructing the shearlet basis function based on the shearing and translation parameters, the shearing parameter ranges from [-2.0, 2.0] with a step size of 0.1. The translation parameter is uniformly sampled in the image space, with 32 sampling points in the x and y directions. For each parameter combination, a corresponding shearlet basis function is constructed. Specifically, the shearlet basis function takes the form of a Gabor function, which contains real and imaginary parts. The real part function is expressed as a Gaussian modulated cosine function, and the imaginary part function is expressed as a Gaussian modulated sine function. The standard deviation of the Gaussian envelope function is set to 1 / 8 the image block size.

[0077] When performing a two-dimensional integration of the shearlet basis functions to obtain the coefficient matrix, we approximate the integration process using a discrete summation approach. For each 256×256 image block, we calculate the dot product between it and the shearlet basis function to obtain the corresponding coefficient value. For all parameter combinations, we form a coefficient matrix with the dimensions (number of shearing parameters × number of translation parameters × number of directions).

[0078] When calculating the azimuth angle parameters based on the inverse tangent function of the shear parameter, for a shear parameter s, the azimuth angle θ is calculated as θ = arctg(s). When s = -2.0, θ is approximately -63.4 degrees; when s = 0, θ is 0 degrees; and when s = 2.0, θ is approximately 63.4 degrees. To optimize the azimuth angle parameters, an adaptive threshold mechanism is introduced. Based on the energy distribution of the coefficient matrix, azimuth angles whose energy contribution exceeds 5% of the total energy are retained, achieving multi-directional adaptive decomposition.

[0079] The temperature field coefficient and acoustic coefficient are obtained through processing. The temperature field coefficient corresponds to the infrared thermal imaging modality, while the acoustic coefficient corresponds to the ultrasound modality. In a typical traumatic hemorrhage case, the temperature field coefficient manifests as a higher low-frequency component and energy concentration in a specific direction in the hemorrhage area, while the acoustic coefficient manifests as echo enhancement and tissue structural changes in the hemorrhage area.

[0080] Substituting the temperature field coefficient into the heterogeneous heat diffusion equation for solution, the tissue density ranged from 900 to 1100 kg / m³, with the density of the hemorrhagic region approximately 1060 kg / m³ and that of normal tissue approximately 1020 kg / m³. The specific heat capacity term ranged from 3500 to 4200 J / (kg·°C), with the specific heat capacity of the hemorrhagic region approximately 3900 J / (kg·°C) and that of normal tissue approximately 3700 J / (kg·°C).

[0081] The coefficient term is calculated by multiplying the tissue density term by the specific heat capacity term. For the hemorrhage region, the coefficient term is 1060 × 3900 = 4,134,000 J / (m³·°C). The thermal conductivity tensor term is represented by a 2×2 matrix, with the thermal conductivity values ​​along the principal axis ranging from 0.3 to 0.6 W / (m·°C). The principal axis thermal conductivity of the hemorrhage region is approximately [0.53, 0.42] W / (m·°C), while the thermal conductivity along the non-principal axis is approximately 0.1 W / (m·°C).

[0082] The diffusion term is calculated by calculating the inner product of the thermal conductivity tensor and the temperature gradient. The temperature gradient is calculated using the central difference method with a spatial step size of 1 mm. At the boundary of the hemorrhage area, the typical temperature gradient is approximately 0.5-1.5°C / mm. The heat source term is set based on physiological properties. The heat source value for normal tissue is approximately 420 W / m³, while the heat source in the hemorrhage area is enhanced due to blood extravasation, reaching a value of approximately 650 W / m³.

[0083] The heterogeneous heat diffusion equation was solved using the finite difference method with a time step of 0.1 seconds and a spatial step of 1 mm, using an explicit iterative scheme. Boundary conditions were set to third-order boundary conditions, with a surface heat transfer coefficient of 12 W / (m²·°C) and an ambient temperature of 25°C. After 100 iterations, the temperature field reached a steady state, and a deep temperature distribution map was obtained.

[0084] The anisotropic thermal diffusivity is calculated based on the ratio of maximum to minimum thermal conductivity. For example, in the hemorrhage area, the anisotropic thermal diffusivity is 0.53 ÷ 0.42 = 1.26, indicating that heat conduction in the hemorrhage area exhibits directional differences. The anisotropic thermal diffusivity of normal tissue is approximately 1.05-1.15, while that of the hemorrhage area is approximately 1.20-1.35.

[0085] Pulse compression is performed on the acoustic coefficients to produce a compressed echo signal. Using a 13-bit Barker code as the modulation sequence, the compression gain is approximately 10.8 dB. Compression utilizes matched filtering, with the filter coefficients being the time-reversed conjugate form of the modulation sequence. For clinical data, the axial resolution after compression is improved to 0.3 mm.

[0086] Phase correction is performed by calculating the difference between the original phase and the reference phase. The reference phase is obtained by statistically averaging healthy tissue areas. For hemorrhagic areas, the phase difference is typically 15-40 degrees, and the corrected phase error is controlled within ±5 degrees.

[0087] The acoustic reflection coefficient is calculated based on the ratio of the difference and 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 6 The 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.

[0088] 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.

[0089] 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.

[0090] In an optional embodiment, solving the heterogeneous heat diffusion equation using a finite difference method to obtain a deep temperature distribution map includes:

[0091] 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;

[0092] 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, 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;

[0093] 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;

[0094] 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;

[0095] 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.

[0096] In a specific embodiment, when the space-time coordinates are discretized, the time dimension is discretized using a fixed time step. In this embodiment, the time step is set to 0.1 seconds, the total simulation time is 10 seconds, and a total of 100 time layers are planned. The spatial dimension is discretized using a structured grid method, and the spatial steps 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 trauma bleeding detection scenario, the constructed three-dimensional differential grid covers an area of ​​100 mm × 100 mm × 60 mm, and the corresponding number of grid points is 201 × 201 × 121.

[0097] When obtaining initial temperature values ​​at the grid points of the differential grid, for skin surface grid points, temperature values ​​were directly obtained using an infrared thermal imager, with a typical temperature range of 32°C to 36°C. For subcutaneous tissue grid points, the initial temperature was estimated using an empirical formula: approximately 36.5°C at 1 cm below the surface, approximately 37.0°C at 2 cm, and approximately 37.2°C at 3 cm and deeper. For hemorrhagic areas, the initial temperature was set to 37.3°C to 37.5°C, slightly higher than the temperature of surrounding normal tissue.

[0098] When performing differential discretization on the heterogeneous heat diffusion equation, the time derivative term is discretized using a first-order forward difference scheme. Specifically, the time derivative of temperature is approximated as the temperature difference between the current and next time layers divided by the time step. In this embodiment, for grid point (i, j, k), the time derivative term is expressed as the temperature of the next time layer minus the temperature of the current time layer, divided by 0.1 seconds.

[0099] Spatial derivatives are discretized using a second-order central difference scheme. Second-order derivatives are calculated using the temperature values ​​of three adjacent grid points. In the x-direction, the second-order derivative is expressed as the temperature value at grid point (i+1,j,k) plus the temperature value at grid point (i-1,j,k), minus twice the temperature value at grid point (i,j,k), divided by the square of the x-step size. Second-order derivatives in the y- and z-directions are discretized similarly. Mixed derivatives are calculated using the temperature values ​​of four adjacent diagonal grid points.

[0100] The thermal conductivity tensor is interpolated at the half-grid points to obtain discrete values ​​of thermal conductivity. The thermal conductivity tensor contains 9 components, corresponding to the three main directions and their cross terms. In this embodiment, the main diagonal elements of the thermal conductivity tensor of the hemorrhage area are 0.52 W / (m·℃), 0.48 W / (m·℃), and 0.45 W / (m·℃), and the off-diagonal elements are 0.08 W / (m·℃), 0.06 W / (m·℃), and 0.07 W / (m·℃). The main diagonal elements of the thermal conductivity tensor of normal tissue are 0.45 W / (m·℃), 0.43 W / (m·℃), and 0.42 W / (m·℃), and the off-diagonal elements are approximately 0.02 W / (m·℃).

[0101] The thermal conductivity at grid point (i+1 / 2, j, k) is calculated by taking the arithmetic mean of the thermal conductivities at grid points (i, j, k) and (i+1, j, k). If grid point (i, j, k) is located at the boundary of the hemorrhage area, the thermal conductivity value is taken as the weighted average of the thermal conductivities of adjacent grid points, with the weight coefficient being inversely proportional to the distance from each grid point to the boundary. This interpolation method ensures a smooth transition of thermal conductivity at tissue boundaries and avoids numerical instability.

[0102] Substitute the discrete terms of the time derivative, the discrete terms of the space derivative, and the discrete values ​​of the thermal conductivity into the heterogeneous heat diffusion equation to establish a differential equation. In this embodiment, the left side of the differential equation is the product of the discrete terms of the time derivative and the density and the specific heat capacity, and the density is 1060 kg / m 3 , the specific heat capacity is 3900 J / (kg·℃). 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 is 650 W / m in the bleeding area. 3 , in normal tissue area the value is 420 W / m 3 .

[0103] Construct boundary conditions at the boundary positions of the differential grid. For the temperature boundary, directly specify the temperature value of the boundary grid point. In this embodiment, the temperature boundary of the skin surface is set based on the infrared thermal imaging measurement results, and the temperature value range is 32°C to 36°C. For the heat flow boundary, specify the heat flux density on the boundary surface. On the deep tissue boundary surface, the heat flux density is set to 20 W / m 2 , pointing to the inside of the body. For the convection boundary, the heat flow is calculated based on the ambient temperature and the heat transfer coefficient. On the skin surface, the ambient temperature is set to 25℃ and the heat transfer coefficient is set to 12 W / (m 2 ·℃).

[0104] 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.

[0105] 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.

[0106] 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.

[0107] 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.

[0108] 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.

[0109] like Figure 2As shown in the figure, the computation time comparison of the three different methods at various grid sizes is shown on a logarithmic scale. The alternating direction implicit method of this technical solution shows significant performance advantages at all grid sizes. 3 On a grid, the explicit difference method takes 1734.5 seconds, the standard implicit method takes 912.8 seconds, and this solution takes only 229.4 seconds, making it 7.6 times and 4.0 times faster, respectively. Of particular note, the performance advantage of this solution becomes more pronounced as the grid size increases, demonstrating that this method has significant computational efficiency advantages when processing high-resolution three-dimensional temperature field reconstructions. This high efficiency enables near-real-time temperature field analysis in clinical settings.

[0110] 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, blood flow velocity vectors are 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:

[0111] 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;

[0112] 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;

[0113] 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;

[0114] 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;

[0115] 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.

[0116] In a specific embodiment, a deep temperature distribution map and a tissue layer map of the patient are obtained, wherein the deep temperature distribution map can be obtained by a thermoacoustic imaging device and represented as a two-dimensional image of 256×256 pixels with a pixel value range of 20°C to 42°C; the tissue layer map is obtained by a medical tomography device and represented as a two-dimensional image of the same size, in which different tissue types are represented by different grayscale values.

[0117] The deep temperature distribution map is normalized, mapping the temperature values ​​to a range of 0 to 1. Specifically, for the original temperature value T(x,y), the normalized temperature value is calculated as (T(x,y) - Tmin) / (Tmax - Tmin), where Tmin = 20°C and Tmax = 42°C. A Gaussian filter is applied to the normalized temperature distribution map, and convolution is performed using a Gaussian kernel with a standard deviation of 1.5 to obtain a smoothed temperature field. The temperature field gradient is then calculated, with the gradient components in the x and y directions calculated separately to form a gradient vector field.

[0118] Adaptive threshold segmentation is performed on the tissue layer image, using the OTSU method to automatically determine the optimal threshold. In this example, thresholds of 85, 120, and 165 (grayscale range 0-255) are set for different tissue types, such as muscle, fat, and blood vessels, respectively. After segmentation, a binary tissue interface image is generated. The Sobel operator is then used to detect tissue interface edges, and the normal vector at the edge is calculated to obtain the tissue interface normal vector field.

[0119] The temperature field gradient and the tissue interface normal vector are combined to form the initial eigenvector. For each pixel (x, y) in the image, the temperature gradient vector (Gx, Gy) and the tissue interface normal vector (Nx, Ny) are extracted and combined to form a four-dimensional eigenvector (Gx, Gy, Nx, Ny). To convert to a complex-valued representation, the temperature gradient is used as the real part and the tissue interface normal vector as the imaginary part, forming a complex-valued eigenvector (Gx+jNx, Gy+jNy).

[0120] A complex-valued neural network was constructed to process the initial eigenvector. This network employs a three-layer structure consisting of an input layer, a hidden layer, and an output layer. The input layer has two nodes, corresponding to the two components of the complex-valued eigenvector; the hidden layer has 16 nodes; and the output layer has four nodes, corresponding to the elements of the output feature matrix. The network uses a piecewise linear unit function as its activation function, defined as follows: when the real 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; and when the input is greater than or equal to 1, the output is 1.

[0121] The network was trained using stochastic gradient descent with a batch size of 64, an initial learning rate of 0.01, and 200 training epochs. A mean squared error (MSE) loss function was used, and training was terminated when an error of 0.05 or less was achieved on the validation set. After training, the initial eigenvectors were input into the network to generate the thermal-acoustic coupling feature matrix, a 2×2 complex-valued matrix representing the thermal-acoustic coupling coefficients.

[0122] The Doppler-shifted signal is processed by Hilbert transforming the acquired Doppler-shifted signal (sampling rate: 10kHz, number of sampling points: 1024) to obtain a complex analytic signal. Specifically, the original frequency-shifted signal s(t) is first Fourier transformed to obtain a spectrum. The negative frequency components are then set to zero. Finally, an inverse Fourier transform is performed to obtain a complex analytic signal s_a(t) = s(t) + j·H[s(t)], where H[s(t)] is the Hilbert transform of s(t).

[0123] 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 that position to obtain the three-dimensional blood flow velocity vector field v(x, y, t).

[0124] The velocity divergence distribution is calculated based on the three-dimensional blood flow velocity vector field. Specifically, the velocity vector divergence value is calculated for each point (x, y). Regions with higher velocity divergence values ​​typically correspond to blood vessel boundaries. A divergence threshold is set at 0.15, and points with divergence values ​​greater than the threshold are extracted as the blood vessel boundary point set. In this example, approximately 2500 boundary points are obtained.

[0125] The set of vascular boundary points was fed into the cubic Marchi cube interpolation algorithm to construct a smooth vascular surface. The algorithm interpolated control points between adjacent boundary points to ensure C2 continuity of the generated surface. Specifically, the tension parameter was set to 0.5 and the deviation parameter to 0.0, resulting in a smoothed vascular surface model containing approximately 15,000 surface meshes.

[0126] In the range of blood vessel surface, radial basis function is used to interpolate the three-dimensional blood velocity vector field. Gaussian radial basis function is selected, and its form is exp(-r 2 / σ 2 ), where r is the spatial distance and σ is the shape parameter, set to 3.0. For any point within a blood vessel, the velocity value is calculated based on the weighted average of surrounding known velocity vectors, with weights determined by radial basis functions. 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 of blood flow within the vessel.

[0127] In this embodiment, a complete process of generating a high-precision blood flow velocity field distribution map is successfully implemented, starting from a deep temperature distribution map and tissue layer map, constructing a thermoacoustic coupling feature matrix via a complex-valued neural network, and ultimately generating a high-precision blood flow velocity field distribution map. By fusing the normalized temperature gradient with the tissue interface normal vector and processing it using a three-layer complex-valued neural network, the weak signal characteristics of thermoacoustic coupling are effectively captured, improving the integrity and discrimination of the feature representation. A tensor product operation is performed on the complex analytical Doppler signal and the thermoacoustic feature matrix to directly output a three-dimensional blood flow velocity vector field, enabling high-precision depiction of complex hemodynamics. By determining the vascular boundary point set based on the velocity divergence distribution and introducing a cubic Marchi cube interpolation algorithm, a smooth and continuous three-dimensional vascular surface model can be reconstructed without human intervention. Radial basis functions are used to interpolate the velocity vector field within the vascular surface, generating a continuous and smooth blood flow velocity field distribution map, which facilitates subsequent hemodynamic analysis and simulation. The full-link algorithm design, from temperature field and tissue layering to blood flow field interpolation, enables seamless data transfer and processing, improving the system's automation and robustness to noise and variation.

[0128] In an optional embodiment, calculating the velocity divergence distribution based on 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 the blood vessel surface includes:

[0129] The three-dimensional blood flow velocity vector field is decomposed into a scale pyramid to obtain velocity vector fields at multiple scale levels. The velocity divergence distribution of the velocity vector field is calculated separately to determine the multi-scale divergence distribution. Based on the multi-scale divergence distribution, feature enhancement is performed through average pooling and multi-layer perceptron to obtain the enhanced velocity divergence distribution.

[0130] The enhanced velocity divergence distribution is constructed as a graph structure including a vertex set, an edge set, and an adjacency matrix, and the graph structure is input into a graph attention network to perform feature iterative update to obtain blood vessel boundary features;

[0131] Extracting a skeleton blood vessel boundary point set based on the blood vessel boundary features, and recursively refining it to obtain a complete blood vessel boundary point set;

[0132] constructing a joint loss function including a topological consistency constraint term and a geometric constraint term, optimizing the complete blood vessel boundary point set based on the joint loss function, and determining an optimized blood vessel boundary point set;

[0133] The optimized vascular boundary point set is divided into local surface patches and a parameterized grid is established. The control vertices of the parameterized grid are calculated. A cubic Marchi basis function is constructed based on the control vertices. A tensor product operation is performed with the position information of the control vertices to obtain a local surface reconstruction result. After applying curvature continuity constraints and boundary continuity constraints to the local surface reconstruction result, the vascular surface is spliced.

[0134] In a specific embodiment, when performing a scale pyramid decomposition on the three-dimensional blood flow velocity vector field, a Gaussian filter and downsampling operation are used to construct a multi-scale representation. In this embodiment, five scale levels are constructed, and the scale parameters are 1.0, 2.0, 4.0, 8.0, and 16.0, respectively. The corresponding spatial resolutions are 1 times, 1 / 2 times, 1 / 4 times, 1 / 8 times, and 1 / 16 times 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 level of the scale pyramid are 0.1 mm, 0.2 mm, 0.4 mm, 0.8 mm, and 1.6 mm, respectively. Through multi-scale decomposition, large-scale vascular structure and small-scale vascular branch information can be captured simultaneously.

[0135] When calculating the velocity divergence distribution of the velocity vector field, the divergence value represents the net outflow of fluid per unit volume. 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 large positive divergence value. For the velocity vector field at each scale level, the velocity difference between adjacent grid points in three directions is calculated and divided by the grid spacing in the corresponding direction to obtain the velocity divergence value. In this embodiment, the divergence value range of the normal blood vessel 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 , up to 2.0 s -1 .

[0136] After determining the multi-scale divergence distribution, feature enhancement was performed using average pooling and a multi-layer perceptron. Average pooling used a 3×3×3 pooling window with a step size of 1 to smooth the divergence distribution and reduce the effects of noise. The multi-layer perceptron consists of a three-layer neural network. The number of nodes in the input layer is equal to 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 the input layer. The Reluctant Unit (ReLU) function was used as the activation function to enhance the network's nonlinear expression capabilities. During training, the learning rate was set to 0.001, the Adam optimizer was used, and the number of training epochs was 100. This feature enhancement improved the contrast of the divergence distribution, enhancing the distinction between hemorrhage and normal areas, and resulting in an enhanced velocity divergence distribution.

[0137] When constructing the enhanced velocity divergence distribution as a graph, each voxel is treated as a vertex. The vertex eigenvector contains the divergence value of that point and its divergence value at various scales, totaling six components. The edges of the graph connect adjacent voxels in space, with a neighborhood radius set to a distance of two voxels. On average, each vertex is connected to approximately 26 neighboring vertices. The adjacency matrix is ​​a sparse matrix, and the non-zero elements represent the similarity between vertices. This similarity is calculated as an exponential function of the cosine similarity of the vertex eigenvectors, with the similarity ranging from 0 to 1.

[0138] The graph structure is input into the graph attention network for iterative feature updates. The graph attention network consists of three layers of graph convolution. Each layer uses a multi-head attention mechanism with eight heads, an input feature dimension of 6, a hidden layer feature dimension of 32, and an output feature dimension of 16. Attention weights are normalized using a softmax function to ensure that the sum of the weights of all neighboring vertices is 1. After each layer of graph convolution, residual connections and batch normalization are added to improve network stability. Through feature propagation through three layers of graph convolution, each vertex integrates information within its local neighborhood, more accurately representing vascular boundary features. For hemorrhagic regions, since their features differ significantly from those of normal vascular regions, a unique representation pattern is formed during the feature update process, which facilitates subsequent boundary extraction.

[0139] To extract the skeleton vessel boundary point set based on the vessel boundary features, a threshold segmentation and thinning algorithm is used. First, the feature map is binarized, and the threshold is set to the feature mean plus 0.8 times the standard deviation to obtain an initial vessel region mask. A distance transform is applied to the mask to determine the shortest distance from each point to the background. Local maxima of the distance transform are marked as skeleton points. In this embodiment, a typical skeleton vessel network contains 500-1000 skeleton points, varying depending on the complexity of the vascular network.

[0140] The skeletonized vascular boundary point set was recursively refined to obtain a complete vascular boundary point set. Recursive refinement employed an adaptive mesh density strategy, with initial boundary point spacing of approximately 1 mm. This density was increased in regions of greater curvature, achieving a minimum inter-point distance of 0.1 mm. During refinement, the positions of newly added points were determined using cubic spline interpolation, while also adjusting for local curvature information. Five recursive refinement iterations were performed, increasing the boundary point density by approximately 10-fold in regions of high curvature, ultimately resulting in a complete vascular boundary point set consisting of 5,000–10,000 points.

[0141] A joint loss function consisting of topological consistency and geometric constraints was constructed to optimize the complete set of vascular boundary points. The topological consistency constraint ensures that the connectivity and branching structure of the vascular network remain unchanged, with a weight coefficient of 0.6. The geometric constraints include a curvature smoothing constraint and a volume preservation constraint, with weight coefficients of 0.3 and 0.1, respectively. The curvature smoothing constraint uses the Laplace operator to calculate the discrete curvature, limiting the gradient of the curvature change. The volume preservation constraint ensures that the total vascular volume does not change by more than 3% before and after optimization. The optimization algorithm uses the L-BFGS method, with a maximum number of iterations of 200 and a convergence threshold of 0.0001. Through optimization, the boundary point set achieves smoothness and precision of geometry while maintaining the topological structure.

[0142] The optimized vascular boundary point set is divided into local surface patches and a parameterized mesh is created. The local surface patch division is based on the region growing method, with the vascular branch points as the boundaries, decomposing the vascular network 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 parameterized mesh is created with a mesh resolution of 16 circumferential points. The number of axial sampling points is adaptively set based on the vessel length, with approximately 8 sampling points per millimeter. Parameterization uses cylindrical coordinate mapping to map the three-dimensional spatial surface into a two-dimensional parameter space.

[0143] The control vertices of the parameterized mesh are calculated. A 4×4×4 grid of control points is set for each local surface patch. The initial positions of the control points are obtained by least-squares fitting the set of points on the parameterized mesh. A cubic Marchi basis function of order 3 is constructed based on the control vertices, using open uniform node vectors. A tensor product operation is performed on the cubic Marchi basis function and the control vertex position information to obtain a parametric representation of the local surface. In this embodiment, each local surface patch is defined by 64 control points, providing sufficient degrees of freedom to express complex vascular shapes.

[0144] Curvature continuity and boundary continuity constraints are imposed on the local surface reconstruction results. Curvature continuity ensures that adjacent patches have continuous curvature changes at their junctions, technically achieved by matching the first two derivatives of the boundaries of adjacent patches. Boundary continuity ensures that adjacent patches completely overlap at their shared boundaries, technically achieved by sharing boundary control points. Both constraints are incorporated into the objective function using the Lagrange multiplier method, with a weight coefficient of 0.5. After constraint optimization, all local patches are spliced ​​together to obtain a complete and smooth vessel surface reconstruction.

[0145] In the existing technology, blood vessel surface reconstruction mainly adopts the March cube or direct triangulation method, which cannot effectively handle the complex topological changes of the bleeding area. The method of this embodiment starts from the hemodynamic characteristics, uses the divergence field to extract the blood vessel boundary, enhances the boundary features through the graph attention network, and combines geometric constraints and topological constraints to reconstruct a smooth and continuous blood vessel surface. Compared with traditional technologies, this method improves the boundary accuracy by 30%, improves the bleeding area detection rate by 25%, and reduces the average error of the reconstructed surface from 0.5mm to 0.15mm. In particular, for the reconstruction of small blood vessel branches and complex topological structures around bleeding points, the method of this embodiment shows significant advantages, providing reliable technical support for the clinical precise positioning of traumatic bleeding.

[0146] like Figure 3 As shown in the figure, the curvature continuity of the vessel surfaces reconstructed using three different methods was evaluated. The horizontal axis represents the normalized position (0-1) of the vessel surface points, and the vertical axis represents the average curvature variation. The curvature variation of the proposed method (blue line) remains stable across the entire surface, ranging from 0.02 to 0.08, with an average of 0.043. The curve is continuous and smooth, with fluctuations of only 0.032 in the critical region (normalized positions 0.4-0.6, corresponding to branch intersections), indicating a high degree of geometric continuity in the reconstructed surface. The curvature variation of the March cube method (red line) averages 0.176, with a sharp peak of 0.412 at the branch intersection (near 0.5), indicating a significant geometric discontinuity. The curvature variation of the direct triangulation method (green line) averages 0.231, with irregular fluctuations at multiple locations, reaching a maximum of 0.538 (at position 0.53). A distinct step-like change is observed at surface junctions (positions 0.25 and 0.75). The mean square errors of the curvature change rates of the three methods are 0.0019 (this technical solution), 0.0187 (March cube method) and 0.0226 (direct triangular mesh method), respectively. This shows that this technical solution effectively ensures the smooth transition of the surface and reduces geometric artifacts by using cubic March basis functions and imposing curvature continuity constraints. It is particularly suitable for medical application scenarios that require high-precision fluid mechanics analysis.

[0147] In an optional embodiment, a variational particle filter method is used to dynamically track the local blood flow in the time series mapping relationship, calculate the corresponding time-varying coefficient, and construct a state equation to obtain the dynamic bleeding parameters, including:

[0148] Initializing multiple particles according to the time series mapping relationship and assigning initial weights, each particle corresponding to a state of local blood flow;

[0149] Establishing an observation likelihood function to describe the relationship between the cumulative bleeding volume and the local blood flow, and establishing a state transition probability function to describe the change pattern of the local blood flow; combining the observation likelihood function and the state transition probability function to construct a variational objective function;

[0150] updating the particle weights through iterative optimization based on the variational objective function; calculating the optimal estimate of the local blood flow according to the updated particle weights; performing importance resampling to eliminate particle degradation and obtain a dynamic change result of the local blood flow;

[0151] Analyzing the blood flow variation patterns in different time periods based on the dynamic variation results of the local blood flow; calculating the scale coefficient, delay coefficient, and attenuation coefficient of the hemodynamic characteristics based on the variation patterns, and performing time series smoothing to determine the smoothing coefficients;

[0152] A state transfer matrix of a nonlinear dynamic model is constructed using a smoothing coefficient, and the system noise term is combined with the state transfer matrix to construct a state equation of the dynamic characteristics of bleeding. The state equation of the dynamic characteristics of bleeding is solved to calculate the bleeding rate, cumulative bleeding volume and degree of tissue damage, and determine the dynamic parameters of bleeding.

[0153] In one specific embodiment, multiple particles are initialized and assigned initial weights based on a time-series mapping relationship. Assuming the local blood flow state is X, 200 particles are initialized, each representing a possible blood flow state. Each particle P(i), where i = 1, 2, ..., 200, is assigned the same initial weight w(i) = 1 / 200. For example, if the initial estimated mean blood flow is 5 ml / min, particles can be generated around this value to follow a normal distribution with a standard deviation of 1 ml / min.

[0154] 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 bleeding volume Z and the local blood flow X. In practical applications, a Gaussian distribution model can be used, in which 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 of blood flow, 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 at the previous moment was 4.5 ml / min, the blood flow at the current moment may fluctuate within the range of 4.3-4.7 ml / min.

[0155] Based on the two functions above, a variational objective function is constructed. This function combines the observation likelihood and state transition probability to evaluate the credibility of each particle. The design of the variational objective function considers the balance between data fit and model prior knowledge, optimizing the particle distribution by minimizing the difference in information entropy.

[0156] Particle weights are updated through iterative optimization. At each time step t, a new weight for each particle is calculated based on the variational objective function. In this implementation, 20 iterations are used, with particle weights adjusted via gradient descent in each iteration, using a learning rate of 0.05. For example, if a particle's blood flow at t = 10 minutes is 6.2 ml / min, and the observed data indicates a more likely value of 6.0 ml / min, the particle's weight will be reduced accordingly.

[0157] The optimal estimate is calculated based on the updated particle weights. The estimated local blood flow is calculated using the weighted average method, which is the sum of all particle values ​​multiplied by their corresponding weights. For example, if five particles have values ​​of 5.8, 6.0, 6.1, 5.9, and 6.2 ml / min, and weights of 0.15, 0.25, 0.3, 0.2, and 0.1, respectively, the optimal estimate is 6.01 ml / min.

[0158] Importance resampling is performed to eliminate particle degeneration. Resampling is triggered when the effective number of particles falls below a threshold of 100. The resampling process generates a new particle set based on particle weights, with particles with larger weights having a higher probability of being selected and replicated. For example, a particle with a weight of 0.3 might be replicated into three identical particles with a weight of 0.1. The resulting particle set better represents the actual distribution of local blood flow.

[0159] Based on the dynamic changes in local blood flow, the patterns of change over different time periods were analyzed. The entire observation period was divided into multiple time windows, each lasting 15 minutes. For example, within 0-15 minutes after trauma, blood flow might rapidly increase from 2 ml / min to 8 ml / min; within 15-30 minutes, it might stabilize at 7-8 ml / min; and within 30-60 minutes, it might gradually decrease to 5 ml / min.

[0160] Hemodynamic parameters were calculated. The scaling factor represents blood flow intensity. For example, a maximum blood flow of 9 ml / min corresponds to a scaling factor of 1.8. The delay factor reflects the blood flow response time. For example, blood flow begins to change significantly 2 minutes after trauma, corresponding to a delay factor of 0.033. The attenuation factor represents the rate of blood flow reduction. For example, a blood flow decreases from 7 ml / min to 5 ml / min between 40 and 60 minutes, corresponding to an attenuation factor of 0.017.

[0161] The above coefficients are smoothed over time using a sliding window averaging method with a window size of 5 time points. For example, for a scaling coefficient of 1.75 at a certain time point, the values ​​of 1.70, 1.72, 1.75, 1.78, and 1.76 at the two preceding and following time points are averaged to obtain a smoothed value of 1.742.

[0162] The smoothed coefficients are used to construct the state transition matrix of the nonlinear dynamic model. This matrix describes how the system state transitions from one time point to the next. For example, the state transition matrix might express the relationship between the current bleeding rate and the previous bleeding rate as follows: current bleeding rate = 0.95 × previous bleeding rate + system noise. The 0.95 is the state transition coefficient, reflecting the system's memory effect.

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

[0164] Solving the dynamic state equation for bleeding characteristics allows calculation of bleeding rate, cumulative bleeding volume, and degree of tissue damage. For example, 30 minutes after trauma, the bleeding rate is 6.5 ml / min, the cumulative bleeding volume is 180 ml, and the degree of tissue damage is moderate (corresponding to a damage index of 0.65). By analyzing these parameters at multiple time points, a comprehensive assessment of bleeding status can be made to guide clinical intervention decisions.

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

[0166] In an optional embodiment, updating the weights of particles through iterative optimization based on the variational objective function includes:

[0167] Calculating the partial derivative of the variational objective function with respect to each particle weight to obtain a sensitivity eigenvalue; performing symbol extraction on the sensitivity eigenvalue and multiplying it with a power term of the sensitivity eigenvalue to obtain a weight update direction;

[0168] Calculate the KL divergence of the particle weight distribution of 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;

[0169] Calculating the maximum and minimum values ​​of the updated particle weights, and performing renormalization processing on the updated particle weights to obtain renormalized particle weights;

[0170] The Euclidean distance between the particle weights of two adjacent iterations is calculated to obtain a stability index. When the stability index is less than a preset stability threshold for a consecutive preset number of times, the particle weight of the current iteration is determined to be the optimization result.

[0171] In one specific embodiment, when calculating the partial derivative of the variational objective function with respect to each particle weight to obtain the sensitivity eigenvalue, the variational objective function consists of two parts: an image similarity measurement term and a spatial smoothing regularization term. The image similarity measurement term uses normalized mutual information, and the value range is generally 0.5 to 1.5. The mutual information value in the bleeding area is approximately 0.8, and the mutual information value in the normal tissue area is approximately 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 trauma bleeding detection scenario, the number of particles is set to 1024, distributed in a detection area of ​​100mm×100mm, and the particle spacing is approximately 3mm. The partial derivative of the variational objective function with respect to each particle weight is calculated using the finite difference method, and the central difference step size is set to 0.01. For each particle, the objective function value is calculated by increasing and decreasing its weight by 0.01, and the difference between the two divided by 0.02 is the sensitivity eigenvalue of the particle weight. The sensitivity eigenvalue usually ranges from -2.0 to 2.0. A positive value indicates that increasing the particle weight will increase the objective function value, and a negative value indicates that decreasing the particle weight will increase the objective function value.

[0172] The sensitivity eigenvalue is sign-extracted and multiplied 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 power of 5 / 3 and then multiplied by its sign. This nonlinear 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 areas. For 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 amplitude of the update direction is approximately 2.2 to 3.2.

[0173] The adaptive step size is constructed by calculating the KL divergence of the particle weight distributions between two consecutive iterations. The KL divergence measures the degree of difference between two probability distributions, with larger values ​​indicating greater disparity. In the early stages of an iteration, the KL divergence is typically large, ranging from 0.5 to 1.0. As the iterations progress, the KL divergence gradually decreases, eventually converging to below 0.01. The adaptive step size is linked to the KL divergence via an exponential decay function: the initial step size multiplied by the exponential function of the KL divergence, with the exponent coefficient set to -2. With an initial step size of 0.2, when the KL divergence is 1.0, the adaptive step size is approximately 0.027; when the KL divergence is 0.1, the adaptive step size is approximately 0.164; and when the KL divergence is 0.01, the adaptive step size is approximately 0.196.

[0174] The adaptive step size is multiplied by the weight update direction to determine the weight update amount. This weight update amount is then added to the particle weight of the current iteration to obtain the updated particle weight. For a typical bleeding detection scenario, the absolute value of the weight update amount ranges from 0.001 to 0.5, with an average of approximately 0.05. The update amount is relatively large in the early stages of the iteration and gradually decreases as the iteration progresses. After 30 iterations, the average absolute value of the update amount is approximately 0.01, indicating that the optimization process is nearing convergence.

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

[0176] The stability index is obtained by calculating the Euclidean distance between the particle weights of two adjacent iterations. The Euclidean distance is calculated as the square root of the sum of the squares of the differences in all particle weights. In the early stages of the iteration, the stability index is usually between 0.1 and 0.2; as the iteration proceeds, the stability index gradually decreases. The preset stability threshold is set to 0.001, and the preset number of times is set to 5. When the stability index of 5 consecutive iterations is less than 0.001, the optimization process is considered to have converged, and the particle weight of the current iteration is the optimization result. In a typical trauma bleeding detection scenario, the optimization usually converges within 40 to 60 iterations.

[0177] This particle weight optimization method accurately identifies hemorrhage areas in infrared thermal imaging and ultrasound images. The optimized particle weights form distinct clusters, with high-weight particles concentrated in the hemorrhage area and low-weight particles distributed in the normal tissue area. At the boundary between the hemorrhage area and normal tissue, the particle weights exhibit a gradient distribution, helping to precisely locate the bleeding boundary. The optimized results are used in subsequent bleeding volume estimation and bleeding velocity calculations, providing important insights for clinical diagnosis.

[0178] Traditional medical image registration methods are mainly based on grayscale similarity or feature matching, and suffer from problems such as insufficient accuracy and poor robustness when processing infrared thermal imaging and ultrasonic multimodal images. Particle optimization methods commonly used 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 apply to real-time trauma bleeding detection scenarios. The particle weight optimization method proposed in this embodiment introduces a nonlinear transformation of the sensitivity eigenvalue, KL divergence adaptive step size adjustment, and Euclidean distance stability judgment mechanism. Starting from the convergence and stability of the optimization algorithm, it solves the key problems in multimodal medical image registration.

[0179] Compared with the existing technology, this embodiment enhances the accuracy of the optimization direction through the power transformation of the sensitivity eigenvalue, so that the weight update pays more attention to important areas; introduces the KL divergence adaptive step size mechanism to balance the convergence speed and stability; and uses the Euclidean distance stability index to judge the convergence condition, thereby improving the reliability of the optimization results.

[0180] like Figure 4 The figure shows a detailed demonstration of the entire multimodal medical image registration process based on particle weight optimization. The process begins with initializing particle weights and registration parameters, followed by computing a variational objective function to assess the current registration quality. The core optimization step calculates sensitivity eigenvalues—the partial derivatives of the variational objective function with respect to each particle weight—and applies sign extraction and power nonlinear transformations to determine the weight update direction. A KL divergence calculation module analyzes the differences in particle weight distributions between two consecutive iterations and constructs an adaptive step-size mechanism to balance convergence speed and stability. The weight update and renormalization steps ensure the effectiveness of the weight update, maintaining the weights within a reasonable range through the calculation of minimum and maximum constraints and renormalization. The parallel computing module on the right side of the flowchart demonstrates the simultaneous processing of infrared thermal and ultrasound images, while the parameter monitoring module tracks key parameters such as the KL divergence (0.0023-0.0158) and Euclidean distance (0.0067-0.0344) in real time. The Euclidean distance stability metric is used to determine algorithm convergence. The optimization results are output when the number of consecutive iterations is less than a preset stability threshold.

[0181] In an alternative embodiment, designed for emergency medical rescue scenarios, a lightweight configuration is employed, consisting of a portable infrared thermal imager (weighing less than 1 kg) and a handheld ultrasound probe (weighing less than 0.5 kg), both wirelessly connected to a processing terminal. The calibration plate is designed to fold, measuring 200 mm x 200 mm when unfolded and only 10 mm thick when folded, making it easy to carry in the field. The entire system can be packed into a medical first aid kit and carried by a single person.

[0182] Specially optimized for field rescue conditions, the infrared thermal imaging acquisition process incorporates an ambient temperature compensation algorithm to address the significant variations in ambient lighting. By setting a temperature reference point on the calibration plate, the imaging parameters are dynamically adjusted to ensure temperature measurement accuracy within ±0.2°C within an ambient temperature range of 0°C to 40°C.

[0183] To address the issue of limited power supply in the wild, the processing algorithm has been optimized for computational complexity. The adaptive shearlet transform uses a sparse sampling strategy, performing calculations only on key frequency points, reducing the computational complexity by approximately 70%. Adaptive meshing technology is used to solve the heterogeneous heat diffusion equation, with mesh density increased in boundary regions and sparse meshes in flat regions. This improves computational efficiency by a factor of three while maintaining accuracy. The structure of the complex-valued neural network has also been streamlined, reducing the number of hidden layer nodes from 16 to 8, reducing the model size by 60%, and reducing the single inference time from 120ms to 45ms. These optimizations enable the system to achieve near-real-time processing on mid-range mobile processing terminals, with a battery life of up to 4 hours.

[0184] In the trauma bleeding detection process in rescue scenarios, the entire detection process is divided into three stages: rapid screening, precise positioning, and continuous monitoring. In the rapid screening stage, medical staff use infrared thermal imagers to scan the victim's entire body and automatically identify areas with abnormal temperatures. The entire process takes no more than 30 seconds. In the precise positioning stage, the ultrasound probe is placed in the area with abnormal temperature, and multimodal registration and bleeding feature extraction are automatically completed to provide the bleeding location, depth, and preliminary quantitative results. This stage takes about 90 seconds. In the continuous monitoring stage, the system automatically updates the dynamic bleeding parameters every 5 minutes, including bleeding rate, cumulative bleeding volume, and tissue damage score, and provides the trend of risk level changes.

[0185] Compared with traditional inspection methods in emergency medical rescue scenarios, the technology of this embodiment does not require large-scale equipment support and is suitable for natural disaster sites such as earthquakes and floods, as well as battlefield environments. The detection process is fast and efficient, taking no more than 2 minutes from discovering the casualty to completing the initial bleeding assessment, saving about 80% of the time compared with traditional methods. The dynamic bleeding parameters provided directly support injury classification and transfer priority judgment, which helps to rationally allocate medical resources in the case of large-scale casualties. The continuous monitoring function enables rescue personnel to promptly detect changes in the injury condition and make medical intervention adjustments during the transfer process.

[0186] This example demonstrates the significant advantages of this technology in emergency medical rescue scenarios. It can provide reliable trauma bleeding detection under resource-constrained conditions, support rapid medical decision-making, and improve rescue efficiency and treatment success rates. Through lightweight design and algorithm optimization, it addresses practical application challenges in field rescue, providing a new technical means to enhance emergency medical rescue capabilities.

[0187] In an optional embodiment, in search and rescue scenarios following natural disasters such as earthquakes and mudslides, rescuers often face the challenge of rapidly detecting and classifying scattered casualties over a large area. The method and system are integrated into a lightweight unmanned aerial vehicle (UAV) platform to construct a flexible and maneuverable aerial rescue detection system. The system includes: 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 ultrasound 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 with a processing power of 8 TOPS and a power consumption of only 15W; and a wireless transmission module that supports 5G communication and has a transmission delay of less than 100ms.

[0188] In rescue applications, the drone first uses infrared thermal imaging technology to conduct a large-scale scan at an altitude of 200-300 meters to identify areas with abnormal surface temperature. For suspected human heat sources identified, the drone lowers 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. Once the person is confirmed to be trapped, the drone lands near the injured person, and rescue personnel control or remotely operate the robotic arm to use an ultrasound probe to examine the possible trauma sites of the injured person. The collected infrared thermal imaging and ultrasound images are processed using the method of the present invention. The deep temperature distribution map can detect internal bleeding 2-6 cm below the body surface with a sensitivity of 92%. The blood flow velocity field distribution map is used to determine the location and range of bleeding, with a spatial accuracy of better than 3mm. The severity of the injury is automatically classified based on the dynamic parameters of bleeding, and the injured are divided into three levels: those requiring urgent treatment (within 30 minutes), those requiring less urgent treatment (within 2 hours), and those who can be delayed (within 6 hours). The estimated bleeding volume is calculated in real time. When the cumulative bleeding volume exceeds 800ml, the system automatically classifies the injured person as critically ill and gives priority treatment.

[0189] In another optional embodiment, in battlefield environments, casualty triage and rapid treatment face multiple challenges, including harsh environments, personnel shortages, and high pressure. This embodiment integrates the method and system into battlefield first aid equipment to develop a portable combat trauma detection system; the system is designed as a rugged structure that is waterproof and dustproof (IP67), shock-resistant (MIL-STD-810G), and weighs no more than 3.5 kg. It includes: a high-resolution infrared thermal imager (384×288 pixels) that can operate stably at ambient temperatures of -20°C to 50°C; a military-grade ultrasound probe that supports B-type, color Doppler, and power Doppler modes, with a penetration depth of up to 15 cm; and an embedded computing unit with a low-power design that can operate continuously for 8 hours on a single charge. In combat trauma treatment applications, rapid assessment of the wounded is achieved. Field medical personnel use infrared thermal imagers to perform full-body scans of the wounded, which takes about 30 seconds and automatically identifies hot spots with abnormal temperatures. For the suspected bleeding areas identified, an ultrasound probe is used for local scanning, with an acquisition time of about 45 seconds. The algorithms of the method of the present invention, such as spatiotemporal registration, adaptive shear wave transform, thermoacoustic coupling feature extraction, and variational particle filtering, are executed in a processing time of less than 15 seconds. A dynamic bleeding parameter report is generated, including: three-dimensional positioning of the bleeding position (accuracy ±5mm), bleeding rate (error ≤10%), cumulative bleeding volume (error ≤15%), estimated bleeding duration, and degree of tissue damage.

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

[0191] In pre-hospital emergency applications, this workflow improves the efficiency of trauma assessment. After arriving at the scene, the system is used to perform a preliminary scan of the injured person, including routine vital sign measurements and infrared thermal imaging of the trauma site. Based on the infrared thermal imaging results, the system automatically guides the rescuer to perform an ultrasonic scan of the suspicious area. After the collected bimodal data 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 (MTS) based on the dynamic parameters of bleeding and provides treatment recommendations.

[0192] The trauma bleeding detection method system based on infrared thermal imaging and ultrasound multimodality in an embodiment of the present invention includes:

[0193] The first unit is used to collect infrared thermal imaging images and ultrasound images of the wound site, and perform spatiotemporal registration of the images through a calibration plate to obtain registered image data;

[0194] 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 thermal diffusion coefficient, and perform pulse compression and phase correction on the acoustic coefficients to obtain the tissue layer map;

[0195] The third unit is used to construct a thermoacoustic coupling characteristic matrix based on the deep temperature distribution map and tissue layer map using a complex-valued neural network, extract the blood flow velocity vector from the Doppler frequency shift signal based on the thermoacoustic coupling characteristic matrix, reconstruct the vascular structure using Marchi cube interpolation, and generate a blood flow velocity field distribution map;

[0196] The fourth unit is used to calculate the local blood flow flux based on the blood flow velocity field distribution map, perform volume integration based on the anisotropic thermal diffusion coefficient to obtain the cumulative bleeding volume, and establish a time series mapping relationship between the local blood flow flux and the cumulative bleeding volume;

[0197] The fifth unit is used to dynamically track the local blood flow in the time series mapping relationship using the variational particle filter method, calculate the corresponding time-varying coefficients, and construct the state equation to obtain the dynamic parameters of bleeding;

[0198] Unit 6 is used to calculate the duration of bleeding and the degree of tissue damage based on the dynamic parameters of bleeding, and to determine the risk level and treatment plan.

[0199] According to a third aspect of an embodiment of the present invention, an electronic device is provided, including:

[0200] processor;

[0201] a memory for storing processor-executable instructions;

[0202] The processor is configured to call the instructions stored in the memory to execute the aforementioned method.

[0203] According to a fourth aspect of an embodiment of the present invention, a computer-readable storage medium is provided, on which computer program instructions are stored. When the computer program instructions are executed by a processor, the method described above is implemented.

[0204] The present invention may be a method, an apparatus, a system and / or a computer program product. The computer program product may include a computer-readable storage medium carrying computer-readable program instructions for executing various aspects of the present invention.

[0205] 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 it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the above embodiments, or replace some or all of the technical features therein with equivalents. However, 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: include: Collect infrared thermal imaging images and ultrasound images of the wound site, perform spatiotemporal registration of the images using a calibration plate, and obtain registered image data; Adaptive shear wave transform is performed on the registered image data to obtain temperature field coefficients and acoustic coefficients. The temperature field coefficients are substituted into the heterogeneous heat diffusion equation to calculate the deep temperature distribution map and obtain the anisotropic thermal diffusion coefficient. Pulse compression and phase correction are performed on the acoustic coefficients to obtain the tissue layer map. Based on 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 cube interpolation to generate a blood flow velocity field distribution map. The local blood flow flux is calculated based on the blood flow velocity field distribution map, and the cumulative bleeding volume is obtained by volume integration combined with the anisotropic thermal diffusion coefficient. The time series mapping relationship between the local blood flow flux and the cumulative bleeding volume is established. 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, and construct the state equation to obtain the dynamic bleeding parameters. The duration of bleeding and the degree of tissue damage are calculated based on the dynamic parameters of bleeding to determine the risk level and treatment plan.

2. The method according to claim 1, characterized in that Adaptive shear wave transform is performed on the registered image data to obtain the temperature field coefficient and acoustic coefficient. The temperature field coefficient is substituted into the heterogeneous heat diffusion equation to calculate the deep temperature distribution map and obtain the anisotropic thermal diffusion coefficient. Pulse compression and phase correction are performed on the acoustic coefficient to obtain the tissue layer map, including: Performing an adaptive shearlet transform on the registered image data, wherein the adaptive shearlet transform jointly decomposes the registered image data by setting a direction-selective filter group, constructs a shearlet basis function according to a shearing parameter and a translation parameter, obtains a coefficient matrix by performing a two-dimensional integral operation on the shearlet basis function, calculates an angular direction parameter based on an inverse tangent function of the shearing parameter, and optimizes the angular direction parameter to achieve multi-directional adaptive decomposition and obtain a temperature field coefficient and an acoustic coefficient; Substituting the temperature field coefficient into a 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; calculating the coefficient term by multiplying the tissue density term by the specific heat capacity term; calculating the diffusion term by the inner product of the thermal conductivity tensor term and the temperature gradient; substituting the coefficient term, the diffusion term, and the heat source term into the heterogeneous heat diffusion equation; solving the heterogeneous heat diffusion equation using a finite difference method to obtain a deep temperature distribution map; and calculating the anisotropic thermal diffusion coefficient based on the ratio of maximum thermal conductivity to 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.

3. The method according to claim 2, characterized in that The finite difference method is used to solve the heterogeneous heat diffusion equation to obtain the deep temperature distribution map including: 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, 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, characterized in that 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 three-dimensional blood flow velocity vector field is decomposed into a scale pyramid to obtain velocity vector fields at multiple scale levels. The velocity divergence distribution of the velocity vector field is calculated separately to determine the multi-scale divergence distribution. Based on the multi-scale divergence distribution, feature enhancement is performed through average pooling and multi-layer perceptron to obtain the enhanced velocity divergence distribution. The enhanced velocity divergence distribution is constructed as a graph structure including a vertex set, an edge set, and an adjacency matrix, and the graph structure is input into a graph attention network to perform feature iterative update to obtain blood vessel boundary features; Extracting a skeleton blood vessel boundary point set based on the blood vessel boundary features, and recursively refining it to obtain a complete blood vessel boundary point set; constructing a joint loss function including a topological consistency constraint term and a geometric constraint term, optimizing the complete blood vessel boundary point set based on the joint loss function, and determining an optimized blood vessel boundary point set; The optimized vascular boundary point set is divided into local surface patches and a parameterized grid is established. The control vertices of the parameterized grid are calculated. A cubic Marchi basis function is constructed based on the control vertices. A tensor product operation is performed with the position information of the control vertices to obtain a local surface reconstruction result. After applying curvature continuity constraints and boundary continuity constraints to the local surface reconstruction result, the vascular surface is spliced.

6. The method according to claim 1, characterized in that 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, and construct the state equation to obtain the dynamic bleeding parameters including: Initializing multiple particles according to the time series mapping relationship and assigning initial weights, each particle corresponding to a state of local blood flow; Establishing an observation likelihood function to describe the relationship between the cumulative bleeding volume and the local blood flow, and establishing a state transition probability function to describe the change pattern of the local blood flow; combining the observation likelihood function and the state transition probability function to construct a variational objective function; updating the particle weights through iterative optimization based on the variational objective function; calculating the optimal estimate of the local blood flow according to the updated particle weights; performing importance resampling to eliminate particle degradation and obtain a dynamic change result of the local blood flow; Analyzing the blood flow variation patterns in different time periods based on the dynamic variation results of the local blood flow; calculating the scale coefficient, delay coefficient, and attenuation coefficient of the hemodynamic characteristics based on the variation patterns, and performing time series smoothing to determine the smoothing coefficients; A state transfer matrix of a nonlinear dynamic model is constructed using a smoothing coefficient, and the system noise term is combined with the state transfer matrix to construct a state equation of the dynamic characteristics of bleeding. The state equation of the dynamic characteristics of bleeding is solved to calculate the bleeding rate, cumulative bleeding volume and degree of tissue damage, and determine the dynamic parameters of bleeding.

7. The method according to claim 6, characterized in that Updating the particle weights 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 a sensitivity eigenvalue; performing symbol extraction on the sensitivity eigenvalue and multiplying it with a power term of the sensitivity eigenvalue to obtain a weight update direction; Calculate the KL divergence of the particle weight distribution of 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; Calculating the maximum and minimum values ​​of the updated particle weights, and performing renormalization processing on the updated particle weights to obtain renormalized particle weights; The Euclidean distance between the particle weights of two adjacent iterations is calculated to obtain a stability index. When the stability index is less than a preset stability threshold for a consecutive preset number of times, the particle weight of the current iteration is determined to be the optimization result.

8. A trauma bleeding detection method system based on infrared thermal imaging and ultrasound multimodality, used to implement the method according to any one of claims 1 to 7, characterized in that: include: The first unit is used to collect infrared thermal imaging images and ultrasound images of the wound site, and perform spatiotemporal registration of the images through a calibration plate to 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 thermal diffusion coefficient, and perform pulse compression and phase correction on the acoustic coefficients to obtain the tissue layer map; The third unit is used to construct a thermoacoustic coupling characteristic matrix based on the deep temperature distribution map and tissue layer map using a complex-valued neural network, extract the blood flow velocity vector from the Doppler frequency shift signal based on the thermoacoustic coupling characteristic matrix, reconstruct the vascular structure using Marchi cube interpolation, and generate a blood flow velocity field distribution map; The fourth unit is used to calculate the local blood flow flux based on the blood flow velocity field distribution map, perform volume integration based on the anisotropic thermal diffusion coefficient to obtain the cumulative bleeding volume, and establish a time series mapping relationship between the local blood flow flux and the cumulative bleeding volume; The fifth unit is used to dynamically track the local blood flow in the time series mapping relationship using the variational particle filter method, calculate the corresponding time-varying coefficients, and construct the state equation to obtain the dynamic parameters of bleeding; Unit 6 is used to calculate the duration of bleeding and the degree of tissue damage based on the dynamic parameters of bleeding, and to determine the risk level and treatment plan.

9. An electronic device, characterized in that: include: processor; a memory for storing processor-executable instructions; 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: When the computer program instructions are executed by a processor, the method according to any one of claims 1 to 7 is implemented.

Images