Method and system for monitoring and controlling pressure in catheter in real time in breast duct endoscopy

Abstract

The invention provides a real-time monitoring and control method and system for pressure in a catheter in breast duct endoscopy, and relates to the technical field of medical equipment monitoring and control, and the method comprises the following steps: forming a multi-dimensional parameter characteristic matrix by collecting temperature, conductivity and image data in the catheter, and carrying out adaptive filtering analysis to extract fluid characteristic parameters; pressure distribution characteristics in the conduit are calculated, fluid flow state parameters are monitored in combination with a microfluid sensor array to identify a local turbulence area, when the turbulence intensity exceeds an early warning threshold value, wavelet transform is adopted to analyze and calculate the turbulence development trend, and the pressure in the conduit is controlled within a safe range through an active pulsation compensation technology; and the safety and the accuracy of breast duct endoscope examination are improved.

Description

Technical Field

[0001] This invention relates to the field of medical device monitoring and control technology, and in particular to a method and system for real-time monitoring and control of intracatheter pressure during ductoscopy. Background Technology

[0002] Ductoscopy is a minimally invasive technique used to diagnose breast diseases. A thin endoscope is inserted through the nipple into the mammary ducts to observe the internal structure and lesions. This technique is widely used in the diagnosis of intraductal papilloma, duct ectasia, and early-stage breast cancer. During ductoscopy, a saline solution or other irrigation fluid is injected into the mammary ducts to dilate them and facilitate observation of the internal structure.

[0003] During ductoscopy, controlling the intraductal pressure is crucial. Appropriate intraductal pressure ensures a good field of vision and accurate examination, while excessive pressure can lead to breast tissue damage, increased pain during the procedure, and even fluid leakage into surrounding tissues, causing complications. Traditional ductoscopy relies heavily on the operator's experience to control the injection rate and pressure, lacking precise pressure monitoring and control methods.

[0004] Existing pressure control methods for ductoscopy have the following shortcomings: They lack real-time, precise intra-catheter pressure monitoring technology, making it difficult for operators to accurately grasp the actual pressure state within the catheter, easily leading to excessively high or low pressure, affecting examination results and patient safety; existing methods cannot effectively identify and predict local turbulence and pressure fluctuations within the catheter, and these abnormal hydrodynamic phenomena can cause tissue damage or decreased image quality; traditional pressure control mostly employs passive regulation, lacking an intelligent active pressure control system, and cannot adaptively adjust pressure according to the real-time fluid state within the catheter, especially when facing complex catheter morphologies or lesions, making it difficult to maintain a stable pressure environment. Summary of the Invention

[0005] This invention provides a method and system for real-time monitoring and control of intraductal pressure during ductoscopy, which can solve the problems in the prior art.

[0006] A first aspect of the present invention provides a method for real-time monitoring and control of intraductal pressure during ductoscopy, comprising: Real-time temperature data, conductivity data, and image data inside the catheter were collected during ductoscopy and assembled into a multidimensional parameter feature matrix. Adaptive filtering analysis is performed on the multidimensional parameter feature matrix to extract the feature parameters of submillimeter fluid particles in the fluid. Based on the feature parameters, the spatial distribution density value and velocity vector of the submillimeter fluid particles are calculated. Based on the spatial distribution density value, the spatial distribution characteristics of the pressure inside the duct are calculated to determine the location of pressure peak points and pressure valley points, and to identify regions of abrupt pressure gradient changes. Based on the velocity vector, the propagation path of the pressure wave is traced to identify the reflection points and attenuation points of the pressure wave on the propagation path, thus obtaining a pressure wave propagation characteristic distribution map. By combining the distribution map of pressure gradient abrupt change regions and pressure wave propagation characteristics, the pressure inside the catheter can be monitored in real time; The fluid flow parameters of each cross section of the catheter are monitored by a microfluidic sensor array built into the catheter. Local turbulent regions are identified based on the fluid flow parameters and the multidimensional parameter feature matrix. When the turbulence intensity of the local turbulent region exceeds the warning intensity threshold, the turbulence development trend is calculated. Based on the turbulence development trend and the frequency domain characteristics of pressure fluctuations, pressure fluctuations are eliminated, and the pressure inside the catheter is controlled to be maintained within a safe pressure range.

[0007] Adaptive filtering analysis is performed on the multidimensional parameter feature matrix to extract the feature parameters of submillimeter fluid particles in the fluid. The spatial distribution density and velocity vector of the submillimeter fluid particles are calculated based on these feature parameters, including: The autocorrelation matrix of the multidimensional parameter feature matrix is ​​calculated by adaptive Wiener filtering to obtain the filtered feature data. In the filtered feature data, edge points are detected using a multi-scale operator, and the edge points are morphologically connected to obtain a complete contour. The equivalent diameter and roundness coefficient are calculated based on the complete contour, and the gray-level gradient value of the contour region is calculated at the same time. The equivalent diameter, roundness coefficient, and gray-level gradient value are combined to form the feature parameters of the submillimeter fluid particle. The initial bandwidth is determined based on the characteristic parameters. The initial bandwidth is then cross-validated and iteratively optimized until the density estimation error is less than the distribution fitting threshold, resulting in the optimized bandwidth parameters. The spatial distribution density of submillimeter fluid particles is then calculated using the optimized bandwidth parameters. Extract images of adjacent frames from the filtered feature data, calculate the spatial and temporal gradients of the adjacent frame images using the feature parameters, and construct a velocity field smoothing constraint. Use the spatial gradient, temporal gradient, and velocity field smoothing constraint to construct an optimization objective function. Iterate through the gradient descent method until the value of the optimization objective function is less than the velocity field convergence threshold, and determine the motion velocity vector of the submillimeter fluid particle.

[0008] Based on the spatial distribution density value, calculate the spatial distribution characteristics of the pressure inside the catheter, determine the locations of pressure peak points and pressure trough points, and identify regions of abrupt pressure gradient changes, including: Calculate the initial pressure distribution based on the spatial density value within the catheter; The Laplacian operator for density is calculated, and a density gradient correction term is obtained based on the Laplacian operator and the density gradient coefficient. The density gradient correction term is then superimposed on the initial pressure distribution to obtain the corrected pressure distribution. Multi-scale analysis is performed on the corrected pressure distribution, 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. Local extreme points where the gradient vector field shows convergence are marked as pressure peak points, and local extreme points where the gradient vector field shows divergence are marked as pressure valley points. Based on the distribution locations of the pressure peak points and pressure valley points, the gradient amplitude of the corrected pressure distribution at each location is calculated, and the gradient direction angle between adjacent locations is calculated to obtain the direction change rate. The gradient amplitude is divided by the average gradient amplitude in the local area to obtain the amplitude ratio. The gradient mutation index is calculated based on the amplitude ratio and the direction change rate. When the gradient mutation index is greater than the mutation judgment threshold, the area between adjacent pressure peak points and pressure valley points is marked as a pressure gradient mutation region.

[0009] Based on the motion velocity vector, the propagation path of the pressure wave is traced, and the reflection points and attenuation points of the pressure wave on the propagation path are identified to obtain the pressure wave propagation characteristic distribution map, including: Extract the velocity components of discrete points in the corrected pressure distribution, calculate the pressure wave acceleration based on the velocity components of adjacent discrete points, and integrate the pressure wave acceleration over time to obtain the pressure wave propagation path. The inner product of the velocity direction at each moment and the velocity direction at the next moment is calculated along the propagation path of the pressure wave. The propagation direction change angle is obtained based on the inner product and the velocity amplitude. When the propagation direction change angle is greater than the direction change threshold, the corresponding coordinate point is marked as the reflection point. 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. The pressure gradient along the path is divided by the pressure value to obtain the attenuation coefficient. When the attenuation coefficient is greater than the attenuation characteristic threshold, the corresponding coordinate point is marked as an attenuation point. Calculate the distance from each reflection point to each attenuation point, pair reflection points and attenuation points whose distance is less than the feature association threshold to form a feature point group, and calculate the spatial distribution density of the feature point group to obtain the propagation feature intensity. Draw the propagation path of the pressure wave, mark the propagation characteristic intensity, mark the location and direction of the reflection point, mark the location and attenuation coefficient of the attenuation point, and generate a distribution map of the propagation characteristics of the pressure wave.

[0010] Fluid flow parameters at various cross-sections of the catheter are monitored using a microfluidic sensor array embedded within the catheter. Local turbulent regions are identified based on these fluid flow parameters and the multidimensional parameter feature matrix. When the turbulence intensity in a local turbulent region exceeds a warning intensity threshold, the turbulence development trend is calculated, including: The pressure fluctuation values ​​and flow velocity distribution values ​​of each section are collected and mapped with the multidimensional parameter feature matrix to obtain the strain rate tensor and vorticity tensor. The second-order invariants of the vorticity tensor and the second-order invariants of the strain rate tensor are calculated and subtracted to obtain the turbulence discrimination factor. The local turbulence region is marked based on the turbulence discrimination factor. Extract the velocity fluctuation component in the local turbulent region, and calculate the ratio of the root mean square value of the velocity fluctuation component to the average velocity to obtain the turbulence intensity. When the turbulence intensity exceeds the warning intensity threshold, the velocity fluctuation signal of the local turbulence region is extracted and wavelet transform is performed. The wavelet coefficient matrix is ​​obtained by integrating at different scales. The square of each element in the wavelet coefficient matrix is ​​calculated to obtain the local wave energy. The local wave energy is summed at the same scale to obtain the scale energy. The energy distribution interval is extracted based on the scale energy. Within the energy distribution range, the energy change rate is obtained by dividing the difference in turbulent energy at adjacent times by the time interval, and the scale change rate is obtained by dividing the difference in characteristic length of adjacent scales by the time interval. Based on the energy change rate and the scale change rate, a turbulence development index is obtained, and the turbulence development trend is judged based on the time series change characteristics of the turbulence development index.

[0011] Based on the aforementioned turbulence development trend and the frequency domain characteristics of pressure fluctuations, eliminating pressure fluctuations and controlling the pressure inside the duct to maintain within a safe pressure range includes: The pressure fluctuation signal inside the catheter is extracted, and the pressure fluctuation signal is decomposed into periodic fluctuation components in the time domain. The main frequency and amplitude of the fluctuation are obtained based on the Fourier series decomposition result of the periodic fluctuation components. Set a target value for catheter pressure, and use the difference between the real-time detected catheter pressure value and the target value for catheter pressure as the pressure deviation; Based on the turbulence development trend, the dominant frequency of the wave and its amplitude, calculate the proportional coefficient, integral time constant and derivative time constant. Multiply the proportional coefficient by the pressure deviation to obtain the proportional control term, multiply the integral time constant by the integral value of the pressure deviation to obtain the integral control term, multiply the derivative time constant by the rate of change of the pressure deviation to obtain the derivative control term, and add the proportional control term, integral control term and derivative control term to obtain the compensation control signal. The compensation control signal is input to the piezoelectric ceramic diaphragm pump to generate compensation pressure. The pressure value of the compensated conduit is detected in real time. When the pressure value of the compensated conduit exceeds the upper or lower safety limit, the rate of change of the pressure deviation is calculated. The compensation pressure is adjusted according to the rate of change of the pressure deviation to keep the pressure value in the conduit within the safe pressure range.

[0012] A second aspect of the present invention provides a real-time monitoring and control system for intraductal pressure during ductoscopy, comprising: The first unit is used to collect real-time temperature data, conductivity data and image data inside the catheter during ductoscopy and form a multi-dimensional parameter feature matrix. The second unit is used to perform adaptive filtering analysis on the multidimensional parameter feature matrix, extract the feature parameters of submillimeter fluid particles in the fluid, calculate the spatial distribution density value and velocity vector of the submillimeter fluid particles based on the feature parameters, calculate the spatial distribution characteristics of the pressure inside the duct based on the spatial distribution density value, determine the location of pressure peak points and pressure valley points, and identify regions of abrupt pressure gradient changes; track the pressure wave propagation path based on the velocity vector, identify the reflection points and attenuation points of the pressure wave on the propagation path, and obtain a pressure wave propagation characteristic distribution map. The third unit is used to monitor the pressure inside the catheter in real time by combining the pressure gradient abrupt region and the pressure wave propagation characteristic distribution map. The fourth unit is used to monitor the fluid flow parameters of each cross section of the catheter through a microfluidic sensor array built into the catheter, identify local turbulent regions based on the fluid flow parameters and the multidimensional parameter feature matrix, calculate the turbulence development trend when the turbulence intensity of the local turbulence region exceeds the warning intensity threshold, and eliminate pressure fluctuations based on the turbulence development trend and the frequency domain characteristics of pressure fluctuations to control the pressure inside the catheter to be maintained within a safe pressure range.

[0013] A third aspect of the embodiments of the present invention, An electronic device is provided, comprising: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0014] Fourth aspect of the present invention, A computer-readable storage medium is provided, having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0015] The beneficial effects of this application are as follows: This invention achieves accurate real-time monitoring of intraductal pressure during ductoscopy through multidimensional parameter feature matrix analysis and adaptive filtering technology, enabling timely detection of abnormal pressure areas and preventing tissue damage risks caused by excessive pressure.

[0016] By using a microfluidic sensor array and fluid flow parameter analysis, this invention can accurately identify local turbulent regions and predict the development trend of turbulence through wavelet transform analysis, thereby achieving proactive early warning of pressure within the duct and improving the safety of ductoscopy.

[0017] By employing aortic pressure compensation technology to control intracatheter pressure, this invention can effectively eliminate intracatheter pressure fluctuations, ensuring that the pressure is always maintained within a safe range. This significantly reduces the risk of tissue damage during ductoscopy, and improves the comfort and diagnostic accuracy of the examination. Attached Figure Description

[0018] Figure 1 This is a flowchart illustrating the method for real-time monitoring and control of intraductal pressure during ductoscopy according to an embodiment of the present invention. Figure 2 This is a diagram showing the performance comparison of adaptive filtering. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

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

[0021] Figure 1 This is a flowchart illustrating the method for real-time monitoring and control of intraductal pressure during ductoscopy according to an embodiment of the present invention. Figure 1 As shown, the method includes: Real-time temperature data, conductivity data, and image data inside the catheter were collected during ductoscopy and assembled into a multidimensional parameter feature matrix. Adaptive filtering analysis is performed on the multidimensional parameter feature matrix to extract the feature parameters of submillimeter fluid particles in the fluid. Based on the feature parameters, the spatial distribution density value and velocity vector of the submillimeter fluid particles are calculated. Based on the spatial distribution density value, the spatial distribution characteristics of the pressure inside the duct are calculated to determine the location of pressure peak points and pressure valley points, and to identify regions of abrupt pressure gradient changes. Based on the velocity vector, the propagation path of the pressure wave is traced to identify the reflection points and attenuation points of the pressure wave on the propagation path, thus obtaining a pressure wave propagation characteristic distribution map. By combining the distribution map of pressure gradient abrupt change regions and pressure wave propagation characteristics, the pressure inside the catheter can be monitored in real time; The fluid flow parameters of each cross section of the catheter are monitored by a microfluidic sensor array built into the catheter. Local turbulent regions are identified based on the fluid flow parameters and the multidimensional parameter feature matrix. When the turbulence intensity of the local turbulent region exceeds the warning intensity threshold, the turbulence development trend is calculated. Based on the turbulence development trend and the frequency domain characteristics of pressure fluctuations, pressure fluctuations are eliminated, and the pressure inside the catheter is controlled to be maintained within a safe pressure range.

[0022] In one optional implementation, adaptive filtering analysis is performed on the multidimensional parameter feature matrix to extract feature parameters of submillimeter fluid particles in the fluid. The spatial distribution density and velocity vector of the submillimeter fluid particles are then calculated based on these feature parameters, including: The autocorrelation matrix of the multidimensional parameter feature matrix is ​​calculated by adaptive Wiener filtering to obtain the filtered feature data. In the filtered feature data, edge points are detected using a multi-scale operator, and the edge points are morphologically connected to obtain a complete contour. The equivalent diameter and roundness coefficient are calculated based on the complete contour, and the gray-level gradient value of the contour region is calculated at the same time. The equivalent diameter, roundness coefficient, and gray-level gradient value are combined to form the feature parameters of the submillimeter fluid particle. The initial bandwidth is determined based on the characteristic parameters. The initial bandwidth is then cross-validated and iteratively optimized until the density estimation error is less than the distribution fitting threshold, resulting in the optimized bandwidth parameters. The spatial distribution density of submillimeter fluid particles is then calculated using the optimized bandwidth parameters. Extract images of adjacent frames from the filtered feature data, calculate the spatial and temporal gradients of the adjacent frame images using the feature parameters, and construct a velocity field smoothing constraint. Use the spatial gradient, temporal gradient, and velocity field smoothing constraint to construct an optimization objective function. Iterate through the gradient descent method until the value of the optimization objective function is less than the velocity field convergence threshold, and determine the motion velocity vector of the submillimeter fluid particle.

[0023] Calculate the autocorrelation matrix R of the feature matrix, which is N×N, where N is the dimension of the feature vector. In practical applications, taking a 5-dimensional feature vector as an example, the autocorrelation matrix is ​​a 5×5 matrix. Calculate the trace Tr(R) of the autocorrelation matrix, which reflects the overall energy of the data. Dynamically adjust the filter step size μ based on the trace value. When the trace value is large (e.g., greater than 100), the step size is set to 0.01; when the trace value is small (e.g., less than 50), the step size is set to 0.05 to ensure algorithm stability and convergence speed. Iteratively update the filter coefficients w, using an iteration formula based on the minimum mean square error criterion. Update the coefficients in each iteration until the change in coefficients between two adjacent iterations is less than a preset threshold (e.g., 0.001) or the maximum number of iterations is reached (e.g., 50 times). Finally, the filtered feature data is obtained, which significantly reduces noise and preserves the edge features of submillimeter bubbles and colloidal clusters.

[0024] Multi-scale processing is achieved by setting different Gaussian kernel standard deviations σ, typically three scale values: σ = 0.8, 1.2, and 1.6. Image gradients are calculated for each scale, and non-maximum suppression is applied. Edge points are determined using a dual-threshold method, with the lower threshold set to 40% of the higher threshold. The higher threshold is automatically determined using the OTSU method, usually at the 85th percentile of the gradient histogram. Edge points detected at different scales are fused to obtain a complete set of edge points. Morphological connectivity processing is applied to the edge points, including thinning, closure, and connection operations, with a connection distance threshold set to 5 pixels, forming a complete target contour. Geometric feature parameters are calculated based on the contour: area S (number of pixels), perimeter L (contour length), and equivalent diameter. And the roundness coefficient C=4πS / L 2 Simultaneously, the average gray-level gradient G and standard deviation σG within the contour region are calculated. These parameters constitute the characteristic parameter set for submillimeter bubbles and colloidal clusters. In practical cases, the typical equivalent diameter of submillimeter bubbles ranges from 0.1mm to 0.9mm, with a roundness coefficient greater than 0.85, while colloidal clusters have similar equivalent diameters but a roundness coefficient typically lower than 0.75, and the gray-level gradient value of bubbles is usually higher than that of clusters.

[0025] The initial Gaussian kernel bandwidth h0 is chosen based on the Silverman criterion: h0 = 0.9 × min(σ, IQR / 1.34) × n -1 / 5Where σ is the sample standard deviation, IQR is the interquartile range, and n is the number of samples. The initial bandwidth is optimized using k-fold cross-validation (k=5). The data is divided into 5 parts, and the density function is estimated using 4 parts each time, with the remaining part used to calculate the estimation error. Within the bandwidth range [0.5h0, 1.5h0], the bandwidth value that minimizes the error is selected with a step size of 0.1h0. Iterative optimization continues until the density estimation error is less than a preset threshold (e.g., 0.01). The optimized bandwidth h* is used to calculate the density value ρ(x,y) at each point (x,y) in two-dimensional space. In practical applications, the bubble distribution density in fluids typically ranges from 0 to 25 bubbles / cm². 2 High-density regions correspond to regions of fluid dynamic instability.

[0026] Extract adjacent frame images I(x,y,t) and I(x,y,t+Δt) from the filtered feature data, typically Δt = 1 / 30 second (corresponding to a 30fps acquisition rate). Combine the position and shape information from the feature parameters to calculate the spatial gradients Ix, Iy, and temporal gradient It of adjacent frame images. Construct an optical flow constraint equation and add a velocity field smoothing constraint to form the optimization objective function J. Iterate the velocity field using the gradient descent method, updating the step size to 0.05, iterating until the change in the objective function is less than a preset threshold (e.g., 0.0005) or the maximum number of iterations (e.g., 100). Finally, obtain the velocity vector (u,v) for each bubble and cluster, representing its velocity in the x and y directions. In practical case analysis, the typical velocity range of bubbles in fluids is 0-15 cm / s, and the velocity field distribution reflects the internal flow state and turbulence characteristics of 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.

[0027] Figure 2This diagram illustrates the performance comparison of adaptive filtering. It shows the performance of the proposed adaptive filtering algorithm compared to four classic filtering methods over 50 iterations. The vertical axis represents the filtering error (in dB), with smaller negative values ​​indicating better noise suppression. The diagram clearly shows that the proposed adaptive Wiener filtering algorithm (solid black line) exhibits the best filtering performance, rapidly converging from an initial -9.2 dB to -30.3 dB. This represents a 12.9 dB improvement compared to the traditional Wiener filter's -17.4 dB, a 6.8 dB improvement compared to the RLS adaptive filter's -23.5 dB, a 14.3 dB improvement compared to the LMS adaptive filter's -16.0 dB, and an 8.7 dB improvement compared to the Kalman filter's -21.6 dB. Particularly noteworthy is the faster convergence speed observed in the first 20 iterations. This is attributed to the dynamic step-size adjustment mechanism based on the autocorrelation matrix trace, which adaptively optimizes the filter parameters according to data characteristics. In the later stages of iteration, the present invention can continue to improve the filtering effect and maintain stable convergence, while other methods reach the performance bottleneck at an earlier stage. This fully verifies the technical advantages and practical value of the present invention in the feature extraction of submillimeter bubbles and colloidal clusters.

[0028] In one optional implementation, the spatial distribution characteristics of the pressure within the catheter are calculated based on the spatial distribution density value, the locations of pressure peak points and pressure trough points are determined, and regions of abrupt pressure gradient changes are identified, including: Calculate the initial pressure distribution based on the spatial density value within the catheter; The Laplacian operator for density is calculated, and a density gradient correction term is obtained based on the Laplacian operator and the density gradient coefficient. The density gradient correction term is then superimposed on the initial pressure distribution to obtain the corrected pressure distribution. Multi-scale analysis is performed on the corrected pressure distribution, 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. Local extreme points where the gradient vector field shows convergence are marked as pressure peak points, and local extreme points where the gradient vector field shows divergence are marked as pressure valley points. Based on the distribution locations of the pressure peak points and pressure valley points, the gradient amplitude of the corrected pressure distribution at each location is calculated, and the gradient direction angle between adjacent locations is calculated to obtain the direction change rate. The gradient amplitude is divided by the average gradient amplitude in the local area to obtain the amplitude ratio. The gradient mutation index is calculated based on the amplitude ratio and the direction change rate. When the gradient mutation index is greater than the mutation judgment threshold, the area between adjacent pressure peak points and pressure valley points is marked as a pressure gradient mutation region.

[0029] In methods for calculating the pressure distribution characteristics within a catheter, it is necessary to obtain the spatial density values ​​within the catheter. These values ​​can be acquired by sensors at different locations within the catheter or obtained through computational fluid dynamics simulations. For example, in a catheter with a length of 200 mm, density values ​​can be measured or calculated at 10 mm intervals, forming a density distribution data set containing 21 measurement points.

[0030] When calculating the initial pressure distribution based on the obtained density values ​​within the conduit, a transformation can be performed using equations of state. For liquids, a relationship similar to p = ρgh can be applied, where ρ is density, g is gravitational acceleration, and h is the height of the liquid column; for gases, a relationship similar to p = ρRT / M can be applied, where R is the gas constant, T is temperature, and M is molecular weight. Taking a gas as an example, if the density at a certain measurement point is 1.2 kg / m³, the temperature is 293 Kelvin, and the molecular weight is 29, then the initial pressure at that point can be calculated to be approximately 101325 Pa.

[0031] In discrete space, this can be achieved using the finite difference method. For the one-dimensional case, the three-point difference formula can be used; for the three-dimensional case, the second derivatives can be calculated separately in the x, y, and z directions, and then summed. For example, in the one-dimensional case, the Laplace operator at a point i can be expressed as (ρi+1 - 2ρi + ρi-1) / Δx 2 Calculate the density, where Δx is the spatial sampling interval. If the density at point i is 1.2 kg / m³, at point i+1 it is 1.25 kg / m³, at point i-1 it is 1.18 kg / m³, and Δx is 10 mm, then the Laplace operator value is (1.25 - 2 × 1.2 + 1.18) / (0.01). 2 =300 kg / m³ / m 2 .

[0032] When multiplying the Laplace operator by the density gradient coefficient to obtain the density gradient correction term, an appropriate density gradient coefficient needs to be selected. This coefficient can be determined based on the fluid properties and the conduit material, with typical values ​​between 0.001 and 0.1. If a density gradient coefficient of 0.01 is chosen, the density gradient correction term obtained from the Laplace operator calculated above will be 300 × 0.01 = 3 Pa. This correction term is then superimposed on the initial pressure distribution; for example, if the initial pressure at this point is 101325 Pa, the corrected pressure will be 101325 + 3 = 101328 Pa.

[0033] Three scales with standard deviations of 1 mm, 2 mm, and 4 mm are set. At each scale, local extrema are extracted by comparing the pressure values ​​of each point with those of its neighbors. If the pressure value of a point is greater than the pressure values ​​of all its neighbors, then that point is a local maximum; if it is less than the pressure values ​​of all its neighbors, then that point is a local minimum.

[0034] For a point (x, y, z) in three-dimensional space, its pressure gradient vector contains three components: x, y, and z, representing the rate of change of pressure in these three directions, respectively. The convergence and divergence of the gradient vector can be determined by calculating its divergence. A negative divergence indicates that the gradient vector field converges at that point, which is a pressure peak; a positive divergence indicates that the gradient vector field diverges at that point, which is a pressure trough. In practical applications, a threshold value can be set, such as ±0.5 Pa / mm². 2 A peak or valley point is considered valid only when the absolute value of the divergence is greater than the threshold.

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

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

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

[0038] When calculating the gradient abrupt change index using the weighted sum of the magnitude ratio and the rate of change of direction, it is necessary to determine the weighting coefficients for both. These weights can be adjusted according to the specific application scenario. For example, the weight of the magnitude ratio can be set to 0.7, and the weight of the rate of change of direction to 0.3. If the magnitude ratio is 1.42 and the rate of change of direction is 0.785 radians, then the gradient abrupt change index is approximately 0.7 × 1.42 + 0.3 × 0.785 ≈ 1.229.

[0039] When the gradient abrupt change index exceeds the abrupt change threshold, the region between adjacent pressure peaks and valleys is marked as a pressure gradient abrupt change region. The abrupt change threshold can be determined based on historical data and expert experience, for example, set to 1.0. In the example above, the gradient abrupt change index is 1.229, which is greater than the threshold of 1.0; therefore, this region is marked as a pressure gradient abrupt change region. These abrupt change regions typically correspond to locations where significant changes occur in the fluid dynamics characteristics within the catheter, such as changes in cross-sectional area, bends, or bifurcations, and are of great significance for catheter design and fluid control.

[0040] In one optional implementation, the pressure wave propagation path is traced based on the motion velocity vector, and the reflection points and attenuation points of the pressure wave on the propagation path are identified to obtain a pressure wave propagation characteristic distribution map, including: Extract the velocity components of discrete points in the corrected pressure distribution, calculate the pressure wave acceleration based on the velocity components of adjacent discrete points, and integrate the pressure wave acceleration over time to obtain the pressure wave propagation path. The inner product of the velocity direction at each moment and the velocity direction at the next moment is calculated along the propagation path of the pressure wave. The propagation direction change angle is obtained based on the inner product and the velocity amplitude. When the propagation direction change angle is greater than the direction change threshold, the corresponding coordinate point is marked as the reflection point. 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. The pressure gradient along the path is divided by the pressure value to obtain the attenuation coefficient. When the attenuation coefficient is greater than the attenuation characteristic threshold, the corresponding coordinate point is marked as an attenuation point. Calculate the distance from each reflection point to each attenuation point, pair reflection points and attenuation points whose distance is less than the feature association threshold to form a feature point group, and calculate the spatial distribution density of the feature point group to obtain the propagation feature intensity. Draw the propagation path of the pressure wave, mark the propagation characteristic intensity, mark the location and direction of the reflection point, mark the location and attenuation coefficient of the attenuation point, and generate a distribution map of the propagation characteristics of the pressure wave.

[0041] For the acquired pressure distribution data, the pressure data of each discrete point is first extracted and corrected. The corrected pressure distribution data includes the position coordinates of each discrete point in the time series and the corresponding pressure value. From this data, the velocity components in the x and y directions of each discrete point are extracted. For example, at a certain time t1, the velocity component in the x direction of discrete point P1 is 3.5 m / s, and the velocity component in the y direction is 2.1 m / s. The pressure wave acceleration is calculated based on the velocity components at adjacent times. For example, for adjacent times t1 and t2, the acceleration ax is calculated as (vx2-vx1) / (t2-t1), and ay is calculated similarly. Specifically, if the x-direction velocity of point P1 at time t1 is 3.5 m / s, and at time t2 it is 4.2 m / s, with a time interval of 0.01 seconds, then the x-direction acceleration is 70 m / s². 2 The propagation path of the pressure wave is obtained by integrating the acceleration of the pressure wave over time, that is, by accumulating the displacement increments at each moment to construct the complete propagation path.

[0042] The change in propagation direction along the pressure wave propagation path is calculated. For each moment along the propagation path, the dot product of the current velocity direction and the next velocity direction is calculated. For example, if the velocity vector at time t1 is (3.5, 2.1) and the velocity vector at time t2 is (3.8, 1.9), then the dot product is 3.5 × 3.8 + 2.1 × 1.9 = 17.29. This dot product is divided by the product of the magnitudes of the two velocity vectors to obtain the direction cosine. In this example, the magnitude of the velocity vector at time t1 is 4.07 and at time t2 it is 4.25, so the direction cosine is 17.29 / (4.07 × 4.25) = 0.9996. Converting the direction cosine to radians, the change angle of the propagation direction is arccos(0.9996) = 0.0283 radians, approximately 1.62 degrees. When the angle of change of direction is greater than a preset threshold for abrupt change of direction (e.g., set to 10 degrees or 0.1745 radians), the corresponding coordinate point is marked as a reflection point. In this example, the point is not marked as a reflection point.

[0043] For pressure attenuation analysis along the propagation path, the pressure difference between adjacent coordinate points is calculated and divided by the distance between the two points to obtain the pressure gradient along the path. For example, if the pressure at point P1 at time t1 is 120 kPa and the pressure at point P2 at time t2 is 115 kPa, and the distance between the two points is 0.05 m, then the pressure gradient along the path is (120-115) / 0.05 = 100 kPa / m. Dividing the pressure gradient along the path by the pressure value yields the attenuation coefficient, which in this example is 100 / 120 = 0.833 m. -1 When the attenuation coefficient is greater than the preset attenuation characteristic threshold (e.g., set to 0.5 m), -1 When the attenuation point is reached, the corresponding coordinate point is marked as the attenuation point. In this example, the point is marked as the attenuation point.

[0044] Feature association analysis is performed on the labeled reflection and attenuation points. The distance from each reflection point to each attenuation point is calculated. When the distance is less than the feature association threshold (e.g., set to 0.2 m), the reflection point and attenuation point are 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... The distance between these two points is less than the feature association threshold of 0.2 m, therefore they form a feature point group. The distribution density of all feature point groups in space is statistically analyzed to obtain the propagation feature intensity. For example, if there are 5 feature point groups in an area of ​​1 square meter, then the propagation feature intensity of that area is 5 points / square meter.

[0045] Finally, the propagation path of the pressure wave is plotted in a two-dimensional coordinate system. The intensity of the propagation characteristics is represented by the shade of color or the thickness of the line; stronger characteristics are represented by darker colors or thicker lines. The locations of reflection points are marked on the graph, and arrows indicate the direction of reflection. At the same time, the locations of attenuation points are marked, and the attenuation coefficient is numerically labeled. For example, reflection point R1 is marked with a red dot, and arrows indicate the directions before reflection (0.8, 0.6) and after reflection (0.6, -0.8); attenuation point A1 is marked with a blue triangle, and the attenuation coefficient 0.833 m is labeled. -1 The pressure wave propagation characteristic distribution map generated in this way visually demonstrates the reflection and attenuation characteristics of the pressure wave during propagation, which helps to analyze the propagation law of pressure waves and predict their propagation behavior.

[0046] The above method achieves comprehensive analysis and visualization of pressure wave propagation characteristics through precise tracking of the pressure wave propagation path and identification of feature points, providing an effective technical means for pressure wave propagation research.

[0047] In one optional implementation, a microfluidic sensor array embedded in the catheter monitors the fluid flow parameters at various cross-sections of the catheter. Local turbulent regions are identified based on these fluid flow parameters and the multidimensional parameter feature matrix. When the turbulence intensity in the detected local turbulent region exceeds a warning intensity threshold, the turbulence development trend is calculated, including: The pressure fluctuation values ​​and flow velocity distribution values ​​of each section are collected and mapped with the multidimensional parameter feature matrix to obtain the strain rate tensor and vorticity tensor. The second-order invariants of the vorticity tensor and the second-order invariants of the strain rate tensor are calculated and subtracted to obtain the turbulence discrimination factor. The local turbulence region is marked based on the turbulence discrimination factor. Extract the velocity fluctuation component in the local turbulent region, and calculate the ratio of the root mean square value of the velocity fluctuation component to the average velocity to obtain the turbulence intensity. When the turbulence intensity exceeds the warning intensity threshold, the velocity fluctuation signal of the local turbulence region is extracted and wavelet transform is performed. The wavelet coefficient matrix is ​​obtained by integrating at different scales. The square of each element in the wavelet coefficient matrix is ​​calculated to obtain the local wave energy. The local wave energy is summed at the same scale to obtain the scale energy. The energy distribution interval is extracted based on the scale energy. Within the energy distribution range, the energy change rate is obtained by dividing the difference in turbulent energy at adjacent times by the time interval, and the scale change rate is obtained by dividing the difference in characteristic length of adjacent scales by the time interval. Based on the energy change rate and the scale change rate, a turbulence development index is obtained, and the turbulence development trend is judged based on the time series change characteristics of the turbulence development index.

[0048] An integrated microfluidic sensor array is deployed at multiple key cross-sectional locations within the conduit. Each cross-section contains multiple sensor nodes arranged in a ring pattern to collect pressure fluctuations, flow velocity distribution, and shear stress values ​​of the fluid within the conduit. The sensor node acquisition frequency is set to 2000Hz to meet the requirements for capturing high-frequency fluctuations in turbulent flow. The sensing unit employs a combination of a hot-film flow velocity sensor and a piezoresistive pressure sensor. Each sensor measures 100 μm × 100 μm and has a thickness of 5 μm, ensuring measurement accuracy without interfering with the fluid flow inside the conduit.

[0049] The multidimensional parameter feature matrix is ​​a reference database constructed based on a large amount of experimental data and numerical simulation results, containing characteristic parameters under various flow regimes. The mapping operation employs the tensor inner product method to calculate the strain rate tensor and vorticity tensor. The strain rate tensor describes the fluid deformation rate, while the vorticity tensor characterizes the fluid rotational properties. In a specific measurement, the velocity distribution detected at the third section of the duct exhibited significant radial non-uniformity, with a central velocity of 1.2 m / s and a velocity of only 0.3 m / s near the wall, indicating the presence of localized turbulence.

[0050] The second-order invariant of the vorticity tensor, which characterizes the vortex intensity, and the second-order invariant of the rate tensor, which characterizes the fluid deformation intensity, are calculated. Subtracting the second-order invariant of the rate tensor from the second-order invariant of the vorticity tensor yields the turbulence discrimination factor Q. When the Q value is greater than zero, it indicates that the vortex effect in that region exceeds the deformation effect, classifying it as a turbulent region. In practical applications, when the Q value exceeds a set threshold of 0.15, the region is marked as a locally turbulent region. In this example, the annular region at the radial position of the third section of the duct, from 0.7R to 0.85R (R is the duct radius), has a Q value of 0.23 and is identified as a locally turbulent region.

[0051] The velocity fluctuation component is extracted from the local turbulent region, which is the actual velocity minus the average velocity over a period of time. The root mean square value of the velocity fluctuation component is calculated, and its ratio to the average velocity of the region yields the turbulence intensity. Turbulence intensity characterizes the degree of fluid turbulence; a higher value indicates stronger turbulence. The warning intensity threshold is set at 0.15. When the turbulence intensity exceeds this threshold, the turbulence is deemed to have reached a level requiring a warning. In this detection, the calculated turbulence intensity of the local turbulent region was 0.18, exceeding the warning intensity threshold and triggering turbulence development trend analysis.

[0052] Velocity fluctuation signals in localized turbulent regions are extracted, and Morlet wavelet transform is applied to these signals. Morlet wavelets possess excellent time-frequency localization properties, making them suitable for analyzing non-stationary turbulent signals. The wavelet transform scale range is set from 1 to 64, with a total of 32 scale points, corresponding to physical scales from 0.5 mm to 25 mm. Wavelet coefficients are calculated at each scale, forming a wavelet coefficient matrix. The number of rows in this matrix represents the number of time sampling points, and the number of columns represents the number of scale points.

[0053] The local wave energy distribution is obtained by squaring each element in the wavelet coefficient matrix. The local wave energies are then summed at the same scale to obtain the energy values ​​at each scale. An energy-scale distribution plot is drawn, with scale on the horizontal axis and energy on the vertical axis. The energy distribution interval is extracted from the plot. In this example, the energy is mainly concentrated in the range of scales 6 to 18, corresponding to physical scales of approximately 3 to 9 millimeters, indicating that the turbulent structure is mainly concentrated in the small to medium scale range.

[0054] Within the energy distribution range, the ratio of the difference in turbulent energy between adjacent time points to the time interval is calculated to obtain the energy change rate. The time interval is set to 0.01 seconds. Simultaneously, the ratio of the difference in characteristic length between adjacent scales to 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 weighted and summed to calculate the turbulence development index, with weights of 0.7 and 0.3, respectively. In this example, the turbulence development indices at the initial 10 time points are 0.12, 0.15, 0.19, 0.24, 0.28, 0.33, 0.38, 0.42, 0.47, and 0.52, showing a clear upward trend.

[0055] When the index continues to increase and the growth rate exceeds 0.03 / second, it is judged as an increasing trend of turbulence; when the index fluctuates but remains generally stable within a range (fluctuation amplitude not exceeding ±0.05), it is judged as a stable trend of turbulence; when the index continues to decrease and the decrease rate exceeds 0.03 / second, it is judged as a decreasing trend of turbulence. According to the above judgment criteria, the turbulence in this case shows a significant increasing trend, with an average growth rate of 0.044 / second, generating a turbulence increasing warning signal, prompting operators to take corresponding measures to reduce the flow velocity or adjust the duct position to prevent the turbulence from developing further and causing potential risks.

[0056] In one optional implementation, based on the turbulence development trend and the frequency domain characteristics of pressure fluctuations, eliminating pressure fluctuations and controlling the pressure inside the duct to maintain within a safe pressure range includes: The pressure fluctuation signal inside the catheter is extracted, and the pressure fluctuation signal is decomposed into periodic fluctuation components in the time domain. The main frequency and amplitude of the fluctuation are obtained based on the Fourier series decomposition result of the periodic fluctuation components. Set a target value for catheter pressure, and use the difference between the real-time detected catheter pressure value and the target value for catheter pressure as the pressure deviation; Based on the turbulence development trend, the dominant frequency of the wave and its amplitude, calculate the proportional coefficient, integral time constant and derivative time constant. Multiply the proportional coefficient by the pressure deviation to obtain the proportional control term, multiply the integral time constant by the integral value of the pressure deviation to obtain the integral control term, multiply the derivative time constant by the rate of change of the pressure deviation to obtain the derivative control term, and add the proportional control term, integral control term and derivative control term to obtain the compensation control signal. The compensation control signal is input to the piezoelectric ceramic diaphragm pump to generate compensation pressure. The pressure value of the compensated conduit is detected in real time. When the pressure value of the compensated conduit exceeds the upper or lower safety limit, the rate of change of the pressure deviation is calculated. The compensation pressure is adjusted according to the rate of change of the pressure deviation to keep the pressure value in the conduit within the safe pressure range.

[0057] A miniature pressure sensor is installed inside the catheter to measure pressure changes in real time. The miniature pressure sensor has a range of -50 to 150 kPa, an accuracy of 0.01 kPa, and a sampling frequency of 1000 Hz. The acquired pressure fluctuation signal is represented as a time series P(t), where t represents the time variable.

[0058] The acquired pressure fluctuation signal P(t) is decomposed in the time domain, and wavelet analysis is used to perform multi-scale decomposition of P(t) to extract the periodic fluctuation component Pp(t). For example, when the values ​​of the pressure fluctuation signal P(t) in the conduit are 92.5 kPa, 93.1 kPa, 94.0 kPa, 93.5 kPa, 92.8 kPa, 92.2 kPa, 91.8 kPa, 92.4 kPa, 93.2 kPa, and 93.9 kPa within 10 seconds, the periodic fluctuation component Pp(t) obtained after wavelet analysis is 0.2 kPa, 0.8 kPa, 1.7 kPa, 1.2 kPa, 0.5 kPa, -0.1 kPa, -0.5 kPa, 0.1 kPa, 0.9 kPa, and 1.6 kPa.

[0059] The Fast Fourier Transform (FFT) algorithm is used to perform spectral analysis on Pp(t) to obtain its frequency components and corresponding amplitudes. In the example above, Fourier analysis shows that the dominant frequency f of the wave is 0.2 Hz, and the corresponding amplitude A is 1.2 kPa.

[0060] A target pressure value Ptarget is set for the catheter; in this embodiment, Ptarget is set to 92.3 kPa. The pressure value Pcurrent inside the catheter is monitored in real time, and the pressure deviation e = Pcurrent - Ptarget is calculated. For example, when the pressure value inside the catheter is 94.0 kPa at a certain moment, the pressure deviation e = 94.0 kPa - 92.3 kPa = 1.7 kPa.

[0061] The PID control parameters are calculated based on the turbulence development trend, the dominant frequency of the wave, and its amplitude. Specifically, the proportional gain Kp is inversely proportional to the wave amplitude A and directly proportional to the turbulence intensity; the integral time constant Ti is inversely proportional to the dominant frequency f; and the derivative time constant Td is directly proportional to the rate of change of the dominant frequency f. In this embodiment, when the turbulence intensity is 0.3, the dominant frequency f is 0.2Hz, and the amplitude A is 1.2kPa, the calculated proportional gain Kp = 0.25, integral time constant Ti = 10s, and derivative time constant Td = 0.05s.

[0062] Calculate the output values ​​of each term in the PID control. The proportional control term is Pout = Kp × e, the integral control term is Iout = e ×dt / Ti, and the derivative control term is 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 example above, when e = 1.7 kPa, the previous time e = 1.2 kPa, and dt = 0.01 s, the calculated values ​​are: 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.

[0063] Adding the proportional control term, integral control term, and derivative control term, we obtain the compensation control signal u = Pout + Iout + Dout. In the example above, u = 0.425kPa + 0.0017kPa + 2.5kPa = 2.9267kPa.

[0064] Input the compensation control signal u into the piezoelectric ceramic diaphragm pump to generate a compensation pressure. The output pressure range of the piezoelectric ceramic 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 ceramic diaphragm pump. In this embodiment, k = 1. In the above example, the compensation pressure Pcomp = 1 × 2.9267 kPa = 2.9267 kPa.

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

[0066] Set the pressure safety upper limit value Pupper = 95.0 kPa and the pressure safety lower limit value Plower = 90.0 kPa. Determine whether Pnew is within the safe pressure range, that is, whether Plower ≤ Pnew ≤ Pupper is satisfied. In the above example, 91.0733 kPa is between 90.0 kPa and 95.0 kPa, meeting the requirements of the safe pressure range.

[0067] When the compensated catheter pressure value exceeds the safe pressure range, calculate the rate of change of the pressure deviation de / dt = (e - eprev) / dt, where eprev is the pressure deviation at the previous moment. Dynamically adjust the compensation pressure value according to the magnitude and sign of de / dt. For example, when Pnew = 89.5 kPa < Plower, calculate de / dt = (e - eprev) / dt = (1.7 kPa - 1.2 kPa) / 0.01 s = 50 kPa / s. 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, adjust the compensation pressure 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, meeting the requirements of the safe pressure range.

[0068] Through the above method, the effective control of the pressure in the catheter is achieved, maintaining it within the safe pressure range and avoiding safety risks caused by too high or too low pressure.

[0069] This invention relates to a real-time monitoring and control system for intraductal pressure during ductoscopy, the system comprising: The first unit is used to collect real-time temperature data, conductivity data and image data inside the catheter during ductoscopy and form a multi-dimensional parameter feature matrix. The second unit is used to perform adaptive filtering analysis on the multidimensional parameter feature matrix, extract the feature parameters of submillimeter fluid particles in the fluid, calculate the spatial distribution density value and velocity vector of the submillimeter fluid particles based on the feature parameters, calculate the spatial distribution characteristics of the pressure inside the duct based on the spatial distribution density value, determine the location of pressure peak points and pressure valley points, and identify regions of abrupt pressure gradient changes; track the pressure wave propagation path based on the velocity vector, identify the reflection points and attenuation points of the pressure wave on the propagation path, and obtain a pressure wave propagation characteristic distribution map. The third unit is used to monitor the pressure inside the catheter in real time by combining the pressure gradient abrupt region and the pressure wave propagation characteristic distribution map. The fourth unit is used to monitor the fluid flow parameters of each cross section of the catheter through a microfluidic sensor array built into the catheter, identify local turbulent regions based on the fluid flow parameters and the multidimensional parameter feature matrix, calculate the turbulence development trend when the turbulence intensity of the local turbulence region exceeds the warning intensity threshold, and eliminate pressure fluctuations based on the turbulence development trend and the frequency domain characteristics of pressure fluctuations to control the pressure inside the catheter to be maintained within a safe pressure range.

[0070] A third aspect of the present invention provides an electronic device, comprising: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0071] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0072] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0073] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions 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 method for real-time monitoring and control of intraductal pressure during ductoscopy, characterized in that, include: Real-time temperature data, conductivity data, and image data inside the catheter were collected during ductoscopy and assembled into a multidimensional parameter feature matrix. Adaptive filtering analysis is performed on the multidimensional parameter feature matrix to extract the feature parameters of submillimeter fluid particles in the fluid. Based on the feature parameters, the spatial distribution density value and velocity vector of the submillimeter fluid particles are calculated. Based on the spatial distribution density value, the spatial distribution characteristics of the pressure inside the duct are calculated to determine the location of pressure peak points and pressure valley points, and to identify regions of abrupt pressure gradient changes. Based on the velocity vector, the propagation path of the pressure wave is traced to identify the reflection points and attenuation points of the pressure wave on the propagation path, thus obtaining a pressure wave propagation characteristic distribution map. By combining the distribution map of pressure gradient abrupt change regions and pressure wave propagation characteristics, the pressure inside the catheter can be monitored in real time; The fluid flow parameters of each cross section of the catheter are monitored by a microfluidic sensor array built into the catheter. Local turbulent regions are identified based on the fluid flow parameters and the multidimensional parameter feature matrix. When the turbulence intensity of the local turbulent region exceeds the warning intensity threshold, the turbulence development trend is calculated. Based on the turbulence development trend and the frequency domain characteristics of pressure fluctuations, pressure fluctuations are eliminated, and the pressure inside the duct is controlled to remain within a safe pressure range.

2. The method according to claim 1, characterized in that, Adaptive filtering analysis is performed on the multidimensional parameter feature matrix to extract the feature parameters of submillimeter fluid particles in the fluid. The spatial distribution density and velocity vector of the submillimeter fluid particles are calculated based on these feature parameters, including: The autocorrelation matrix of the multidimensional parameter feature matrix is ​​calculated by adaptive Wiener filtering to obtain the filtered feature data. In the filtered feature data, edge points are detected using a multi-scale operator, and the edge points are morphologically connected to obtain a complete contour. The equivalent diameter and roundness coefficient are calculated based on the complete contour, and the gray-level gradient value of the contour region is calculated at the same time. The equivalent diameter, roundness coefficient, and gray-level gradient value are combined to form the feature parameters of the submillimeter fluid particle. The initial bandwidth is determined based on the characteristic parameters. The initial bandwidth is then cross-validated and iteratively optimized until the density estimation error is less than the distribution fitting threshold, resulting in the optimized bandwidth parameters. The spatial distribution density of submillimeter fluid particles is then calculated using the optimized bandwidth parameters. Extract images of adjacent frames from the filtered feature data, calculate the spatial and temporal gradients of the adjacent frame images using the feature parameters, and construct a velocity field smoothing constraint. Use the spatial gradient, temporal gradient, and velocity field smoothing constraint to construct an optimization objective function. Iterate through the gradient descent method until the value of the optimization objective function is less than the velocity field convergence threshold, and determine the motion velocity vector of the submillimeter fluid particle.

3. The method according to claim 1, characterized in that, Based on the spatial distribution density value, calculate the spatial distribution characteristics of the pressure inside the catheter, determine the locations of pressure peak points and pressure trough points, and identify regions of abrupt pressure gradient changes, including: Calculate the initial pressure distribution based on the spatial density value within the catheter; The Laplacian operator for density is calculated, and a density gradient correction term is obtained based on the Laplacian operator and the density gradient coefficient. The density gradient correction term is then superimposed on the initial pressure distribution to obtain the corrected pressure distribution. Multi-scale analysis is performed on the corrected pressure distribution, 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. Local extreme points where the gradient vector field shows convergence are marked as pressure peak points, and local extreme points where the gradient vector field shows divergence are marked as pressure valley points. Based on the distribution locations of the pressure peak points and pressure valley points, the gradient amplitude of the corrected pressure distribution at each location is calculated, and the gradient direction angle between adjacent locations is calculated to obtain the direction change rate. The gradient amplitude is divided by the average gradient amplitude in the local area to obtain the amplitude ratio. The gradient mutation index is calculated based on the amplitude ratio and the direction change rate. When the gradient mutation index is greater than the mutation judgment threshold, the area between adjacent pressure peak points and pressure valley points is marked as a pressure gradient mutation region.

4. The method according to claim 1, characterized in that, Based on the motion velocity vector, the propagation path of the pressure wave is traced, and the reflection points and attenuation points of the pressure wave on the propagation path are identified to obtain the pressure wave propagation characteristic distribution map, including: Extract the velocity components of discrete points in the corrected pressure distribution, calculate the pressure wave acceleration based on the velocity components of adjacent discrete points, and integrate the pressure wave acceleration over time to obtain the pressure wave propagation path. The inner product of the velocity direction at each moment and the velocity direction at the next moment is calculated along the propagation path of the pressure wave. The propagation direction change angle is obtained based on the inner product and the velocity amplitude. When the propagation direction change angle is greater than the direction change threshold, the corresponding coordinate point is marked as the reflection point. 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. The pressure gradient along the path is divided by the pressure value to obtain the attenuation coefficient. When the attenuation coefficient is greater than the attenuation characteristic threshold, the corresponding coordinate point is marked as an attenuation point. Calculate the distance from each reflection point to each attenuation point, pair reflection points and attenuation points whose distance is less than the feature association threshold to form a feature point group, and calculate the spatial distribution density of the feature point group to obtain the propagation feature intensity. Draw the propagation path of the pressure wave, mark the propagation characteristic intensity, mark the location and direction of the reflection point, mark the location and attenuation coefficient of the attenuation point, and generate a distribution map of the propagation characteristics of the pressure wave.

5. The method according to claim 1, characterized in that, Fluid flow parameters at various cross-sections of the catheter are monitored using a microfluidic sensor array embedded within the catheter. Local turbulent regions are identified based on these fluid flow parameters and the multidimensional parameter feature matrix. When the turbulence intensity in a local turbulent region exceeds a warning intensity threshold, the turbulence development trend is calculated, including: The pressure fluctuation values ​​and flow velocity distribution values ​​of each section are collected and mapped with the multidimensional parameter feature matrix to obtain the strain rate tensor and vorticity tensor. The second-order invariants of the vorticity tensor and the second-order invariants of the strain rate tensor are calculated and subtracted to obtain the turbulence discrimination factor. The local turbulence region is marked based on the turbulence discrimination factor. Extract the velocity fluctuation component in the local turbulent region, and calculate the ratio of the root mean square value of the velocity fluctuation component to the average velocity to obtain the turbulence intensity. When the turbulence intensity exceeds the warning intensity threshold, the velocity fluctuation signal of the local turbulence region is extracted and wavelet transform is performed. The wavelet coefficient matrix is ​​obtained by integrating at different scales. The square of each element in the wavelet coefficient matrix is ​​calculated to obtain the local wave energy. The local wave energy is summed at the same scale to obtain the scale energy. The energy distribution interval is extracted based on the scale energy. Within the energy distribution range, the energy change rate is obtained by dividing the difference in turbulent energy at adjacent times by the time interval, and the scale change rate is obtained by dividing the difference in characteristic length of adjacent scales by the time interval. Based on the energy change rate and the scale change rate, a turbulence development index is obtained, and the turbulence development trend is judged based on the time series change characteristics of the turbulence development index.

6. The method according to claim 1, characterized in that, Based on the aforementioned turbulence development trend and the frequency domain characteristics of pressure fluctuations, eliminating pressure fluctuations and controlling the pressure inside the duct to maintain within a safe pressure range includes: The pressure fluctuation signal inside the catheter is extracted, and the pressure fluctuation signal is decomposed into periodic fluctuation components in the time domain. The main frequency and amplitude of the fluctuation are obtained based on the Fourier series decomposition result of the periodic fluctuation components. Set a target value for catheter pressure, and use the difference between the real-time detected catheter pressure value and the target value for catheter pressure as the pressure deviation; Based on the turbulence development trend, the dominant frequency of the wave and its amplitude, calculate the proportional coefficient, integral time constant and derivative time constant. Multiply the proportional coefficient by the pressure deviation to obtain the proportional control term, multiply the integral time constant by the integral value of the pressure deviation to obtain the integral control term, multiply the derivative time constant by the rate of change of the pressure deviation to obtain the derivative control term, and add the proportional control term, integral control term and derivative control term to obtain the compensation control signal. The compensation control signal is input to the piezoelectric ceramic diaphragm pump to generate compensation pressure. The pressure value of the compensated conduit is detected in real time. When the pressure value of the compensated conduit exceeds the upper or lower safety limit, the rate of change of the pressure deviation is calculated. The compensation pressure is adjusted according to the rate of change of the pressure deviation to keep the pressure value in the conduit within the safe pressure range.

7. A real-time monitoring and control system for ductal pressure during ductoscopy, used to implement the method as described in any one of claims 1-6, characterized in that, include: The first unit is used to collect real-time temperature data, conductivity data and image data inside the catheter during ductoscopy and form a multi-dimensional parameter feature matrix. The second unit is used to perform adaptive filtering analysis on the multidimensional parameter feature matrix, extract the feature parameters of submillimeter fluid particles in the fluid, calculate the spatial distribution density value and velocity vector of the submillimeter fluid particles based on the feature parameters, calculate the spatial distribution characteristics of the pressure inside the duct based on the spatial distribution density value, determine the location of pressure peak points and pressure valley points, and identify regions of abrupt pressure gradient changes; track the pressure wave propagation path based on the velocity vector, identify the reflection points and attenuation points of the pressure wave on the propagation path, and obtain a pressure wave propagation characteristic distribution map. The third unit is used to monitor the pressure inside the catheter in real time by combining the pressure gradient abrupt region and the pressure wave propagation characteristic distribution map. The fourth unit is used to monitor the fluid flow parameters of each cross section of the catheter through a microfluidic sensor array built into the catheter, identify local turbulent regions based on the fluid flow parameters and the multidimensional parameter feature matrix, and calculate the turbulence development trend when the turbulence intensity of the local turbulent region exceeds the warning intensity threshold. Based on the turbulence development trend and the frequency domain characteristics of pressure fluctuations, pressure fluctuations are eliminated, and the pressure inside the duct is controlled to remain within a safe pressure range.

8. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 6.

9. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 6.

Images