Real-time monitoring and control method and system for intraductal pressure in ductoscopy

Abstract

The application provides a duct-in-pressure real-time monitoring and control method and system in a mammary ductoscopy, relates to the technical field of medical equipment monitoring and control, and comprises the following steps: collecting temperature, conductivity and image data in a duct to form a multi-dimensional parameter feature matrix, performing adaptive filtering analysis to extract fluid characteristic parameters, calculating duct-in-pressure distribution features, combining with microfluid sensor array monitoring fluid flow state parameters to identify local turbulent flow areas, when the turbulent flow intensity exceeds a warning threshold, using wavelet transform analysis to calculate the turbulent flow development trend, and controlling the duct-in-pressure in a safe range through active aortic pulsation compensation technology, so that the safety and accuracy of the mammary ductoscopy are improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of medical device monitoring and control, and particularly relates to a catheter internal pressure real-time monitoring and control method and system in a breast ductoscopy examination. BACKGROUND

[0002] Breast ductoscopy is a minimally invasive technique for diagnosing breast diseases. A small endoscope is inserted through the nipple into the breast duct to observe the internal structure and lesions of the breast duct. This technique is widely used in the diagnosis of intraductal papilloma, duct ectasia and early breast cancer. During the breast ductoscopy examination, physiological saline or other perfusion fluid is injected into the breast duct to dilate the duct, facilitating the observation of the internal structure of the duct.

[0003] During the breast ductoscopy examination, the control of the internal pressure of the catheter is crucial. Proper internal pressure of the catheter can ensure a good view and examination effect, while excessive internal pressure of the catheter can cause damage to the breast tissue, exacerbate pain during the examination, and even cause fluid leakage to the surrounding tissue, resulting in complications. Traditional breast ductoscopy mainly relies on the experience of the operator to control the speed and pressure of the injected fluid, lacking precise pressure monitoring and control means.

[0004] The existing breast ductoscopy pressure control method has the following disadvantages: lack of real-time and accurate internal pressure monitoring technology of the catheter, making it difficult for the operator to accurately grasp the actual pressure state in the catheter, which can easily cause excessive or insufficient pressure, affecting the examination effect and patient safety; the existing method cannot effectively identify and predict local turbulence and pressure fluctuations in the catheter, which can cause tissue damage or decrease the quality of the examination image; traditional pressure control uses passive adjustment, lacking an intelligent active pressure regulation system, which cannot adaptively adjust the pressure according to the real-time fluid state in the catheter, especially when facing complex duct morphology or lesion conditions, it is difficult to maintain a stable pressure environment. SUMMARY

[0005] The present application provides a breast ductoscopy examination catheter internal pressure real-time monitoring and control method and system, which can solve the problems in the prior art.

[0006] In a first aspect of the present application, a breast ductoscopy examination catheter internal pressure real-time monitoring and control method is provided, comprising:

[0007] Collecting real-time temperature data, conductivity data and image data in the catheter during the breast ductoscopy examination and forming a multi-dimensional parameter feature matrix;

[0008] Adaptive filtering analysis is performed on the multi-dimensional parameter feature matrix to extract characteristic parameters of the sub-millimeter fluid particles in the fluid, and a spatial distribution density value and a motion velocity vector of the sub-millimeter fluid particles are calculated according to the characteristic parameters; a spatial distribution feature of the pressure in the conduit is calculated according to the spatial distribution density value, positions of a pressure peak point and a pressure valley point are determined, and a pressure gradient mutation region is identified; a pressure wave propagation path is tracked according to the motion velocity vector, a reflection point and an attenuation point of the pressure wave on the propagation path are identified, and a pressure wave propagation characteristic distribution map is obtained;

[0009] The pressure gradient mutation region and the pressure wave propagation characteristic distribution map are combined to monitor the pressure in the conduit in real time.

[0010] The fluid flow state parameters of each cross section of the conduit are monitored by the micro-fluid sensor array built in the conduit, the local turbulent flow region is identified according to the fluid flow state parameters and the multi-dimensional parameter feature matrix, when the turbulent flow intensity of the local turbulent flow region exceeds a warning intensity threshold, the turbulent flow development trend is calculated; according to the turbulent flow development trend and the frequency domain feature of the pressure fluctuation, the pressure fluctuation is eliminated, and the pressure in the conduit is controlled to be maintained within a safe pressure range.

[0011] The adaptive filtering analysis is performed on the multi-dimensional parameter feature matrix to extract characteristic parameters of the sub-millimeter fluid particles in the fluid, and a spatial distribution density value and a motion velocity vector of the sub-millimeter fluid particles are calculated according to the characteristic parameters, including:

[0012] The autocorrelation matrix of the multi-dimensional parameter feature matrix is calculated by adaptive Wiener filtering to obtain filtered feature data;

[0013] In the filtered feature data, an edge point is detected by using a multi-scale operator, a complete contour is obtained by performing a morphological connection process on the edge point, an equivalent diameter and a circularity coefficient are calculated according to the complete contour, and a gray scale gradient value of the contour region is calculated, and the equivalent diameter, the circularity coefficient, and the gray scale gradient value are combined to form the characteristic parameters of the sub-millimeter fluid particles;

[0014] An initial bandwidth is determined according to the characteristic parameters, the initial bandwidth is iteratively optimized by cross-validation until the density estimation error is less than a distribution fitting threshold, an optimized bandwidth parameter is obtained, and the spatial distribution density value of the sub-millimeter fluid particles is calculated by using the optimized bandwidth parameter;

[0015] Images of adjacent frames in the filtered feature data are extracted, spatial and temporal gradients of the adjacent frame images are calculated and a velocity field smooth constraint is constructed in combination with the characteristic parameters, an optimization objective function is constructed by using the spatial and temporal gradients and the velocity field smooth constraint, the motion velocity vector of the sub-millimeter fluid particles is determined by iteratively optimizing the optimization objective function by a gradient descent method until the optimization objective function value is less than a velocity field convergence threshold.

[0016] calculating a spatial distribution feature of the pressure in the conduit according to the spatial distribution density value, determining positions of the pressure peak points and the pressure valley points, and identifying the pressure gradient mutation region include:

[0017] calculating an initial pressure distribution based on the spatial distribution density value in the conduit;

[0018] calculating a Laplacian of the density, obtaining a density gradient correction term based on the Laplacian and a density gradient coefficient, superimposing the density gradient correction term on the initial pressure distribution to obtain a corrected pressure distribution;

[0019] performing a multi-scale analysis on the corrected pressure distribution, extracting local extreme points at each scale, analyzing convergence and divergence of a pressure gradient vector field corresponding to each local extreme point, marking a local extreme point with a convergent gradient vector field as a pressure peak point, and marking a local extreme point with a divergent gradient vector field as a pressure valley point;

[0020] calculating a gradient amplitude of the corrected pressure distribution at each position according to the distribution positions of the pressure peak points and the pressure valley points, calculating a direction change rate of a gradient direction between adjacent positions, dividing the gradient amplitude by an average value of gradient amplitudes in a local region to obtain an amplitude ratio, and calculating a gradient mutation index according to the amplitude ratio and the direction change rate, marking a region between adjacent pressure peak points and pressure valley points as a pressure gradient mutation region when the gradient mutation index is greater than a mutation determination threshold.

[0021] tracking a pressure wave propagation path according to the motion velocity vector, identifying a reflection point and an attenuation point of the pressure wave on the propagation path, and obtaining a pressure wave propagation characteristic distribution map include:

[0022] extracting a velocity component of a discrete point in the corrected pressure distribution, calculating a pressure wave acceleration according to velocity components of adjacent discrete points, and integrating the pressure wave acceleration along time to obtain a pressure wave propagation path;

[0023] calculating an inner product of a velocity direction at each time and a velocity direction at a next time along the pressure wave propagation path, obtaining a propagation direction change angle according to the inner product and a velocity amplitude, and marking a corresponding coordinate point as a reflection point when the propagation direction change angle is greater than a direction mutation threshold;

[0024] dividing a pressure difference of adjacent coordinate points on the pressure wave propagation path by a distance between the two points to obtain an along-path pressure gradient, dividing the along-path pressure gradient by a pressure value to obtain an attenuation coefficient, and marking a corresponding coordinate point as an attenuation point when the attenuation coefficient is greater than an attenuation characteristic threshold;

[0025] Calculate the distance from each reflection point to each attenuation point, pair the reflection points and attenuation points with distance less than a feature correlation threshold to form a feature point group, and calculate the spatial distribution density of the feature point group to obtain a propagation feature intensity;

[0026] Draw a pressure wave propagation path, mark the propagation feature intensity, mark the reflection point position and reflection direction, mark the attenuation point position and attenuation coefficient, and generate a pressure wave propagation characteristic distribution map.

[0027] Monitor the fluid flow state parameters of each cross section of the catheter through the microfluid sensor array built in the catheter, identify the local turbulent flow area according to the fluid flow state parameters and the multi-dimensional parameter feature matrix, and calculate the turbulent flow development trend when the turbulent flow intensity of the local turbulent flow area exceeds the warning intensity threshold, including:

[0028] Collect the pressure fluctuation value and flow velocity distribution value of each cross section and perform mapping operation with the multi-dimensional parameter feature matrix to obtain the strain rate tensor and vorticity tensor, calculate the second invariant of the vorticity tensor and the second invariant of the rate of change tensor respectively, and obtain the turbulent flow discriminant factor by difference, and mark the local turbulent flow area based on the turbulent flow discriminant factor;

[0029] Extract the flow velocity fluctuation component in the local turbulent flow area, calculate the ratio of the root mean square value of the flow velocity fluctuation component to the average flow velocity to obtain the turbulent flow intensity;

[0030] When the turbulent flow intensity exceeds the warning intensity threshold, extract the flow velocity fluctuation signal of the local turbulent flow area and perform wavelet transform, integrate to obtain a wavelet coefficient matrix under different scales, calculate the square of each element in the wavelet coefficient matrix to obtain the local wave energy, sum the local wave energy under the same scale to obtain the scale energy, and extract the energy distribution interval based on the scale energy;

[0031] In the energy distribution interval, the difference between the turbulent energy of adjacent time points divided by the time interval is the energy change rate, the difference between the characteristic length of adjacent scales divided by the time interval is the scale change rate, and the turbulent flow development index is obtained based on the energy change rate and the scale change rate. According to the time sequence change characteristics of the turbulent flow development index, the turbulent flow development trend is judged.

[0032] According to the turbulent flow development trend and the frequency domain characteristics of pressure fluctuation, eliminate the pressure fluctuation, and control the pressure in the catheter to maintain within a safe pressure range, including:

[0033] Extract the pressure fluctuation signal in the catheter, decompose the pressure fluctuation signal into periodic fluctuation components in the time domain, and obtain the fluctuation main frequency and its amplitude according to the Fourier series decomposition result of the periodic fluctuation component;

[0034] Set the catheter pressure target value, and the difference between the real-time detected catheter pressure value and the catheter pressure target value is taken as a pressure deviation;

[0035] According to the turbulence development trend, the main frequency and the amplitude of the fluctuation, a proportional coefficient, an integral time constant and a differential time constant are calculated, the proportional coefficient is multiplied by the pressure deviation to obtain a proportional control term, the integral time constant is multiplied by the pressure deviation integral value to obtain an integral control term, the differential time constant is multiplied by the pressure deviation change rate to obtain a differential control term, and the proportional control term, the integral control term and the differential control term are added to obtain a compensation control signal;

[0036] The compensation control signal is input into a piezoelectric ceramic diaphragm pump to generate a compensation pressure, and the catheter pressure value after compensation is detected in real time; when the catheter pressure value after compensation exceeds the upper limit value or the lower limit value of the safe pressure range, the change rate of the pressure deviation is calculated, and the compensation pressure is adjusted according to the change rate of the pressure deviation, so that the pressure value in the catheter is maintained within the safe pressure range.

[0037] In a second aspect of the embodiment of the present application, a catheter pressure real-time monitoring and control system in a breast ductoscopy examination is provided, comprising:

[0038] A first unit is configured to collect real-time temperature data, conductivity data and image data in the catheter during the breast ductoscopy examination and form a multi-dimensional parameter feature matrix;

[0039] A second unit is configured to perform adaptive filtering analysis on the multi-dimensional parameter feature matrix, extract characteristic parameters of sub-millimeter fluid particles in the fluid, calculate the spatial distribution density value and the motion velocity vector of the sub-millimeter fluid particles according to the characteristic parameters, calculate the spatial distribution characteristics of the catheter pressure according to the spatial distribution density value, determine the positions of the pressure peak point and the pressure valley point, and identify the pressure gradient mutation region; track the pressure wave propagation path according to the motion velocity vector, identify the reflection point and the attenuation point of the pressure wave on the propagation path, and obtain a pressure wave propagation characteristic distribution map;

[0040] A third unit is configured to combine the pressure gradient mutation region and the pressure wave propagation characteristic distribution map to monitor the pressure in the catheter in real time.

[0041] A fourth unit is configured to monitor the fluid flow state parameters of each cross section of the catheter through a microfluidic sensor array built in the catheter, identify a local turbulence region according to the fluid flow state parameters and the multi-dimensional parameter feature matrix, calculate a turbulence development trend when the turbulence intensity of the local turbulence region exceeds a warning intensity threshold, eliminate pressure fluctuation according to the turbulence development trend and the frequency domain characteristics of pressure fluctuation, and control the pressure in the catheter to be maintained within a safe pressure range.

[0042] In a third aspect of the embodiment of the present application,

[0043] An electronic device is provided, comprising:

[0044] a processor;

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

[0046] wherein the processor is configured to invoke the instructions stored by the memory to perform the method as described above.

[0047] A fourth aspect of the embodiments of the present application,

[0048] A computer-readable storage medium is provided, which stores computer program instructions, and the computer program instructions are executed by a processor to implement the method as described above.

[0049] The beneficial effects of the present application are as follows:

[0050] The present application realizes accurate real-time monitoring of the pressure in the catheter during the breast ductoscopy process through multi-dimensional parameter feature matrix analysis and adaptive filtering technology, can timely find abnormal pressure areas, and prevent tissue damage risk caused by excessive pressure.

[0051] Through microfluidic sensor array and fluid flow state parameter analysis, the present application can accurately identify local turbulent flow areas, and predict the development trend of turbulent flow through wavelet transform analysis method, realize active early warning of the pressure in the catheter, and improve the safety of breast ductoscopy.

[0052] The present application can effectively eliminate the pressure fluctuation in the catheter by adopting active pulsation compensation technology to control the pressure in the catheter, so that the pressure is always maintained within a safe range, significantly reduces the risk of tissue damage during breast ductoscopy, and improves the comfort and diagnostic accuracy of the examination. BRIEF DESCRIPTION OF DRAWINGS

[0053] Figure 1 It is a flowchart of the method for real-time monitoring and control of the pressure in the catheter during breast ductoscopy of the embodiments of the present application;

[0054] Figure 2 It is a schematic diagram for comparison of adaptive filtering performance. DETAILED DESCRIPTION

[0055] In order to make the purpose, technical scheme and advantages of the embodiments of the present application clearer, the technical scheme in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, not all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor are within the scope of protection of the present application.

[0056] The technical solutions of the present application will be described in detail below with specific examples. The following specific examples can be combined with each other, and the same or similar concepts or processes may not be described in detail in some examples.

[0057] Figure 1 The flowchart of the duct pressure real-time monitoring and control method in the breast ductoscopy examination of the embodiment of the present application is shown in FIG. 1, which comprises the following steps: Figure 1

[0058] Collecting real-time temperature data, conductivity data and image data in the duct during the breast ductoscopy examination and forming a multi-dimensional parameter feature matrix;

[0059] Performing adaptive filtering analysis on the multi-dimensional parameter feature matrix, extracting the characteristic parameters of the sub-millimeter fluid particles in the fluid, calculating the spatial distribution density value and the motion velocity vector of the sub-millimeter fluid particles according to the characteristic parameters, calculating the spatial distribution characteristics of the duct pressure according to the spatial distribution density value, determining the positions of the pressure peak point and the pressure valley point, and identifying the pressure gradient mutation region; tracking the pressure wave propagation path according to the motion velocity vector, identifying the reflection point and the attenuation point of the pressure wave on the propagation path, and obtaining the pressure wave propagation characteristic distribution map;

[0060] Combining the pressure gradient mutation region and the pressure wave propagation characteristic distribution map, the pressure in the duct is monitored in real time;

[0061] Monitoring the fluid flow state parameters of each cross section of the duct through the microfluidic sensor array built in the duct, identifying the local turbulent flow region according to the fluid flow state parameters and the multi-dimensional parameter feature matrix, calculating the turbulent flow development trend when the turbulent flow intensity of the local turbulent flow region exceeds the early warning intensity threshold, eliminating the pressure fluctuation according to the turbulent flow development trend and the frequency domain characteristics of the pressure fluctuation, and controlling the pressure in the duct to maintain within the safe pressure range.

[0062] In an optional implementation, the adaptive filtering analysis on the multi-dimensional parameter feature matrix, the extraction of the characteristic parameters of the sub-millimeter fluid particles in the fluid, and the calculation of the spatial distribution density value and the motion velocity vector of the sub-millimeter fluid particles according to the characteristic parameters comprise:

[0063] Calculating the autocorrelation matrix of the multi-dimensional parameter feature matrix through adaptive Wiener filtering to obtain the filtered feature data;

[0064] In the filtered feature data, the edge points are detected by using a multi-scale operator, the complete contour is obtained by performing morphological connection processing on the edge points, the equivalent diameter and the circularity coefficient are calculated according to the complete contour, the gray scale gradient value of the contour region is calculated, and the equivalent diameter, the circularity coefficient and the gray scale gradient value are combined to form the characteristic parameters of the sub-millimeter fluid particles; ​

[0065] determining an initial bandwidth according to the feature parameters, iteratively optimizing the initial bandwidth through cross-validation until a density estimation error is less than a distribution fitting threshold to obtain an optimized bandwidth parameter, and calculating a spatial distribution density value of the sub-millimeter fluid particles using the optimized bandwidth parameter;

[0066] extracting images of adjacent frames from the filtered feature data, calculating spatial and temporal gradients of the adjacent frame images in combination with the feature parameters and constructing a velocity field smoothness constraint, constructing an optimization objective function using the spatial and temporal gradients and the velocity field smoothness constraint, and determining a motion velocity vector of the sub-millimeter fluid particles by iteratively optimizing the optimization objective function through a gradient descent method until a value of the optimization objective function is less than a velocity field convergence threshold.

[0067] An autocorrelation matrix R of the feature matrix is calculated, which has a size of N x N, where N is the dimension of the feature vector. In actual applications, taking a 5-dimensional feature vector as an example, the autocorrelation matrix is a 5 x 5 matrix. The trace value Tr(R) of the autocorrelation matrix is calculated, which reflects the overall energy of the data. The filter step size μ is dynamically adjusted according to the trace value, that is, when the trace value is large (such as greater than 100), the step size is set to 0.01; when the trace value is small (such as less than 50), the step size is set to 0.05, to ensure the stability and convergence speed of the algorithm. The filter coefficient w is iteratively updated, and the iteration formula is based on the least mean square error criterion. The coefficient is updated each time the iteration is performed, until the variation of the coefficient between two adjacent iterations is less than a preset threshold (such as 0.001) or the maximum number of iterations (such as 50 times) is reached. Finally, the filtered feature data is obtained, and the noise of the data is significantly reduced, and the edge features of the sub-millimeter bubbles and colloidal clusters are retained.

[0068] Multi-scale processing is achieved by setting different Gaussian kernel standard deviations σ. Generally, three scale values are set: σ = 0.8, 1.2 and 1.6. The image gradient is calculated for each scale, and non-maximum suppression is performed. The double-threshold method is used to determine the edge points, the low threshold is set to 40% of the high threshold, and the high threshold is automatically determined by the OTSU method, which is usually the 85% quantile point of the gradient histogram. The edge points detected at different scales are fused to obtain a complete set of edge points. Morphological connection processing is applied to the edge points, including thinning, closing and connection operations, and the connection distance threshold is set to 5 pixels to form a complete target contour. The geometric feature parameters are calculated according to the contour: area S (number of pixels), perimeter L (contour length), equivalent diameter and circularity coefficient C = 4πS / L 2The average value G and standard deviation σG of the gray level gradient within the contour region are calculated simultaneously. These parameters form the feature parameter set of sub-millimeter bubbles and colloidal clusters. In practical cases, typical sub-millimeter bubbles have an equivalent diameter ranging from 0.1 mm to 0.9 mm and a circularity coefficient greater than 0.85, while colloidal clusters have a similar equivalent diameter but a circularity coefficient usually lower than 0.75, and the gray level gradient value of bubbles is usually higher than that of clusters.

[0069] The initial Gaussian kernel function bandwidth h0 is selected based on the Silverman criterion, h0 = 0.9 × min(σ, IQR / 1.34) × n -1 / 5 , where σ is the sample standard deviation, IQR is the interquartile range, and n is the sample size. The initial bandwidth is optimized by k-fold cross-validation (k = 5), dividing the data into 5 parts, using 4 parts to estimate the density function each time, and using the remaining 1 part to calculate the estimation error. Within the bandwidth range [0.5h0, 1.5h0], select the bandwidth value that minimizes the error with a step size of 0.1h0. Iterate the optimization until the density estimation error is less than a preset threshold (e.g., 0.01). Use the optimized bandwidth h* to calculate the density value ρ(x, y) at each point (x, y) in the two-dimensional space. In practical applications, the bubble distribution density value in the fluid generally ranges from 0 to 25 per cm 2 , and high-density areas correspond to fluid dynamics instability areas.

[0070] From the filtered feature data, adjacent frame images I(x, y, t) and I(x, y, t + Δt) are extracted, typically Δt = 1 / 30 seconds (corresponding to a 30 fps acquisition rate). Combined with the position and shape information in the feature parameters, the spatial gradients Ix, Iy and the temporal gradient It of the adjacent frame images are calculated. The optical flow constraint equation is constructed, and a velocity field smoothing constraint is added to form the optimization objective function J. The velocity field is iteratively solved by the gradient descent method, with an update step size of 0.05, and the iteration is stopped when the objective function changes less than a preset threshold (e.g., 0.0005) or reaches the maximum iteration number (e.g., 100 times). The final velocity vector (u, v) of each bubble and cluster is obtained, representing its movement speed in the x and y directions. In practical case analysis, the typical velocity of bubbles in the fluid ranges from 0 to 15 cm / s, and the velocity field distribution reflects the flow state and turbulent characteristics inside the fluid. The velocity vector diagram can visually display the vortex structure and shear region in the fluid, providing important basis for fluid dynamics analysis.

[0071] Figure 2For the adaptive filtering performance comparison diagram, the diagram shows the performance comparison of the adaptive filtering algorithm proposed in the application and four kinds of classical filtering methods in the process of 50 iterations. The vertical axis represents the filtering error (unit: dB), and the smaller the negative value, the better the noise suppression effect. From the figure, it can be clearly observed that the adaptive Wiener filtering algorithm of the application (black solid line) shows the best filtering performance, which quickly converges from the initial -9.2 dB to -30.3 dB, which is 12.9 dB higher than the traditional Wiener filter of -17.4 dB, 6.8 dB higher than the RLS adaptive filter of -23.5 dB, 14.3 dB higher than the LMS adaptive filter of -16.0 dB, and 8.7 dB higher than the Kalman filter of -21.6 dB. It is particularly worth noting that the application shows faster convergence speed in the first 20 iterations, which is due to the dynamic step size adjustment mechanism based on the trace value of the autocorrelation matrix, which can adaptively optimize the filter parameters according to the data characteristics. In the later stage of iteration, the application can still continuously improve the filtering effect and keep stable convergence, while other methods reach the performance bottleneck at an early stage, which fully verifies the technical advantages and practical value of the application in submillimeter bubble and colloid cluster feature extraction.

[0072] In an optional embodiment, the spatial distribution feature of the pressure in the catheter is calculated according to the spatial distribution density value, the positions of the pressure peak point and the pressure valley point are determined, and the pressure gradient mutation region is identified, including:

[0073] The initial pressure distribution is calculated based on the spatial distribution density value in the catheter;

[0074] The Laplacian of the density is calculated, the density gradient correction term is obtained based on the Laplacian and the density gradient coefficient, and the density gradient correction term is superimposed on the initial pressure distribution to obtain the corrected pressure distribution;

[0075] The corrected pressure distribution is subjected to multi-scale analysis, and local extreme points are extracted at each scale; the convergence and divergence of the pressure gradient vector field corresponding to each local extreme point are analyzed, the local extreme point with the convergent gradient vector field is marked as the pressure peak point, and the local extreme point with the divergent gradient vector field is marked as the pressure valley point;

[0076] According to the distribution positions of the pressure peak point and the pressure valley point, the gradient amplitude of the corrected pressure distribution at each position is calculated, the gradient direction angle between adjacent positions is calculated to obtain the direction change rate, the gradient amplitude is divided by the average value of the gradient amplitude in the local region to obtain the amplitude ratio, and the gradient mutation index is calculated according to the amplitude ratio and the direction change rate. When the gradient mutation index is greater than the mutation judgment threshold, the region between the adjacent pressure peak point and the pressure valley point is marked as the pressure gradient mutation region.

[0077] In the method for calculating the pressure distribution characteristics in the conduit, the spatial distribution density values in the conduit need to be obtained, which can be collected by sensors at different positions in the conduit or simulated by computational fluid dynamics. For example, in a conduit with a length of 200 millimeters, density values can be measured or calculated at every 10 millimeter interval, forming a density distribution data set containing 21 measurement points.

[0078] When calculating the initial pressure distribution based on the obtained spatial distribution density values in the conduit, the state equation can be used for conversion. For liquids, a relationship similar to p = pg h can be applied, where p is the density, g is the acceleration of gravity, and h is the height of the liquid column; for gases, a relationship similar to p = pRT / M can be applied, where R is the gas constant, T is the temperature, and M is the molecular weight. Taking the gas as an example, if the density at a certain measurement point is 1.2 kilograms per cubic meter, the temperature is 293 kelvins, and the molecular weight is 29, the initial pressure at this point can be calculated to be approximately 101325 pascals.

[0079] In discrete space, it can be realized by finite difference method. For one-dimensional case, three-point difference formula can be used; for three-dimensional case, second-order derivatives in x, y, z directions can be calculated respectively, and then summed. For example, in one-dimensional case, the Laplacian of point i can be calculated by (pi+1 - 2pi + pi-1) / Ax 2 , where Ax is the spatial sampling interval. If the density of point i is 1.2 kilograms per cubic meter, the density of point i+1 is 1.25 kilograms per cubic meter, the density of point i-1 is 1.18 kilograms per cubic meter, and Ax is 10 millimeters, the Laplacian value is (1.25 - 2 x 1.2 + 1.18) / (0.01) 2 = 300 kilograms per cubic meter per meter 2 .

[0080] When multiplying the Laplacian with the density gradient coefficient to obtain the density gradient correction term, an appropriate density gradient coefficient needs to be selected. The coefficient can be determined according to the fluid characteristics and conduit material, and the typical value is between 0.001 and 0.1. If the density gradient coefficient is selected as 0.01, the density gradient correction term obtained by the Laplacian in the above calculation is 300 x 0.01 = 3 pascals. Adding this correction term to the initial pressure distribution, for example, the initial pressure at this point is 101325 pascals, then the corrected pressure is 101325 + 3 = 101328 pascals.

[0081] The standard deviation is set to 1 mm, 2 mm, and 4 mm. At each scale, the local extreme points are extracted by comparing the pressure value of each point with its neighborhood points. If the pressure value of a point is greater than the pressure values of all its neighborhood points, the point is a local maximum point; if it is less than the pressure values of all its neighborhood points, the point is a local minimum point.

[0082] For a point (x, y, z) in three-dimensional space, its pressure gradient vector contains three components, x, y, and z, which represent the rate of change of pressure in these three directions. The convergence and divergence of the gradient vector can be determined by calculating its divergence. If the divergence is negative, it means that the gradient vector field converges at this point, and the point is a pressure peak point; if the divergence is positive, it means that the gradient vector field diverges at this point, and the point is a pressure valley point. In practical applications, a threshold value can be set, for example, ±0.5 Pa / mm 2 , and only when the absolute value of the divergence is greater than the threshold value is it considered an effective peak point or valley point.

[0083] For a one-dimensional case, the gradient amplitude of point i can be calculated by |(pi+1 - pi-1) / (2Δx)|. For a three-dimensional case, the gradient amplitude is the square root of the sum of the squares of the gradient components in each direction. For example, if the gradient of a point in the x direction is 0.5 Pa / mm, in the y direction is 0.3 Pa / mm, and in the z direction is 0.4 Pa / mm, then the gradient amplitude of the point is Pa / mm.

[0084] The angle θ between two gradient vectors v1 and v2 can be calculated by cos(θ)=(v1·v2) / (|v1|·|v2|), and then the angle value is converted to radians or degrees as the direction change rate. For example, if the gradient vectors of two adjacent points are (0.5, 0.3, 0.4) and (0.2, 0.6, 0.1) respectively, their angle is about 45 degrees, and the corresponding direction change rate is 0.785 radians.

[0085] A spherical or cubic region centered on the current point can usually be selected, such as a spherical region with a radius of 5 mm. The average of the gradient amplitudes of all points in the region is calculated, and then the gradient amplitude of the current point is divided by this average to obtain the amplitude ratio. For example, if the gradient amplitude of the current point is 0.71 Pa / mm and the average of the gradient amplitudes in the local region is 0.5 Pa / mm, then the amplitude ratio is 0.71 / 0.5=1.42.

[0086] When calculating the gradient mutation index with the weighted sum of the amplitude ratio and the direction change rate, the weight coefficients of the two need to be determined. These weights can be adjusted according to specific application scenarios, for example, setting the weight of the amplitude ratio to 0.7 and the weight of the direction change rate to 0.3. If the amplitude ratio is 1.42 and the direction change rate is 0.785 rad, then the gradient mutation index is 0.7 x 1.42 + 0.3 x 0.785 ≈ 1.229.

[0087] When the gradient mutation index is greater than the mutation determination threshold, the region between the adjacent pressure peak point and the pressure valley point is marked as a pressure gradient mutation region. The mutation determination threshold can be determined according to historical data and expert experience, for example, set to 1.0. In the above example, the gradient mutation index is 1.229, which is greater than the threshold 1.0, so the region is marked as a pressure gradient mutation region. These mutation regions usually correspond to positions where the fluid dynamics characteristics in the catheter change significantly, such as changes in cross-sectional area, bends, or bifurcations, which are of great significance to catheter design and fluid control.

[0088] In an alternative embodiment, the pressure wave propagation path is tracked according to the motion velocity vector, the reflection points and the attenuation points of the pressure wave on the propagation path are identified, and the pressure wave propagation characteristic distribution map is obtained, comprising:

[0089] The velocity component of the discrete point in the corrected pressure distribution is extracted, the pressure wave acceleration is calculated according to the velocity components of adjacent discrete points, and the pressure wave propagation path is obtained by integrating the pressure wave acceleration along time;

[0090] The inner product of the velocity direction at each time and the velocity direction at the next time along the pressure wave propagation path is calculated, and the propagation direction change angle is obtained according to the inner product and the velocity amplitude. When the propagation direction change angle is greater than the direction mutation threshold, the corresponding coordinate point is marked as a reflection point;

[0091] The pressure difference between adjacent coordinate points on the pressure wave propagation path is divided by the distance between the two points to obtain the pressure gradient along the path, and the attenuation coefficient is obtained by dividing the pressure gradient along the path by the pressure value. When the attenuation coefficient is greater than the attenuation feature threshold, the corresponding coordinate point is marked as an attenuation point;

[0092] The distance from each reflection point to each attenuation point is calculated, the reflection points and attenuation points with a distance less than the feature correlation threshold are paired to form a feature point group, and the spatial distribution density of the feature point group is calculated to obtain the propagation feature strength;

[0093] The pressure wave propagation path is drawn, the propagation feature strength is labeled, the reflection point position and reflection direction are marked, and the attenuation point position and attenuation coefficient are marked to generate a pressure wave propagation characteristic distribution map.

[0094] For the acquired pressure distribution data, first the pressure data of each discrete point is extracted and corrected. The corrected pressure distribution data contains the position coordinates of each discrete point in time series and the corresponding pressure value. From these data, the velocity components of each discrete point in x and y directions are extracted, for example, at a certain time t1, the x direction velocity component of discrete point P1 is 3.5 m / s, and the y direction velocity component is 2.1 m / s. According to the velocity components of adjacent time, the pressure wave acceleration is calculated, for example, for adjacent time t1 and t2, the acceleration ax is calculated as (vx2-vx1) / (t2-t1), and ay is calculated in the same way. Specifically, if the x direction velocity of P1 point at t1 time is 3.5 m / s, and the x direction velocity at t2 time is 4.2 m / s, the time interval is 0.01 seconds, then the x direction acceleration is 70 m / s 2 . The pressure wave propagation path is obtained by integrating the pressure wave acceleration along time, that is, the complete propagation path is constructed by accumulating the displacement increments at each time.

[0095] The change of propagation direction along the pressure wave propagation path is calculated. For each time on the propagation path, the inner product of the velocity direction at the current time and the velocity direction at the next time is calculated. For example, the velocity vector at t1 time is (3.5, 2.1), and the velocity vector at t2 time is (3.8, 1.9), then the inner product is calculated as 3.5x3.8+2.1x1.9=17.29. The inner product is divided by the product of the module lengths of the two velocity vectors to obtain the direction cosine. In this example, the module length of the velocity vector at t1 time is 4.07, and the module length at t2 time is 4.25, and the direction cosine is 17.29 / (4.07x4.25)=0.9996. The direction cosine is converted to radian value to obtain the propagation direction change angle arccos(0.9996)=0.0283 radian, which is about 1.62 degrees. When the direction change angle is greater than the preset direction mutation threshold (for example, set to 10 degrees or 0.1745 radian), the corresponding coordinate point is marked as a reflection point. In this example, the point is not marked as a reflection point.

[0096] For the pressure attenuation analysis on the propagation path, the pressure difference between adjacent coordinate points is divided by the distance between the two points to obtain the along-path pressure gradient. For example, the pressure of P1 point at t1 time is 120 kPa, and the pressure of P2 point at t2 time is 115 kPa, the distance between the two points is 0.05 m, then the along-path pressure gradient is (120-115) / 0.05=100 kPa / m. The attenuation coefficient is obtained by dividing the along-path pressure gradient by the pressure value, which is 100 / 120=0.833 m -1 . When the attenuation coefficient is greater than the preset attenuation feature threshold (for example, set to 0.5 m -1 ), the corresponding coordinate point is marked as an attenuation point. In this example, the point is marked as an attenuation point.

[0097] For the marked reflection points and attenuation points, feature correlation analysis is performed. The distance of each reflection point to each attenuation point is calculated, and when the distance is less than a feature correlation threshold (for example, set to 0.2 m), the pair of reflection point and attenuation point is paired to form a feature point group. For example, the distance between reflection point R1(2.5, 1.8) and attenuation point A1(2.6, 1.9) is less than the feature correlation threshold 0.2 m, so the two points form a feature point group. The distribution density of all feature point groups in space is counted to obtain the propagation feature intensity. For example, there are 5 feature point groups in an area of 1 square meter, so the propagation feature intensity of the area is 5 / m2.

[0098] Finally, the pressure wave propagation path is plotted in a two-dimensional coordinate system, and the propagation feature intensity is represented by color depth or line thickness. Stronger features are represented by darker colors or thicker lines. In the figure, the positions of the reflection points are marked, and the reflection directions are indicated by arrows. At the same time, the positions of the attenuation points are marked, and the attenuation coefficients are labeled with numerical values. For example, reflection point R1 is marked with a red dot and the pre-reflection direction (0.8, 0.6) and the post-reflection direction (0.6, -0.8) are indicated by arrows; attenuation point A1 is marked with a blue triangle and the attenuation coefficient 0.833 m -1 is labeled. In this way, the pressure wave propagation characteristic distribution map visually displays the reflection and attenuation characteristics of the pressure wave during propagation, which helps to analyze the propagation law of the pressure wave and predict the propagation behavior.

[0099] The above method realizes comprehensive analysis and visual display of the propagation characteristics of the pressure wave by accurately tracking the propagation path of the pressure wave and identifying feature points, providing an effective technical means for pressure wave propagation research.

[0100] In an alternative embodiment, the fluid flow state parameters of each cross section of the catheter are monitored by a microfluid sensor array built into the catheter, and the local turbulent flow region is identified based on the fluid flow state parameters and the multi-dimensional parameter feature matrix. When the turbulent intensity of the local turbulent flow region is detected to exceed the warning intensity threshold, the turbulent development trend is calculated, including:

[0101] The pressure fluctuation values and flow velocity distribution values of each cross section are collected and mapped with the multi-dimensional parameter feature matrix to obtain the strain rate tensor and vorticity tensor. The second invariant of the vorticity tensor and the second invariant of the rate of change tensor are calculated and subtracted to obtain the turbulent discriminant factor. The local turbulent flow region is marked based on the turbulent discriminant factor.

[0102] The flow velocity fluctuation component in the local turbulent flow region is extracted, and the ratio of the root mean square value of the flow velocity fluctuation component to the average flow velocity is calculated to obtain the turbulent intensity.

[0103] When the turbulence intensity exceeds the early warning intensity threshold, the flow velocity fluctuation signal of the local turbulence region is extracted and wavelet transformed, the wavelet coefficient matrix is integrated at different scales, the square of each element in the wavelet coefficient matrix is calculated to obtain the local wave energy, the scale energy is obtained by summing the local wave energy at the same scale, and the energy distribution interval is extracted based on the scale energy;

[0104] In the energy distribution interval, the difference between the turbulence energies of adjacent time instants is divided by the time interval to obtain the energy change rate, the difference between the characteristic lengths of adjacent scales is divided by the time interval to obtain the scale change rate, the turbulence development index is obtained based on the energy change rate and the scale change rate, and the turbulence development trend is judged according to the time sequence change characteristics of the turbulence development index.

[0105] The in-catheter microfluid sensor array is arranged at a plurality of key cross-section positions of the catheter, a plurality of sensor nodes are arranged at each cross-section to form a ring-shaped distribution, and are used to collect pressure fluctuation values, flow velocity distribution values and shear stress values of the fluid in the catheter. The sensor node collection frequency is set to 2000 Hz, which meets the capture demand of the high-frequency fluctuation characteristics of turbulence. The sensor unit adopts a structure combining a hot film type flow velocity sensor and a piezoresistive pressure sensor, the size of a single sensor is 100 microns x 100 microns, and the thickness is 5 microns, which guarantees the measurement accuracy and does not interfere with the fluid flow in the catheter.

[0106] The multi-dimensional parameter feature matrix is a reference database constructed based on a large amount of experimental data and numerical simulation results, and contains characteristic parameters in various flow states. The mapping operation adopts a tensor inner product method to calculate a strain rate tensor and a vorticity tensor. The strain rate tensor describes the deformation rate of the fluid, and the vorticity tensor represents the rotation characteristics of the fluid. For a certain actual measurement, the flow velocity distribution value detected at the third cross-section of the catheter presents obvious non-uniformity in the radial direction, the central flow velocity is 1.2 m / s, and the flow velocity close to the wall is only 0.3 m / s, which indicates that there is local turbulence.

[0107] The second invariant of the vorticity tensor is calculated, which represents the vortex intensity; the second invariant of the strain rate tensor is calculated, which represents the deformation intensity of the fluid. The second invariant of the vorticity tensor is subtracted from the second invariant of the strain rate tensor to obtain a turbulence discrimination factor Q. When the Q value is greater than zero, it indicates that the vortex effect exceeds the deformation effect in the region, which is a turbulence region. In actual application, when the Q value is greater than a set threshold value 0.15, the region is marked as a local turbulence region. In this example, the ring-shaped region with a radial position of 0.7R to 0.85R (R is the radius of the catheter) at the third cross-section of the catheter has a Q value of 0.23, which is identified as a local turbulence region.

[0108] The flow velocity fluctuation component is extracted from the local turbulent region, i.e. the actual flow velocity minus the average flow velocity in a period of time. The root mean square value of the flow velocity fluctuation component is calculated, and the ratio of the turbulent intensity is obtained. The turbulent intensity represents the degree of fluid turbulence, and the larger the value, the more intense the turbulence. The early warning intensity threshold is set to 0.15, and when the turbulent intensity exceeds the threshold, it is determined that the turbulence has reached a level that requires early warning. In this test, the turbulent intensity of the local turbulent region is calculated to be 0.18, which exceeds the early warning intensity threshold, triggering the turbulence development trend analysis.

[0109] The flow velocity fluctuation signal of the local turbulent region is extracted, and Morlet wavelet transform is performed on the signal. Morlet wavelet has good time-frequency localization characteristics and is suitable for analyzing non-stationary turbulent signals. The wavelet transform scale range is set to 1 to 64, with a total of 32 scale points, corresponding to a physical scale from 0.5 mm to 25 mm. The wavelet coefficients are calculated at each scale to form a wavelet coefficient matrix. The number of rows in the matrix is the number of time sampling points, and the number of columns is the number of scale points.

[0110] The square of each element in the wavelet coefficient matrix is calculated to obtain the local wave energy distribution. The local wave energy is summed at the same scale to obtain the energy value at each scale. The energy-scale distribution graph is drawn, with the horizontal axis representing the scale and the vertical axis representing the energy. The energy distribution interval is extracted from the distribution graph, and in this example, the energy is mainly concentrated in the scale interval of 6 to 18, corresponding to a physical scale of about 3 mm to 9 mm, indicating that the turbulent structure is mainly concentrated in the small scale range.

[0111] In the energy distribution interval, the difference between the turbulent energy of adjacent time points and the time interval is calculated to obtain the energy change rate. The time interval is set to 0.01 seconds. At the same time, the difference between the characteristic length of adjacent scales and the time interval is calculated to obtain the scale change rate. The characteristic length corresponds to the physical length of the scale. The energy change rate and the scale change rate are used to calculate the turbulence development index by weighted summation, with weights of 0.7 and 0.3 respectively. In this example, the turbulence development index of the initial 10 time points is 0.12, 0.15, 0.19, 0.24, 0.28, 0.33, 0.38, 0.42, 0.47, 0.52, showing a clear upward trend.

[0112] When the index continues to grow and the growth rate exceeds 0.03 / s, it is determined that the turbulence has an enhancing trend; when the index fluctuates but is stable within a range (with a fluctuation amplitude not exceeding ±0.05), it is determined that the turbulence has a stable trend; when the index continues to decline and the decline rate exceeds 0.03 / s, it is determined that the turbulence has a weakening trend. According to the above judgment standard, in this example, the turbulence presents a clear enhancing trend, with an average growth rate of 0.044 / s, generating a turbulence enhancement early warning signal, prompting the operator to take appropriate measures to reduce the flow velocity or adjust the catheter position to prevent the turbulence from further developing and causing potential risks.

[0113] In an alternative embodiment, eliminating pressure fluctuation, controlling the pressure in the conduit to maintain within a safe pressure range according to the frequency domain characteristics of the pressure fluctuation and the turbulence development trend comprises:

[0114] extracting the pressure fluctuation signal in the conduit, decomposing the pressure fluctuation signal in the time domain into periodic fluctuation components, obtaining the fluctuation main frequency and its amplitude according to the Fourier series decomposition result of the periodic fluctuation components;

[0115] setting a conduit pressure target value, taking the difference between the real-time detected conduit pressure value and the conduit pressure target value as a pressure deviation;

[0116] calculating a proportional coefficient, an integral time constant and a differential time constant according to the turbulence development trend, the fluctuation main frequency and its amplitude, multiplying the proportional coefficient with the pressure deviation to obtain a proportional control term, multiplying the integral time constant with the pressure deviation integral value to obtain an integral control term, multiplying the differential time constant with the pressure deviation change rate to obtain a differential control term, adding the proportional control term, the integral control term and the differential control term to obtain a compensation control signal;

[0117] inputting the compensation control signal into a piezoelectric ceramic diaphragm pump to generate a compensation pressure, real-time detecting the compensated conduit pressure value, when the compensated conduit pressure value exceeds the upper limit value or the lower limit value of the safe pressure range, calculating the change rate of the pressure deviation, adjusting the compensation pressure according to the change rate of the pressure deviation to maintain the pressure value in the conduit within the safe pressure range.

[0118] By installing a micro pressure sensor in the conduit, the pressure change in the conduit is measured in real time. The micro pressure sensor has a range of -50~150kPa, an accuracy of 0.01kPa and a sampling frequency of 1000Hz. The collected pressure fluctuation signal is represented by a time series P(t), where t represents the time variable.

[0119] The collected pressure fluctuation signal P(t) is decomposed in the time domain, and the periodic fluctuation component Pp(t) is extracted by using wavelet analysis method. For example, when the pressure fluctuation signal P(t) in the conduit has values of 92.5kPa, 93.1kPa, 94.0kPa, 93.5kPa, 92.8kPa, 92.2kPa, 91.8kPa, 92.4kPa, 93.2kPa and 93.9kPa within 10 seconds, the periodic fluctuation component Pp(t) obtained after wavelet analysis is 0.2kPa, 0.8kPa, 1.7kPa, 1.2kPa, 0.5kPa, -0.1kPa, -0.5kPa, 0.1kPa, 0.9kPa and 1.6kPa.

[0120] The Pp(t) is analyzed by using the fast Fourier transform algorithm to obtain the frequency components and the corresponding amplitudes. In the above example, the Fourier analysis can obtain the main frequency f of the fluctuation as 0.2 Hz, and the corresponding amplitude A as 1.2 kPa.

[0121] The conduit pressure target value Ptarget is set, which is set as 92.3 kPa in the embodiment. The pressure value Pcurrent in the conduit is detected in real time, and the pressure deviation e = Pcurrent - Ptarget is calculated. For example, when the pressure value in the conduit at a certain time is 94.0 kPa, the pressure deviation e = 94.0 kPa - 92.3 kPa = 1.7 kPa.

[0122] The PID control parameters are calculated according to the turbulence development trend, the main frequency of the fluctuation and the amplitude. Specifically, the proportional coefficient Kp is inversely proportional to the fluctuation amplitude A and proportional to the turbulence intensity; the integral time constant Ti is inversely proportional to the main frequency f of the fluctuation; and the differential time constant Td is proportional to the rate of change of the main frequency f of the fluctuation. In the embodiment, when the turbulence intensity is 0.3, the main frequency f of the fluctuation is 0.2 Hz, and the amplitude A is 1.2 kPa, the proportional coefficient Kp = 0.25, the integral time constant Ti = 10 s, and the differential time constant Td = 0.05 s are calculated.

[0123] The PID control output values are calculated. The proportional control term Pout = Kp × e, the integral control term Iout = e × dt / Ti, and the differential control term Dout = Td × de / dt, where dt is the sampling time interval, and de / dt is the rate of change of the pressure deviation. In the above example, when e = 1.7 kPa, the previous e = 1.2 kPa, and dt = 0.01 s, Pout = 0.25 × 1.7 kPa = 0.425 kPa, Iout = 1.7 kPa × 0.01 s / 10 s = 0.0017 kPa, and Dout = 0.05 s × (1.7 kPa - 1.2 kPa) / 0.01 s = 2.5 kPa are calculated.

[0124] The proportional control term, the integral control term and the differential control term are added to obtain the compensation control signal u = Pout + Iout + Dout. In the above example, u = 0.425 kPa + 0.0017 kPa + 2.5 kPa = 2.9267 kPa.

[0125] The compensation control signal u is input to the piezoelectric diaphragm pump to generate a compensation pressure. The output pressure range of the piezoelectric diaphragm pump is 0-50 kPa, and the response time is less than 5 ms. The compensation pressure Pcomp = k × u, where k is the pressure conversion coefficient of the piezoelectric diaphragm pump, and in this embodiment k = 1. In the above example, the compensation pressure Pcomp = 1 × 2.9267 kPa = 2.9267 kPa.

[0126] The compensated catheter pressure value Pnew = Pcurrent - Pcomp is detected in real time by the micro pressure sensor. In the above example, Pnew = 94.0 kPa - 2.9267 kPa = 91.0733 kPa.

[0127] The upper limit value Pupper = 95.0 kPa and the lower limit value Plower = 90.0 kPa of the safety pressure range are set. It is determined whether Pnew is within the safety pressure range, i.e., whether Plower ≤ Pnew ≤ Pupper is satisfied. In the above example, 91.0733 kPa is between 90.0 kPa and 95.0 kPa, satisfying the safety pressure range requirement.

[0128] When the compensated catheter pressure value exceeds the safety pressure range, the rate of change of the pressure deviation de / dt = (e - eprev) / dt is calculated, where eprev is the pressure deviation at the previous time. According to the size and sign of de / dt, the compensation pressure value is dynamically adjusted. For example, when Pnew = 89.5 kPa < Plower, de / dt = (e - eprev) / dt = (1.7 kPa - 1.2 kPa) / 0.01 s = 50 kPa / s is calculated. Since de / dt > 0 and Pnew < Plower, it indicates that the pressure is decreasing too fast, and the compensation pressure needs to be reduced. At this time, the compensation pressure is adjusted to Pcomp_adj = Pcomp × (1 - |de / dt| / 100) = 2.9267 kPa × (1 - 50 / 100) = 1.46335 kPa. The adjusted catheter pressure value is Pnew_adj = Pcurrent - Pcomp_adj = 94.0 kPa - 1.46335 kPa = 92.53665 kPa, satisfying the safety pressure range requirement.

[0129] Through the above method, effective control of the pressure in the catheter is achieved, which maintains it within the safety pressure range, avoiding safety risks caused by excessively high or low pressure.

[0130] The duct pressure real-time monitoring and control system in the breast ductoscope examination of the embodiment of the application comprises:

[0131] The first unit is used for collecting real-time temperature data, conductivity data and image data in the duct during the breast ductoscope examination and forming a multi-dimensional parameter feature matrix;

[0132] The second unit is used for performing adaptive filtering analysis on the multi-dimensional parameter feature matrix, extracting the characteristic parameters of the sub-millimeter fluid particles in the fluid, calculating the spatial distribution density value and the motion velocity vector of the sub-millimeter fluid particles according to the characteristic parameters, calculating the spatial distribution characteristics of the pressure in the duct according to the spatial distribution density value, determining the positions of the pressure peak point and the pressure valley point, identifying the pressure gradient mutation region, tracking the pressure wave propagation path according to the motion velocity vector, identifying the reflection point and the attenuation point of the pressure wave on the propagation path, and obtaining a pressure wave propagation characteristic distribution map;

[0133] The third unit is used for combining the pressure gradient mutation region and the pressure wave propagation characteristic distribution map to monitor the pressure in the duct in real time.

[0134] The fourth unit is used for monitoring the fluid flow state parameters of each section of the duct through the micro-fluid sensor array built in the duct, identifying the local turbulent flow region according to the fluid flow state parameters and the multi-dimensional parameter feature matrix, calculating the turbulent flow development trend when the turbulent flow intensity of the local turbulent flow region exceeds the early warning intensity threshold, eliminating the pressure fluctuation according to the turbulent flow development trend and the frequency domain characteristics of the pressure fluctuation, and controlling the pressure in the duct to be maintained within a safe pressure range.

[0135] In a third aspect, an electronic device is provided, comprising:

[0136] a processor;

[0137] a memory storing processor-executable instructions;

[0138] The processor is configured to invoke the instructions stored in the memory to execute the method described above.

[0139] In a fourth aspect, a computer-readable storage medium is provided, which stores computer program instructions, and the computer program instructions are executed by a processor to implement the method described above.

[0140] The present application can be a method, device, system and / or computer program product. The computer program product can include a computer readable storage medium having computer readable program instructions loaded thereon for performing various aspects of the present application.

[0141] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present application, and are not intended to limit the present application; although the present application has been described in detail with reference to the above embodiments, those skilled in the art should understand that the technical solutions recorded in the above embodiments can be modified, or some or all of the technical features can be replaced by equivalents; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present application.

Claims

1. A method for real-time monitoring and control of pressure within a catheter during a ductoscopy, characterized in that, The application relates to a real-time monitoring method for fluid flow state in a catheter, comprising the following steps: Collecting real-time temperature data, conductivity data and image data in the catheter during a ductoscopy process and forming a multi-dimensional parameter feature matrix; Adaptive filtering analysis is performed on the multi-dimensional parameter feature matrix to extract characteristic parameters of sub-millimeter fluid particles in the fluid, the spatial distribution density value and the motion velocity vector of the sub-millimeter fluid particles are calculated according to the characteristic parameters, the spatial distribution characteristics of the pressure in the catheter are calculated according to the spatial distribution density value, the positions of the pressure peak point and the pressure valley point are determined, and the pressure gradient mutation region is identified; the pressure wave propagation path is tracked according to the motion velocity vector, the reflection point and the attenuation point of the pressure wave on the propagation path are identified, and a pressure wave propagation characteristic distribution map is obtained; The pressure in the catheter is monitored in real time by combining the pressure gradient mutation region and the pressure wave propagation characteristic distribution map; The fluid flow state parameters of each section of the catheter are monitored through a micro-fluid sensor array arranged in the catheter, the local turbulent flow region is identified according to the fluid flow state parameters and the multi-dimensional parameter feature matrix, the turbulent flow development trend is calculated when the turbulent flow intensity of the local turbulent flow region exceeds a warning intensity threshold value; According to the turbulent flow development trend and the frequency domain characteristics of the pressure fluctuation, the pressure fluctuation is eliminated, and the pressure in the catheter is controlled to be maintained within a safe pressure range.

2. The method of claim 1, wherein, The adaptive filtering analysis of the multi-dimensional parameter feature matrix to extract the characteristic parameters of the sub-millimeter fluid particles in the fluid, and the calculation of the spatial distribution density value and the motion velocity vector of the sub-millimeter fluid particles according to the characteristic parameters comprise the following steps: The autocorrelation matrix of the multi-dimensional parameter feature matrix is calculated through adaptive Wiener filtering to obtain filtered feature data; In the filtered feature data, edge points are detected by using a multi-scale operator, the edge points are subjected to morphological connection processing to obtain complete contours, the equivalent diameter and the circularity coefficient are calculated according to the complete contours, the gray gradient value of the contour region is calculated, and the equivalent diameter, the circularity coefficient and the gray gradient value are combined to form the characteristic parameters of the sub-millimeter fluid particles; The initial bandwidth is determined according to the characteristic parameters, the initial bandwidth is subjected to cross-validation iterative optimization until the density estimation error is less than a distribution fitting threshold value, the optimized bandwidth parameter is obtained, and the spatial distribution density value of the sub-millimeter fluid particles is calculated by using the optimized bandwidth parameter; The images of adjacent frames in the filtered feature data are extracted, the spatial gradient and the time gradient of the adjacent frame images are calculated according to the characteristic parameters, and a velocity field smooth constraint is constructed, an optimization objective function is constructed by using the spatial gradient, the time gradient and the velocity field smooth constraint, the optimization objective function value is iteratively optimized through a gradient descent method until the optimization objective function value is less than a velocity field convergence threshold value, and the motion velocity vector of the sub-millimeter fluid particles is determined.

3. The method of claim 1, wherein, The calculation of the spatial distribution characteristics of the pressure in the catheter according to the spatial distribution density value, the determination of the positions of the pressure peak point and the pressure valley point, and the identification of the pressure gradient mutation region comprise the following steps: An initial pressure distribution is calculated based on the spatial distribution density value in the catheter; A Laplacian operator of the density is calculated, a density gradient correction term is obtained based on the Laplacian operator and a density gradient coefficient, and the density gradient correction term is superimposed on the initial pressure distribution to obtain a corrected pressure distribution; Performing multi-scale analysis on the modified pressure distribution, extracting local extreme points at each scale; analyzing the convergence and divergence of the pressure gradient vector field corresponding to each local extreme point, marking the local extreme points with convergent gradient vector field as pressure peak points, and marking the local extreme points with divergent gradient vector field as pressure valley points; According to the distribution positions of the pressure peak points and the pressure valley points, calculating the gradient amplitude of the modified pressure distribution at each position, and calculating the gradient direction angle between adjacent positions to obtain the direction change rate, dividing the gradient amplitude by the average value of the gradient amplitude in the local area to obtain the amplitude ratio, and calculating the gradient mutation index according to the amplitude ratio and the direction change rate, when the gradient mutation index is greater than the mutation judgment threshold, marking the region between adjacent pressure peak points and pressure valley points as the pressure gradient mutation region.

4. The method of claim 1, wherein, According to the motion velocity vector, the pressure wave propagation path is tracked, the reflection points and the attenuation points of the pressure wave on the propagation path are identified, and a pressure wave propagation characteristic distribution map is obtained, including: Extracting the velocity component of the discrete points in the modified pressure distribution, calculating the pressure wave acceleration according to the velocity components of adjacent discrete points, and integrating the pressure wave acceleration along time to obtain the pressure wave propagation path; Calculating the inner product of the velocity direction at each time and the velocity direction at the next time along the pressure wave propagation path, and obtaining the propagation direction change angle according to the inner product and the velocity amplitude, when the propagation direction change angle is greater than the direction mutation threshold, marking the corresponding coordinate point as a reflection point; Dividing the pressure difference between adjacent coordinate points on the pressure wave propagation path by the distance between the two points to obtain the along-path pressure gradient, and dividing the along-path pressure gradient by the pressure value to obtain the attenuation coefficient, when the attenuation coefficient is greater than the attenuation characteristic threshold, marking the corresponding coordinate point as an attenuation point; Calculating the distance from each reflection point to each attenuation point, pairing the reflection points and the attenuation points with distances less than the characteristic correlation threshold to form a characteristic point group, and calculating the spatial distribution density of the characteristic point group to obtain the propagation characteristic intensity; Drawing the pressure wave propagation path, labeling the propagation characteristic intensity, marking the reflection point position and reflection direction, marking the attenuation point position and attenuation coefficient, and generating a pressure wave propagation characteristic distribution map.

5. The method of claim 1, wherein, Monitoring the fluid flow state parameters of each cross section of the catheter through the microfluid sensor array built in the catheter, identifying the local turbulent flow region according to the fluid flow state parameters and the multi-dimensional parameter feature matrix, and calculating the turbulent flow development trend when the turbulent intensity of the local turbulent flow region exceeds the early warning intensity threshold, including: Collecting the pressure fluctuation value and the flow velocity distribution value of each cross section and performing mapping operation with the multi-dimensional parameter feature matrix to obtain the strain rate tensor and the vorticity tensor, calculating the second invariant of the vorticity tensor and the second invariant of the rate of change tensor respectively and making a difference to obtain the turbulent flow discrimination factor, and marking the local turbulent flow region based on the turbulent flow discrimination factor; Extracting the flow velocity fluctuation component in the local turbulent flow region, calculating the ratio of the root mean square value of the flow velocity fluctuation component to the average flow velocity to obtain the turbulent intensity; When the turbulence intensity exceeds the early warning intensity threshold, the flow velocity fluctuation signal of the local turbulence region is extracted and wavelet transformed, the wavelet coefficient matrix is integrated at different scales, the square of each element in the wavelet coefficient matrix is calculated to obtain the local wave energy, the scale energy is obtained by summing the local wave energy at the same scale, and the energy distribution interval is extracted based on the scale energy; In the energy distribution interval, the difference between the turbulence energies of adjacent time instants is divided by the time interval to obtain the energy change rate, the difference between the characteristic lengths of adjacent scales is divided by the time interval to obtain the scale change rate, the turbulence development index is obtained based on the energy change rate and the scale change rate, and the turbulence development trend is determined according to the time sequence change characteristics of the turbulence development index.

6. The method of claim 1, wherein, According to the turbulence development trend and the frequency domain characteristics of the pressure fluctuation, the pressure fluctuation is eliminated, and the pressure in the catheter is controlled to maintain within a safe pressure range, including: Extracting the pressure fluctuation signal in the catheter, decomposing the pressure fluctuation signal into periodic fluctuation components in the time domain, and obtaining the fluctuation main frequency and its amplitude according to the Fourier series decomposition result of the periodic fluctuation components; Setting a catheter pressure target value, and taking the difference between the real-time detected catheter pressure value and the catheter pressure target value as a pressure deviation; According to the turbulence development trend, the fluctuation main frequency and its amplitude, a proportional coefficient, an integral time constant and a differential time constant are calculated, the proportional coefficient is multiplied by the pressure deviation to obtain a proportional control term, the integral time constant is multiplied by the pressure deviation integral value to obtain an integral control term, the differential time constant is multiplied by the pressure deviation change rate to obtain a differential control term, and the proportional control term, the integral control term and the differential control term are added to obtain a compensation control signal; The compensation control signal is input into a piezoelectric ceramic diaphragm pump to generate a compensation pressure, and the catheter pressure value after compensation is detected in real time. When the catheter pressure value after compensation exceeds the upper limit value or the lower limit value of the safe pressure range, the change rate of the pressure deviation is calculated, and the compensation pressure is adjusted according to the change rate of the pressure deviation, so that the pressure value in the catheter is maintained within the safe pressure range.

7. A real-time monitoring and control system for intraductal pressure in a ductoscopy procedure for implementing the method according to any one of claims 1 to 6, characterized in that, Including: A first unit for collecting real-time temperature data, conductivity data and image data in the catheter during a breast ductoscopy process and forming a multi-dimensional parameter feature matrix; A second unit for performing adaptive filtering analysis on the multi-dimensional parameter feature matrix, extracting the characteristic parameters of sub-millimeter fluid particles in the fluid, calculating the spatial distribution density value and motion velocity vector of the sub-millimeter fluid particles according to the characteristic parameters, calculating the spatial distribution characteristics of the pressure in the catheter according to the spatial distribution density value, determining the positions of the pressure peak point and the pressure valley point, and identifying the pressure gradient mutation region; tracking the pressure wave propagation path according to the motion velocity vector to identify the reflection point and the attenuation point of the pressure wave on the propagation path, and obtaining a pressure wave propagation characteristic distribution map; A third unit for real-time monitoring of the pressure in the catheter in combination with the pressure gradient mutation region and the pressure wave propagation characteristic distribution map. A fourth unit is configured to monitor fluid flow state parameters of each section of the conduit by a microfluid sensor array built in the conduit, identify a local turbulent flow region according to the fluid flow state parameters and the multi-dimensional parameter feature matrix, and calculate a turbulent flow development trend when detecting that a turbulent flow intensity of the local turbulent flow region exceeds a warning intensity threshold. According to the turbulent flow development trend and a frequency domain feature of pressure fluctuation, the pressure fluctuation is eliminated, and the pressure in the conduit is controlled to be maintained within a safe pressure range.

8. An electronic device, comprising: The system comprises: a processor; a memory for storing processor-executable instructions; wherein the processor is configured to invoke the instructions stored in the memory to execute the method of any one of claims 1 to 6.

9. A computer-readable storage medium having stored thereon computer program instructions, wherein, The computer program instructions, when executed by the processor, implement the method of any one of claims 1 to 6.

Images