Transfusion dripping speed abnormity early warning adaptive algorithm method and system

Abstract

The invention provides an infusion dripping speed abnormity early warning self-adaptive algorithm method and system, and relates to the technical field of medical monitoring, and the method comprises the steps: obtaining a dripping speed data flow, constructing a phase space trajectory matrix, carrying out the singular value decomposition, determining a mode demarcation point, separating a dominant mode from a disturbance mode, and calculating a dripping speed deviation measurement value; and extracting a topological feature vector by using the recurrence plot matrix to judge abnormity, and determining an early warning level based on a feature time scale and an isolation degree change rate. According to the invention, early recognition and early warning of the abnormal dripping speed in the infusion process can be realized, and the infusion safety is improved.

Description

Technical Field

[0001] This invention relates to the field of medical monitoring technology, and in particular to an adaptive algorithm method and system for early warning of abnormal infusion drip rate. Background Technology

[0002] Intravenous infusion is a common method of drug administration in clinical medicine, and accurate control of the infusion drip rate is crucial to ensuring patient medication safety. Traditionally, medical staff monitor the infusion drip rate through regular rounds, but this method suffers from long monitoring intervals and poor real-time performance. With the development of medical information technology, intelligent infusion monitoring systems are gradually being applied in clinical practice. These systems use sensors to collect infusion drip rate data in real time and utilize algorithms to provide early warnings of abnormal drip rates.

[0003] Existing technologies for monitoring abnormal infusion drip rates mainly employ fixed threshold judgment, statistical analysis, and simple machine learning methods. Fixed threshold judgment uses preset upper and lower limits for drip rate to identify anomalies; statistical analysis detects anomalies based on the statistical characteristics of drip rate data, such as mean and variance; and machine learning methods identify abnormal drip rate patterns by training models.

[0004] However, existing technologies still have some shortcomings and deficiencies. Most existing methods for detecting abnormal drip rates use static threshold settings, which are difficult to adapt to individual differences among patients and the characteristics of drip rate changes at different stages of infusion, resulting in a high false positive rate and affecting the accuracy of clinical decision-making. Traditional algorithms have difficulty effectively distinguishing between transient disturbances and persistent abnormalities during infusion, lack the ability to deeply analyze the dynamic trend of drip rate changes, and cannot predict possible infusion abnormalities in advance, only passively identifying them after they occur. Existing technologies usually only focus on the change of the drip rate value itself, ignoring the temporal structure information and nonlinear dynamic characteristics contained in the drip rate sequence, making it difficult to capture the complex state transition process in the infusion system, thus limiting the accuracy and foresight of abnormality warnings. Summary of the Invention

[0005] This invention provides an adaptive algorithm method and system for early warning of abnormal infusion drip rate, which can solve the problems in the prior art.

[0006] A first aspect of this invention provides an adaptive algorithm method for early warning of abnormal infusion drip rate, comprising: Acquire a drip rate monitoring data stream during the infusion process, wherein the drip rate monitoring data stream includes time-series drip rate values ​​and corresponding time identifiers; A phase space trajectory matrix is ​​constructed for the time series drop rate values ​​within the historical time window, and singular value decomposition is performed to obtain the singular value spectrum. The mode boundary point is determined by the bidirectional decay gradient field, and the singular value spectrum is decomposed into the dominant mode and the perturbation mode. The drop rate center baseline is reconstructed and the drop rate fluctuation boundary is determined. The difference between the real-time drip rate value and the drip rate center baseline is calculated and compared with the drip rate fluctuation boundary to obtain the drip rate deviation metric value. The drip rate deviation metric values ​​within a continuous time window are organized into a drip rate deviation sequence. The phase space of the drop rate deviation metric sequence is reconstructed and a recursive graph matrix is ​​generated. Based on the recursive graph matrix, a topological connectivity graph is constructed and recursive topological feature vectors are extracted. The isolation degree is determined by the decay rate of the recursive structure stability index. When the isolation degree exceeds the preset isolation threshold, the anomaly identification result is output. For the sequence of drop rate deviation measurements before the corresponding time point of the anomaly identification result, calculate the state transition information transmission amount at different time intervals, determine the characteristic time scale, calculate the comprehensive anomaly index based on the frequency domain features of the characteristic time scale and the isolation degree change rate, determine the warning level and generate warning information.

[0007] In one optional embodiment, a phase space trajectory matrix is ​​constructed for the time-series drop rate values ​​within a historical time window, and singular value decomposition is performed to obtain a singular value spectrum. The mode boundary point is determined through a bidirectional decaying gradient field. The singular value spectrum is decomposed into dominant modes and perturbation modes. Reconstructing the drop rate center baseline and determining the drop rate fluctuation boundary includes: The mutual information function under different time delays is calculated for the time-series drop rate values ​​within the historical time window. The time delay at which the mutual information function first reaches its minimum value is identified, and the embedding delay time is determined. The pseudo-nearest neighbor ratio is calculated for the time-series drop rate values ​​under different embedding dimensions. The embedding dimension at which the pseudo-nearest neighbor ratio drops below a preset stable threshold is identified, and the phase space embedding dimension is determined. Based on the embedding delay time and the phase space embedding dimension, the time-series drop rate values ​​are expanded according to the time sliding window to generate the phase space trajectory matrix. Singular value decomposition is performed on the phase space trajectory matrix to obtain the singular value spectrum arranged in descending order of energy and the corresponding left and right singular vector sets; A bidirectional decaying gradient field is constructed on the singular value spectrum to obtain a decaying acceleration sequence. Based on the decaying acceleration sequence, a candidate mode boundary set is extracted, and the mode separation index is calculated to determine the mode boundary point. Singular values ​​and their corresponding left and right singular vectors before the mode boundary are classified as dominant modes, while singular values ​​and their corresponding left and right singular vectors after the mode boundary are classified as perturbation modes. The drop rate center baseline is reconstructed by performing matrix multiplication on the singular values ​​and left and right singular vectors in the dominant mode and projecting them back to the original time dimension. The mean and standard deviation of the singular values ​​in the perturbation mode are calculated, and confidence intervals are constructed to determine the drop rate fluctuation boundary.

[0008] In one optional embodiment, constructing a bidirectional decaying gradient field on the singular value spectrum to obtain a decaying acceleration sequence, extracting a set of candidate mode boundaries based on the decaying acceleration sequence, and calculating the mode separation index to determine the mode boundary points include: A two-way decay gradient field is constructed for the singular value spectrum. The forward energy decay gradient and the backward energy decay gradient are calculated at each singular value position. The difference between the forward energy decay gradient and the backward energy decay gradient is used to obtain the gradient symmetry deviation sequence. The decay acceleration sequence is obtained by stepwise recursive difference. An adaptive detection threshold is constructed based on the statistical distribution of the decay acceleration sequence. The positions of singular values ​​in the decay acceleration sequence that exceed the adaptive detection threshold are extracted to form a set of candidate mode boundaries. For each candidate modal boundary in the candidate modal boundary set, an energy segmentation topology is constructed. The energy distribution entropy difference between the front and back subspaces of the candidate modal boundary is calculated. The energy distribution entropy difference is coupled with the energy proportion variance to obtain the modal separation index. The candidate modal boundary with the largest modal separation index is selected to determine the modal boundary point.

[0009] In one optional embodiment, the drop rate deviation metric sequence is reconstructed in phase space to generate a recursive graph matrix. A topological connectivity graph is constructed based on the recursive graph matrix, and recursive topological feature vectors are extracted. The isolation degree is determined by the decay rate of the recursive structure stability index, including: The phase space of the drip rate deviation measurement sequence is reconstructed, the phase space distance between each time point in the reconstructed phase space and the historical time point is calculated, and the effective recursion points are determined based on the dual screening of adaptive recursion radius threshold and trajectory similarity, and a recursion graph matrix is ​​generated. A topological connectivity graph is constructed based on the distribution of recursive points at each time point in the recursive graph matrix. Clusters and isolated nodes are identified by graph clustering algorithm. The cohesion of clusters is obtained by calculating the connection density of clusters, and the separation degree is obtained by calculating the shortest path length between isolated nodes and the nearest cluster. A recursive topological feature vector is constructed based on the cohesion and separation degree. A temporal sliding window extraction is performed on the recursive topological feature vectors, and the temporal interval entropy and spatial distribution entropy of the recursive points within the window are calculated to construct a recursive structure stability index. The decay rate sequence on the time axis of the recursive structure stability index is calculated, the moment when the decay rate exceeds the preset mutation threshold is identified, and the product of the decay amplitude and decay rate at that moment is extracted to obtain the isolation degree.

[0010] In one optional embodiment, the drop rate deviation metric sequence is reconstructed in phase space. The phase space distance between each time point in the reconstructed phase space and the historical time points is calculated. Valid recursive points are determined based on a dual screening of adaptive recursion radius threshold and trajectory similarity. The resulting recursion graph matrix includes: The phase space of the drop rate deviation metric sequence is reconstructed. In the reconstructed phase space, a time backtracking window is constructed for each time point. Historical time points are extracted within the time backtracking window, and the Euclidean distance between the time point and the historical time point is calculated to obtain the phase space distance. The sequence of drip rate deviation metrics is segmented and statistically analyzed according to time windows. The standard deviation of the drip rate deviation metrics within each time window is calculated, and the product of the standard deviation and the preset sensitivity coefficient is used as the adaptive recursive radius threshold. Historical time points with phase space distances less than the adaptive recursive radius threshold are selected to construct a candidate recursive point set. For each candidate recursive point, the corresponding subsequent evolution fragment is extracted and a point-to-point distance matrix is ​​formed with the subsequent evolution fragment of the current time point. Based on the distance matrix, dynamic programming is performed to calculate the cumulative distance matrix. The alignment path that minimizes the total distance is extracted from the cumulative distance matrix. The cumulative distance value of the endpoint of the alignment path is extracted to determine the dynamic time warping distance. The dynamic time warping distance is transformed by the inverse to obtain the trajectory similarity. Candidate recursive points whose trajectory similarity exceeds a preset trajectory similarity threshold are marked as valid recursive points. Based on the distribution of valid recursive points on the time axis, a recursion graph matrix is ​​generated.

[0011] In one optional embodiment, for the drop rate deviation measurement sequence before the corresponding time of the anomaly identification result, the state transition information transmission amount at different time intervals is calculated to determine the characteristic time scale. Based on the frequency domain features of the characteristic time scale and the isolation degree change rate, a comprehensive anomaly index is calculated to determine the warning level and generate warning information, including: Discretize the sequence of drop rate deviation metrics before the time corresponding to the anomaly identification result, calculate the transition probability between different state values ​​in the sequence of drop rate deviation metrics, calculate the information transmission amount of state transition according to different time intervals, identify the time interval with the largest information transmission amount, and determine the characteristic time scale. The isolation degree change rate sequence is obtained by calculating the difference between adjacent time points of the isolation degree sequence. The frequency domain spectrum is obtained by performing a Fourier transform on the isolation degree change rate sequence. The frequency component with the largest amplitude is identified from the frequency domain spectrum as the dominant frequency. The synchronization deviation degree is obtained by calculating the deviation ratio between the period corresponding to the dominant frequency and the characteristic time scale. The energy values ​​of the frequency components other than the dominant frequency in the frequency domain spectrum are accumulated. The ratio of the accumulated energy value to the energy value of the dominant frequency is used to obtain the spectral dispersion. The synchronization deviation and the spectral dispersion are weighted and summed to obtain the comprehensive anomaly index. The warning level is determined according to the numerical range of the comprehensive anomaly index. The warning level is encapsulated with the corresponding time of the anomaly identification result to generate warning information.

[0012] In one optional embodiment, the discretization process includes: A state space is constructed for the drop rate deviation metric sequence. Attractor centers with numerical clusters are identified in the state space. An attraction domain boundary is formed by expanding outward from each attractor center. The attraction strength from any point in the state space to each attractor center is calculated. The state label corresponding to the attractor with the largest attraction strength is assigned to the corresponding point. Based on the attractor affiliation relationship, the drop rate deviation metric sequence is mapped to a discrete state sequence.

[0013] A second aspect of this invention provides an adaptive algorithm system for early warning of abnormal infusion drip rate, comprising: The data acquisition unit is used to acquire the drip rate monitoring data stream during the infusion process. The drip rate monitoring data stream includes time-series drip rate values ​​and corresponding time identifiers. The baseline reconstruction unit is used to construct the phase space trajectory matrix of the time-series drop rate values ​​within the historical time window, and to perform singular value decomposition to obtain the singular value spectrum. The mode boundary point is determined by the bidirectional decay gradient field, and the singular value spectrum is decomposed into the dominant mode and the perturbation mode. The drop rate center baseline is reconstructed and the drop rate fluctuation boundary is determined. The deviation calculation unit is used to calculate the difference between the real-time drip rate value and the drip rate center baseline, and to perform a ratio operation with the drip rate fluctuation boundary to obtain the drip rate deviation metric value. The drip rate deviation metric values ​​within a continuous time window are organized into a drip rate deviation sequence. An anomaly identification unit is used to reconstruct the phase space of the drop rate deviation measurement value sequence and generate a recursive graph matrix. Based on the recursive graph matrix, a topological connectivity graph is constructed and recursive topological feature vectors are extracted. The isolation degree is determined by the decay rate of the recursive structure stability index. When the isolation degree exceeds the preset isolation threshold, the anomaly identification result is output. The early warning generation unit is used to calculate the state transition information transmission amount at different time intervals from the drop rate deviation measurement value sequence before the corresponding time of the anomaly identification result, determine the characteristic time scale, calculate the comprehensive anomaly index based on the frequency domain features of the characteristic time scale and the isolation degree change rate, determine the early warning level, and generate early warning information.

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

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

[0016] In this embodiment of the invention, by constructing a phase space trajectory matrix from the time-series drip rate values ​​within a historical time window and performing singular value decomposition, combined with a bidirectional decay gradient field to determine the mode boundary point, the dominant mode and the perturbation mode of the drip rate data are accurately separated. This enables more accurate reconstruction of the drip rate center baseline and determination of the drip rate fluctuation boundary, improving the sensitivity to drip rate anomalies. A method for calculating the drip rate deviation metric is introduced, quantifying the degree of drip rate deviation by calculating the ratio of the difference between the real-time drip rate value and the center baseline to the fluctuation boundary. By employing phase space reconstruction and a recursive graph matrix, combined with topological connectivity graph extraction of recursive topological feature vectors, the nonlinear dynamic characteristics of drip rate anomalies can be effectively captured. The isolation degree is determined by the decay rate of the recursive structure stability index, achieving high-accuracy identification of drip rate anomalies. A method for calculating a comprehensive anomaly index based on the frequency domain features of the characteristic time scale and the rate of change of isolation degree enables adaptive determination of the anomaly warning level, making the warning information more accurate and reliable, and significantly improving the safety of the infusion process and the work efficiency of medical staff. Attached Figure Description

[0017] Figure 1 This is a flowchart illustrating the adaptive algorithm method for early warning of abnormal infusion drip rate according to an embodiment of the present invention. Figure 2 This is a flowchart of the abnormal warning logic. Detailed Implementation

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

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

[0020] Figure 1 This is a flowchart illustrating the adaptive algorithm method for early warning of abnormal infusion drip rate according to an embodiment of the present invention. Figure 1 As shown, the method includes: Acquire a drip rate monitoring data stream during the infusion process, wherein the drip rate monitoring data stream includes time-series drip rate values ​​and corresponding time identifiers; A phase space trajectory matrix is ​​constructed for the time series drop rate values ​​within the historical time window, and singular value decomposition is performed to obtain the singular value spectrum. The mode boundary point is determined by the bidirectional decay gradient field, and the singular value spectrum is decomposed into the dominant mode and the perturbation mode. The drop rate center baseline is reconstructed and the drop rate fluctuation boundary is determined. The difference between the real-time drip rate value and the drip rate center baseline is calculated and compared with the drip rate fluctuation boundary to obtain the drip rate deviation metric value. The drip rate deviation metric values ​​within a continuous time window are organized into a drip rate deviation sequence. The phase space of the drop rate deviation metric sequence is reconstructed and a recursive graph matrix is ​​generated. Based on the recursive graph matrix, a topological connectivity graph is constructed and recursive topological feature vectors are extracted. The isolation degree is determined by the decay rate of the recursive structure stability index. When the isolation degree exceeds the preset isolation threshold, the anomaly identification result is output. For the sequence of drop rate deviation measurements before the corresponding time point of the anomaly identification result, calculate the state transition information transmission amount at different time intervals, determine the characteristic time scale, calculate the comprehensive anomaly index based on the frequency domain features of the characteristic time scale and the isolation degree change rate, determine the warning level and generate warning information.

[0021] In one optional implementation, a phase space trajectory matrix is ​​constructed for the time-series drop rate values ​​within a historical time window, and singular value decomposition is performed to obtain a singular value spectrum. The mode boundary point is determined using a bidirectional decaying gradient field. The singular value spectrum is decomposed into dominant modes and perturbation modes. Reconstructing the drop rate center baseline and determining the drop rate fluctuation boundary includes: The mutual information function under different time delays is calculated for the time-series drop rate values ​​within the historical time window. The time delay at which the mutual information function first reaches its minimum value is identified, and the embedding delay time is determined. The pseudo-nearest neighbor ratio is calculated for the time-series drop rate values ​​under different embedding dimensions. The embedding dimension at which the pseudo-nearest neighbor ratio drops below a preset stable threshold is identified, and the phase space embedding dimension is determined. Based on the embedding delay time and the phase space embedding dimension, the time-series drop rate values ​​are expanded according to the time sliding window to generate the phase space trajectory matrix. Singular value decomposition is performed on the phase space trajectory matrix to obtain the singular value spectrum arranged in descending order of energy and the corresponding left and right singular vector sets; A bidirectional decaying gradient field is constructed on the singular value spectrum to obtain a decaying acceleration sequence. Based on the decaying acceleration sequence, a candidate mode boundary set is extracted, and the mode separation index is calculated to determine the mode boundary point. Singular values ​​and their corresponding left and right singular vectors before the mode boundary are classified as dominant modes, while singular values ​​and their corresponding left and right singular vectors after the mode boundary are classified as perturbation modes. The drop rate center baseline is reconstructed by performing matrix multiplication on the singular values ​​and left and right singular vectors in the dominant mode and projecting them back to the original time dimension. The mean and standard deviation of the singular values ​​in the perturbation mode are calculated, and confidence intervals are constructed to determine the drop rate fluctuation boundary.

[0022] In one specific implementation, parameters are determined before phase space reconstruction of the time-series drop rate values ​​collected within a historical time window. A drop rate data sequence of appropriate length is selected, such as drop rate records for the past six hours, denoted as X(t), t=1, 2, ..., N. The mutual information function I(τ) under different time delays τ is calculated, specifically by constructing the joint probability distribution and marginal probability distribution of the drop rate values ​​at time t and t+τ, and then calculating their information entropy difference. The mutual information function decreases first and then oscillates as τ increases. When I(τ) reaches its first local minimum, the corresponding τ value is determined as the optimal embedding delay time. For example, when τ=8, I(τ) first reaches its minimum value, so the embedding delay time is set to 8.

[0023] The phase space embedding dimension *d* is determined, and a pseudo-nearest neighbor algorithm is used to calculate the pseudo-nearest neighbor ratio *P(d)* for each possible embedding dimension *d*. Specifically, for each point in the sequence, its nearest neighbor is found in the *d*-dimensional space, and then it is checked whether this pair of points still maintains a proximity relationship in the *d+1*-dimensional space. When the pseudo-nearest neighbor ratio decreases to below a preset threshold (e.g., 0.01) after increasing the dimension, the corresponding dimension is the optimal embedding dimension. Assuming that when *d*=5, *P(d)*<0.01, then the phase space embedding dimension is determined to be 5.

[0024] Based on a defined embedding delay time τ and embedding dimension d, a phase space trajectory matrix A is constructed. Each row of matrix A represents a point in the phase space, in the form [X(i), X(i+τ), X(i+2τ), ..., X(i+(d-1)τ)], where i=1, 2, ..., N-(d-1)τ. For example, if τ=8 and d=5, then the first row of A is [X(1), X(9), X(17), X(25), X(33)], the second row is [X(2), X(10), X(18), X(26), X(34)], and so on.

[0025] Singular value decomposition is performed on the constructed phase space trajectory matrix A. The decomposition yields A = UΣV T Σ is a diagonal matrix, and the elements σ1≥σ2≥...≥σd≥0 on the diagonal are singular values, arranged in descending order of energy; U and V are the left and right singular vector matrices, respectively.

[0026] To determine the mode boundary, the bidirectional decaying gradient field of the singular value spectrum is calculated, and the difference sequence between adjacent singular values, ΔΣ(i) = σi - σi+1, i = 1, 2, ..., d-1, is calculated. The rate of change of the difference sequence, i.e., the decaying acceleration sequence, ΔΔΣ(i) = ΔΣ(i) - ΔΣ(i+1), i = 1, 2, ..., d-2, is calculated. Significant fluctuation points in the decaying acceleration sequence are identified as candidate mode boundaries. By calculating the difference index (such as entropy difference) of the singular value distribution before and after each candidate point, the point with the largest difference is selected as the mode boundary point k. For example, if the decaying acceleration sequence changes significantly at point 3, and the entropy difference is the largest at that point, then k = 3 is determined as the mode boundary point.

[0027] After determining the mode boundary points, the singular value spectrum is decomposed into dominant modes and perturbation modes. The dominant mode contains the first k singular values ​​and their corresponding left and right singular vectors, denoted as Σ1=[σ1, σ2, ..., σ...]. k ],U1=[u1,u2,...,u k V1 = [v1, v2, ..., v] k The perturbation mode contains residual singular values ​​and a vector, represented as Σ2=[σ]. k+1 , σ k+2 , ..., σ d ],U2=[u k+1 u k+2 , ..., u d ], V2=[v k+1 v k+2 , ..., v d ].

[0028] Reconstruct the drop rate center baseline using the dominant mode. This is achieved through matrix multiplication A1=U1Σ1V1. T The phase space trajectory matrix after noise filtering is obtained. A1 is projected back to the original time dimension by calculating the average value of the positions of all trajectory points containing X(t) at each time point t, thus obtaining the reconstructed drop velocity center baseline X'(t).

[0029] Determine the drip rate fluctuation boundary and calculate the mean μ and standard deviation σ of the singular values ​​of the perturbation mode. Based on the normal distribution assumption, construct the confidence interval [X'(t)-zα / 2·σ, X'(t)+zα / 2·σ] as the drip rate fluctuation boundary, where zα / 2 is the critical value of the standard normal distribution corresponding to the significance level α, usually zα / 2=1.96 corresponding to the 95% confidence level.

[0030] For example, drip rate data was collected for four hours (4 hours) from a patient receiving controlled intravenous infusion, totaling 240 data points. The embedding delay τ=6 and embedding dimension d=4 were determined using the method described above, and the modal boundary point k=2 was calculated. The reconstructed drip rate center baseline can smoothly filter out sudden disturbances, while the fluctuation boundary can identify abnormal drip rate changes exceeding the normal range, providing healthcare professionals with accurate monitoring data.

[0031] In this embodiment, deterministic dynamic behavior and random disturbances in drip rate time series data can be effectively distinguished, improving the accuracy of drip rate monitoring and providing a more refined control strategy for the intravenous infusion process.

[0032] In one optional implementation, a bidirectional decaying gradient field is constructed on the singular value spectrum to obtain a decaying acceleration sequence. Based on the decaying acceleration sequence, a set of candidate mode boundaries is extracted, and the mode separation index is calculated to determine the mode boundary points, including: A two-way decay gradient field is constructed for the singular value spectrum. The forward energy decay gradient and the backward energy decay gradient are calculated at each singular value position. The difference between the forward energy decay gradient and the backward energy decay gradient is used to obtain the gradient symmetry deviation sequence. The decay acceleration sequence is obtained by stepwise recursive difference. An adaptive detection threshold is constructed based on the statistical distribution of the decay acceleration sequence. The positions of singular values ​​in the decay acceleration sequence that exceed the adaptive detection threshold are extracted to form a set of candidate mode boundaries. For each candidate modal boundary in the candidate modal boundary set, an energy segmentation topology is constructed. The energy distribution entropy difference between the front and back subspaces of the candidate modal boundary is calculated. The energy distribution entropy difference is coupled with the energy proportion variance to obtain the modal separation index. The candidate modal boundary with the largest modal separation index is selected to determine the modal boundary point.

[0033] In one specific implementation, a bidirectional decaying gradient field is constructed on the singular value spectrum. The singular value spectrum sequence S=[s1, s2, ..., s...] of the signal to be analyzed is obtained. n ], where n is the total number of singular values, s1≥s2≥...≥s n For each singular value position i, calculate the forward energy decay gradient G. forward(i) and backward energy decay gradient G backward(i) The forward energy decay gradient represents the rate of energy change when looking forward from the current position, while the backward energy decay gradient represents the rate of energy change when looking backward from the current position.

[0034] The formula for calculating the forward energy decay gradient is: G forward(i = (s i - s i-w) / w, where w is the width of the gradient calculation window, usually taking values from 3 to 5. When i < w, boundary handling methods are adopted to take the nearest available singular value.

[0035] The formula for the backward energy decay gradient is: G backward(i) = (s i - s i+w ) / w. When i + w > n, boundary handling methods are also adopted.

[0036] Calculate the gradient symmetry deviation sequence D = [d1, d2,..., d n , where d i = G forward(i) - G backward(i) . The gradient symmetry deviation sequence reflects the asymmetry of the singular value decay at each position. The larger the value, the more significant the change in the singular value decay characteristics at that position.

[0037] Perform successive order recursive differences on the gradient symmetry deviation sequence to obtain the decay acceleration sequence A = [a1, a2,..., a n-1 , where a i = d i+1 - d i . The decay acceleration sequence characterizes the change in the singular value decay rate, and significant peaks usually appear at the modal boundaries.

[0038] Construct an adaptive detection threshold based on the statistical distribution of the decay acceleration sequence, and calculate the mean μ and standard deviation σ of the decay acceleration sequence A. Then set the adaptive detection threshold T = μ + k×σ, where k is an adjustment coefficient, which can take values from 2 to 3 according to the actual application scenario. A smaller k value will generate more candidate boundary points, with high sensitivity but possibly introducing false boundaries; a larger k value is more conservative and only retains significant boundary points.

[0039] Extract the positions of the singular values in the decay acceleration sequence that exceed the adaptive detection threshold to form a candidate modal boundary set B = {b1, b2,..., b m}, where for any b j , it satisfies a_{b j} > T.

[0040] Evaluate each candidate point in the candidate modal boundary set to determine the final modal boundary point. For each candidate boundary b j , construct an energy segmentation topological structure to divide the singular value spectrum into a front subspace S front = {s1, s2,..., s_{b j}} and a rear subspace S back = {s_{bj+1},s_{b j+2},...,s n}

[0041] Calculate the energy distribution entropy H of the front subspace. front and the energy distribution entropy H of the back subspace back Entropy of energy distribution reflects the uniformity of the distribution of singular values ​​within a subspace. During calculation, the singular values ​​within the subspace are normalized, and the entropy definition formula is applied.

[0042] Calculate the entropy difference in energy distribution ΔH = |H front - H back The larger the entropy difference, the more significant the difference in energy distribution between the two subspaces, and the more likely it is to be a true modal boundary.

[0043] Simultaneously, calculate the energy percentage variance V. First, calculate the energy percentage P of the front subspace. front = sum(S front The energy percentage P of the back subspace is calculated as ) / sum(S). back = sum(S back ) / sum(S). Then calculate the energy percentage variance V = (P front - 0.5) 2 + (P back - 0.5) 2 The variance of energy proportion reflects the balance of energy distribution between the two subspaces; the smaller the value, the more balanced the division.

[0044] By coupling the entropy difference of energy distribution with the variance of energy proportion, the modal separation index M = ΔH ×(1 - V) is obtained. This index comprehensively considers the distribution differences and energy balance between subspaces, requiring significant differences between subspaces while also favoring a relatively balanced energy partition.

[0045] The candidate mode boundary with the largest modal separation index is selected as the final modal boundary point.

[0046] In this embodiment, by constructing a bidirectional decay gradient field and decay acceleration sequence, combined with adaptive threshold detection and multi-index comprehensive evaluation, the mode boundary points in the singular value spectrum can be automatically and accurately identified, providing strong support for signal processing, pattern recognition and other fields.

[0047] In one optional implementation, the drop rate deviation metric sequence is reconstructed in phase space to generate a recursive graph matrix. A topological connectivity graph is constructed based on the recursive graph matrix, and recursive topological feature vectors are extracted. The isolation degree is determined by the decay rate of the recursive structure stability index, including: The phase space of the drip rate deviation measurement sequence is reconstructed, the phase space distance between each time point in the reconstructed phase space and the historical time point is calculated, and the effective recursion points are determined based on the dual screening of adaptive recursion radius threshold and trajectory similarity, and a recursion graph matrix is ​​generated. A topological connectivity graph is constructed based on the distribution of recursive points at each time point in the recursive graph matrix. Clusters and isolated nodes are identified by graph clustering algorithm. The cohesion of clusters is obtained by calculating the connection density of clusters, and the separation degree is obtained by calculating the shortest path length between isolated nodes and the nearest cluster. A recursive topological feature vector is constructed based on the cohesion and separation degree. A temporal sliding window extraction is performed on the recursive topological feature vectors, and the temporal interval entropy and spatial distribution entropy of the recursive points within the window are calculated to construct a recursive structure stability index. The decay rate sequence on the time axis of the recursive structure stability index is calculated, the moment when the decay rate exceeds the preset mutation threshold is identified, and the product of the decay amplitude and decay rate at that moment is extracted to obtain the isolation degree.

[0048] In one specific implementation, real-time monitoring data of the drip rate from the infusion device is acquired, and a drip rate deviation metric sequence is calculated. The drip rate deviation metric refers to the percentage difference between the actual drip rate and the target drip rate relative to the target drip rate. For example, when the target drip rate is 60 drops / minute and the actual drip rate is 54 drops / minute, the drip rate deviation metric is -10%. For drip rate monitoring data collected over a period of 12 hours with a sampling interval of 5 seconds, a drip rate deviation metric sequence containing 8640 sampling points is formed.

[0049] Phase space reconstruction is performed on the drop rate deviation metric sequence. The phase space reconstruction employs a time delay method, mapping a one-dimensional time series to a multi-dimensional phase space to reveal the system's dynamic characteristics. Specifically, an embedding dimension of 3 and a time delay of 10 sampling points (i.e., 50 seconds) are chosen. For the i-th point in the sequence, a phase space vector (xi) is constructed. i x i+10 x i+20 Each phase space point is composed of three consecutive drop rate deviation metrics spaced 50 seconds apart from the original sequence. After reconstruction, a three-dimensional phase space trajectory of 8600 points is obtained.

[0050] Calculate the phase distance between each time point in phase space. Use Euclidean distance to calculate the distance between any two phase space points. For example, the phase space distance between the 100th and 200th sampling points is calculated to be 0.15, indicating a high similarity in drop velocity states at these two times. For all 8600 phase space points, calculate the pairwise distances to form an 8600×8600 distance matrix.

[0051] Valid recursive points are determined using a dual screening method based on an adaptive recursion radius threshold and trajectory similarity. The recursion radius threshold employs an adaptive method, with an initial value set to 15% of the average distance between phase space points (0.08 in the example). For each time point i, historical time points whose distance is less than the recursion radius threshold are selected as candidate recursive points. Subsequently, trajectory similarity screening is performed, calculating the dynamic time-warped distance between the trajectory segment formed by the five sampling points before and after each candidate recursive point and the trajectory segment formed by the five sampling points before and after the current time point. Candidate points with a distance less than 0.12 are determined as valid recursive points. For example, the phase space distances between the 3600th sampling point and the 1200th and 2450th sampling points are 0.05 and 0.07 respectively, and the trajectory segment similarity is high; therefore, these two historical points are determined as valid recursive points for the 3600th point.

[0052] Generate a recursive graph matrix. For each time point i, if time point j is a valid recursive point, then the recursive graph matrix R(i, j) = 1; otherwise, R(i, j) = 0. This forms an 8600×8600 binary recursive graph matrix. In the example, the sparsity of the matrix is ​​approximately 92%, indicating that there are no recursive relationships between most time points.

[0053] A topological connectivity graph is constructed for the distribution of recursive points at each time point in the recursion graph matrix. Centered on the current time point, all valid recursive points are considered nodes in the graph. If the time interval between two recursive points is less than 100 sampling points (500 seconds) and the phase space distance is less than 0.1, a connection is established between these two nodes. For example, the recursive point set of the 5000th sampling point includes the 1250th, 1350th, 1750th, 3200th, and 3850th sampling points. A connection is established between the 1250th and 1350th points, and between the 3200th and 3850th points, forming two connected components.

[0054] Clusters and isolated nodes were identified using a graph clustering algorithm. A density-based spatial clustering algorithm was employed, with the core parameters set to a density threshold of 3 and a distance threshold of 0.15. Recursive points were divided into clusters and isolated nodes. For example, in the recursive graph of the 5000th sample point, points 1250, 1350, and 1750 formed one cluster, points 3200 and 3850 formed another cluster, and there were no isolated nodes. However, in the recursive graph of the 7200th sample point, there were 12 recursive points, forming 3 clusters and 2 isolated nodes.

[0055] The cohesion of clusters is calculated by determining the connection density of the clusters. Cohesion is defined as the ratio of the number of connections within a cluster to the number of nodes within the cluster. For example, the first cluster at sampling point 5000 contains 3 nodes and 2 connections, with a cohesion of 2 / 3 = 0.67; the second cluster contains 2 nodes and 1 connection, with a cohesion of 1 / 2 = 0.5. The weighted average of the cohesion of all clusters, with the weights being the proportion of nodes in each cluster, yields an overall cohesion of 0.6 at time point 5000.

[0056] The separation degree is obtained by calculating the shortest path length between an isolated node and its nearest cluster. For each isolated node, the phase space distance from it to the center of each cluster is calculated, and the minimum value is taken as the separation degree. For example, in the recursive graph of the 7200th sampling point, the distances from two isolated nodes to their nearest clusters are 0.22 and 0.25, respectively, and the average value of 0.235 is taken as the separation degree at that time point.

[0057] A recursive topological feature vector is constructed based on cohesion and separation. For each time point i, the recursive topological feature vector F(i) = (cohesion, separation, number of clusters, proportion of isolated nodes). For example, the feature vector of the 5000th sampling point is (0.6, 0, 2, 0), and the feature vector of the 7200th sampling point is (0.75, 0.235, 3, 0.167).

[0058] A temporal sliding window extraction method was used to extract the recursive topological feature vectors. The window size was set to 60 sampling points (5 minutes), and the step size was 10 sampling points (50 seconds). Within each window, the temporal interval entropy and spatial distribution entropy of the recursive points were calculated. The temporal interval entropy reflects the uncertainty of the temporal distribution of the recursive points, while the spatial distribution entropy reflects the uncertainty of the distribution of the recursive points in the phase space. For example, within the window of sampling points 3000 to 3060, the temporal interval entropy was 1.82, and the spatial distribution entropy was 2.13.

[0059] A recursive structure stability index is constructed. The stability index S is defined as the weighted sum of the temporal interval entropy and the spatial distribution entropy, with weights set to 0.4 and 0.6, respectively. For example, the stability index for the sampling window from the 3000th to the 3060th sampling point is 1.82 × 0.4 + 2.13 × 0.6 = 2.006. The stability index is calculated for all windows, forming a stability index sequence of length 857.

[0060] Calculate the decay rate sequence on the time axis for the stability index of the recursive structure. The decay rate is defined as the rate of change of the stability index between adjacent windows. For example, if the stability index of the 30th window is 2.006 and the stability index of the 31st window is 1.876, then the decay rate of the 31st window is (2.006-1.876) / 2.006=6.48%.

[0061] Identify time points where the decay rate exceeds a preset mutation threshold. The mutation threshold is set to 15%, and windows where the decay rate exceeds this threshold are marked as potential outliers. For example, in the 78th window (corresponding to sampling points 780 to 840 of the original sequence), the decay rate is 17.2%, exceeding the mutation threshold and thus identified as a potential outlier.

[0062] The isolation degree is obtained by extracting the product of the decay amplitude and decay rate at abnormal time points. The isolation degree is defined as the product of the decay amplitude (the absolute difference between the stability indices of adjacent windows) and the decay rate. For example, if the stability index of window 78 decreases from 2.324 to 1.924, the decay amplitude is 0.4, the decay rate is 17.2%, and the isolation degree is 0.4 × 17.2% = 0.0688.

[0063] By setting an isolation threshold, the final abnormal time points are selected. The isolation threshold is set to 0.05, and potential abnormal points exceeding this threshold are identified as drip rate abnormal points.

[0064] In one optional implementation, the drop rate deviation metric sequence is reconstructed in phase space. The phase space distance between each time point in the reconstructed phase space and the historical time points is calculated. Valid recursive points are determined based on a dual screening method of adaptive recursion radius threshold and trajectory similarity. The resulting recursion graph matrix includes: The phase space of the drop rate deviation metric sequence is reconstructed. In the reconstructed phase space, a time backtracking window is constructed for each time point. Historical time points are extracted within the time backtracking window, and the Euclidean distance between the time point and the historical time point is calculated to obtain the phase space distance. The sequence of drip rate deviation metrics is segmented and statistically analyzed according to time windows. The standard deviation of the drip rate deviation metrics within each time window is calculated, and the product of the standard deviation and the preset sensitivity coefficient is used as the adaptive recursive radius threshold. Historical time points with phase space distances less than the adaptive recursive radius threshold are selected to construct a candidate recursive point set. For each candidate recursive point, the corresponding subsequent evolution fragment is extracted and a point-to-point distance matrix is ​​formed with the subsequent evolution fragment of the current time point. Based on the distance matrix, dynamic programming is performed to calculate the cumulative distance matrix. The alignment path that minimizes the total distance is extracted from the cumulative distance matrix. The cumulative distance value of the endpoint of the alignment path is extracted to determine the dynamic time warping distance. The dynamic time warping distance is transformed by the inverse to obtain the trajectory similarity. Candidate recursive points whose trajectory similarity exceeds a preset trajectory similarity threshold are marked as valid recursive points. Based on the distribution of valid recursive points on the time axis, a recursion graph matrix is ​​generated.

[0065] In one specific implementation, drip rate data during the infusion process is acquired, and the drip rate value per minute is recorded in real time by a sensor. The infusion drip rate is usually set to a fixed value, such as 60 drops / minute, but the actual drip rate can fluctuate due to various factors. The drip rate deviation metric sequence refers to the degree of difference between the actual measured drip rate and the target drip rate. It is calculated by subtracting the target drip rate from the actual drip rate, dividing by the target drip rate, and then multiplying by 100% to obtain a percentage. For example, if the target drip rate is 60 drops / minute and the actual measured drip rate is 54 drops / minute, then the drip rate deviation metric value is -10%. By continuously collecting drip rate data, a time series of drip rate deviation metrics is formed. In a practical case, the sampling frequency was 10 seconds / time, and continuous monitoring was conducted for 8 hours, resulting in a drip rate deviation metric sequence of 2880 sampling points.

[0066] Phase space reconstruction of the drop rate deviation metric sequence is a fundamental step in identifying the dynamic characteristics of the system. Phase space reconstruction employs a time-delay embedding method to map the one-dimensional time series to a high-dimensional phase space. Specifically, the embedding dimension m=3 and the time delay τ=6 sampling points are chosen. For any time point i in the sequence, a three-dimensional vector is constructed as the phase space point, containing the drop rate deviation metric at time i, time i+6, and time i+12. For example, if the drop rate deviation metric values ​​at times 100, 106, and 112 are 5%, 7%, and 6% respectively, then the phase space point corresponding to time 100 is a three-dimensional point (5%, 7%, 6%). In this way, the original 2880-point one-dimensional time series is mapped to a 2868-point three-dimensional phase space trajectory.

[0067] In the reconstructed phase space, a time backtracking window is constructed for each time point. The length of the time backtracking window is set to 300 sampling points, approximately 50 minutes. For the current time point i in the phase space, the time backtracking window includes all historical time points from time i-300 to time i-1. Historical time points are extracted within the time backtracking window, and the phase space distance between the current time point and historical time points is calculated. The phase space distance is calculated using Euclidean distance, which is the square root of the sum of the squares of the differences in each dimension between two three-dimensional phase space points. For example, if the phase space coordinates of time 1500 are (3%, 2%, 4%) and the phase space coordinates of time 1300 are (2%, 4%, 3%), then the Euclidean distance between the two points is 2.45%. The phase space distance between each historical point within the time backtracking window and the current point is calculated, forming a distance set.

[0068] The drip rate deviation sequence is segmented and statistically analyzed according to time windows, and an adaptive recursive radius threshold is calculated. The time window length is set to 30 sampling points, approximately 5 minutes. The standard deviation of the drip rate deviation within each time window is calculated. The standard deviation reflects the severity of drip rate fluctuations. The adaptive recursive radius threshold is obtained by multiplying the standard deviation by a preset sensitivity coefficient. The sensitivity coefficient is set according to clinical needs; the larger the sensitivity coefficient, the greater the tolerance for drip rate fluctuations. In a specific case, the sensitivity coefficient is set to 1.5. For example, if the standard deviation of the drip rate deviation within the time window [1500, 1529] is 2.3%, then the corresponding adaptive recursive radius threshold for this window is 3.45%. Using an adaptive recursive radius threshold instead of a fixed threshold can better adapt to the drip rate fluctuation characteristics of different patients and at different times, improving the accuracy of early warning.

[0069] Historical time points with a phase space distance less than the adaptive recursion radius threshold are selected to construct a candidate recursion point set. For the current time point i, all historical time points j within the time backtracking window are traversed. If the phase space distance between time point i and j is less than the adaptive recursion radius threshold of the time window containing time point i, then time point j is added to the candidate recursion point set. For example, the adaptive recursion radius threshold for time 1500 is 3.45%, and the phase space distance between time 1300 and time 1500 is 2.45%, which is less than the threshold of 3.45%. Therefore, time 1300 is included in the candidate recursion point set for time 1500.

[0070] For each candidate recursive point, extract the corresponding subsequent evolutionary segment and evaluate the trajectory similarity. The length of the subsequent evolutionary segment is set to 10 sampling points, approximately 1.7 minutes. For the current time point i and the candidate recursive point j, extract 10 consecutive points starting from i and j respectively, forming two evolutionary trajectories. Calculate the dynamic time warping distance between these two trajectories and evaluate the trajectory similarity. The dynamic time warping algorithm first constructs a point-to-point distance matrix. In the distance matrix, the element (a, b) represents the Euclidean distance between the phase space points corresponding to time point i+a and time point j+b. Based on the distance matrix, calculate the cumulative distance matrix through dynamic programming. In the cumulative distance matrix, the element (a, b) represents the minimum cumulative distance from the trajectory starting point (0, 0) to point (a, b). Backtrack from the cumulative distance matrix to extract the alignment path that minimizes the total distance, and extract the cumulative distance value of the endpoint of the alignment path, which is the dynamic time warping distance. Perform a reciprocal transformation on the dynamic time warping distance to obtain the trajectory similarity, specifically 1 divided by (1 plus the dynamic time warping distance). For example, if the dynamic time warped distance between time 1500 and time 1300 is 0.25, then the trajectory similarity is 0.8.

[0071] Candidate recursive points whose trajectory similarity exceeds a preset trajectory similarity threshold are selected and marked as valid recursive points. The trajectory similarity threshold is set to 0.7. For each candidate recursive point j at the current time point i, if the trajectory similarity exceeds 0.7, j is marked as a valid recursive point of i. For example, the trajectory similarity between time 1500 and time 1300 is 0.8, which is greater than the threshold of 0.7, so time 1300 is marked as a valid recursive point of time 1500. By calculating both phase space distance and trajectory similarity for dual selection, it is ensured that the found recursive points are not only close in phase space but also have similar subsequent evolutionary behaviors, thus improving the quality of recursive points.

[0072] A recursion graph matrix is ​​generated based on the distribution of valid recursive points along the time axis. The recursion graph matrix is ​​a binary matrix with a dimension equal to the length of the time series. A value of 1 for the matrix element (i, j) indicates that time j is a valid recursive point at time i; otherwise, it is 0. For example, if time 1300 is a valid recursive point at time 1500, then the element (1500, 1300) in the recursion graph matrix has a value of 1. In practical calculations, for a time series with n sampling points, the recursion graph matrix is ​​an n×n sparse matrix, typically stored in a compressed format, recording only the positions of elements with a value of 1. The recursion graph matrix visually illustrates the recursive relationships between time points in the time series, forming the basis for further analysis of the system's dynamic characteristics.

[0073] In one optional implementation, for the drop rate deviation measurement sequence before the corresponding time of the anomaly identification result, the state transition information transmission amount at different time intervals is calculated to determine the characteristic time scale. Based on the frequency domain features of the characteristic time scale and the isolation degree change rate, a comprehensive anomaly index is calculated to determine the warning level and generate warning information, including: Discretize the sequence of drop rate deviation metrics before the time corresponding to the anomaly identification result, calculate the transition probability between different state values ​​in the sequence of drop rate deviation metrics, calculate the information transmission amount of state transition according to different time intervals, identify the time interval with the largest information transmission amount, and determine the characteristic time scale. The isolation degree change rate sequence is obtained by calculating the difference between adjacent time points of the isolation degree sequence. The frequency domain spectrum is obtained by performing a Fourier transform on the isolation degree change rate sequence. The frequency component with the largest amplitude is identified from the frequency domain spectrum as the dominant frequency. The synchronization deviation degree is obtained by calculating the deviation ratio between the period corresponding to the dominant frequency and the characteristic time scale. The energy values ​​of the frequency components other than the dominant frequency in the frequency domain spectrum are accumulated. The ratio of the accumulated energy value to the energy value of the dominant frequency is used to obtain the spectral dispersion. The synchronization deviation and the spectral dispersion are weighted and summed to obtain the comprehensive anomaly index. The warning level is determined according to the numerical range of the comprehensive anomaly index. The warning level is encapsulated with the corresponding time of the anomaly identification result to generate warning information.

[0074] In one specific implementation, a sequence of drop rate deviation metrics is obtained up to the time corresponding to the anomaly identification result. The drop rate deviation metric represents the degree of difference between the actual drop rate and the set drop rate, and is usually calculated as (actual drop rate - set drop rate) / set drop rate. For example, if the set drop rate is 60 drops / minute, and the actual measured drop rate is 66 drops / minute, then the drop rate deviation metric is 0.1, indicating that the actual drop rate is 10% higher than the set value.

[0075] The acquired drop rate deviation metric sequence is discretized, dividing the continuous deviation metric values ​​into a finite number of states. Discretization can employ equal-width or equal-frequency binning methods, or it can be based on a threshold set according to domain knowledge. In this embodiment, the drop rate deviation metric values ​​are divided into 5 states: severely low (below -0.2), slightly low (-0.2 to -0.05), normal (-0.05 to 0.05), slightly high (0.05 to 0.2), and severely high (above 0.2).

[0076] The probability of transitions between different state values ​​in the sequence of drop rate deviations is calculated. For a state transition with a time interval of τ, the conditional probability P(j|i, τ) of transitioning from state i to state j is calculated. For example, for the case of τ=1, the state transitions at adjacent time points in the sequence are statistically analyzed; for the case of τ=2, the state transitions at a time interval of one time point are statistically analyzed, and so on.

[0077] Based on the state transition probability, the information transfer amount is calculated for different time intervals. The information transfer amount is used to quantify the degree of uncertainty reduction in state transitions at a specific time interval. For each time interval τ, its information transfer amount T(τ) is calculated. As the time interval τ gradually increases, the information transfer amount usually increases first and then decreases, reaching a maximum value at a certain specific time interval.

[0078] By comparing the amount of information transmitted at different time intervals, the time interval τmax with the largest amount of information transmitted is identified and determined as the characteristic time scale. The characteristic time scale reflects the inherent periodicity of system state changes and is of great significance for anomaly early warning. For example, if the amount of information transmitted is largest when τ=3, then the characteristic time scale is 3 time units.

[0079] The isolation degree sequence is analyzed. Isolation degree represents the degree of dissimilarity between a data point and its surrounding data points; outliers typically have high isolation degrees. The difference in isolation degree between adjacent time points is calculated to obtain the isolation degree change rate sequence ΔI(t). The isolation degree change rate sequence reflects the dynamic trend of the degree of anomaly.

[0080] A Fourier transform is performed on the isolation degree change rate sequence to convert it to the frequency domain for analysis. The Fourier transform decomposes the time series into periodic components of different frequencies. The frequency component with the largest amplitude is identified from the frequency domain spectrum; this frequency is labeled as the dominant frequency (fdom). The period corresponding to the dominant frequency is T = 1 / fdom, representing the main periodic pattern of isolation degree change.

[0081] The synchronization deviation degree Dsync is obtained by calculating the ratio of the deviation between the dominant frequency corresponding to the period T and the characteristic time scale τmax. The formula for calculating the synchronization deviation degree is Dsync = |T - τmax| / max(T, τmax). The smaller the synchronization deviation degree, the closer the periodicity of the isolation degree change is to the inherent periodicity of the system state transition, and the lower the probability of anomalies.

[0082] The energy values ​​of all frequency components in the frequency domain spectrum except the dominant frequency are accumulated to obtain Enon-dom. The ratio of the accumulated energy value to the energy value Edom of the dominant frequency is used to obtain the spectral dispersion Dispec = Enon-dom / Edom. The larger the spectral dispersion, the more unstable the periodicity of the isolation degree change, and the higher the probability of anomalies.

[0083] The comprehensive anomaly index (CAI) is obtained by weighted summation of synchronization deviation and spectral dispersion. The calculation formula is CAI = w1 × Dsync + w2 × Dspec, where w1 and w2 are weighting coefficients, and w1 + w2 = 1. In this embodiment, w1 = 0.4 and w2 = 0.6 are set.

[0084] The warning level is determined based on the numerical range of the Comprehensive Anomaly Index (CAI). Warning levels are typically divided into four levels: Normal (CAI < 0.3), Mild Anomaly (0.3 ≤ CAI < 0.5), Moderate Anomaly (0.5 ≤ CAI < 0.7), and Severe Anomaly (CAI ≥ 0.7). Different levels correspond to different handling strategies; for example, severe anomalies require immediate manual intervention.

[0085] The warning level is encapsulated with the corresponding time of the anomaly identification result to generate warning information. The warning information includes key information such as warning time, warning level, anomaly type, drip rate deviation value, and characteristic time scale, which facilitates medical staff to respond quickly and take appropriate measures.

[0086] In this embodiment, combining characteristic time scales and frequency domain analysis enables a more accurate assessment of the severity of anomalies, reduces false alarms, and improves the reliability of the early warning system. Especially in infusion systems, timely and accurate early warnings can significantly reduce medical risks.

[0087] like Figure 2 The diagram shown illustrates the logic flowchart for anomaly warning.

[0088] In one optional implementation, the discretization process includes: A state space is constructed for the drop rate deviation metric sequence. Attractor centers with numerical clusters are identified in the state space. An attraction domain boundary is formed by expanding outward from each attractor center. The attraction strength from any point in the state space to each attractor center is calculated. The state label corresponding to the attractor with the largest attraction strength is assigned to the corresponding point. Based on the attractor affiliation relationship, the drop rate deviation metric sequence is mapped to a discrete state sequence.

[0089] In one specific implementation, drip rate data is collected using an infusion monitoring device, typically sampling every 10 seconds. The drip rate deviation metric refers to the percentage difference between the actual drip rate and the target drip rate. For example, if the target drip rate is set to 60 drops / minute and the actual measured value is 54 drops / minute, the drip rate deviation metric is -10%. In clinical applications, continuous monitoring of an 8-hour infusion process yields 2880 sampling points, forming a sequence of drip rate deviation metrics.

[0090] A state space is constructed from the sequence of drip rate deviation metrics. The state space is a multi-dimensional space describing the dynamic characteristics of drip rate changes, constructed from a one-dimensional time series using a time-delay embedding method. The embedding dimension is chosen to be 3, and the time delay is 6 sampling points (i.e., 1 minute). For the i-th point in the sequence, a three-dimensional state space vector is constructed, containing the drip rate deviation metrics at points i-th, i+6-th, and i+12-th. For example, if the drip rate deviation metrics at points 100, 106, and 112 in the sequence are -5%, -3%, and -4% respectively, then the state space vector corresponding to point 100 is (-5%, -3%, -4%). ​​In this way, the original 2880-point one-dimensional sequence is mapped to a 2868-point three-dimensional state space point set. The state space visually displays the dynamic trajectory of drip rate changes, with different drip rate states represented as point sets in different regions of the state space.

[0091] In the state space, attractor centers representing numerical clusters are identified. Attractors are regions with high point density in the state space, representing steady-state drop rates. Density clustering is used to identify attractor centers. Algorithm parameters include a neighborhood radius ε of 5% and a minimum number of points MinPts of 30. For each point in the state space, the number of points in its ε-neighborhood is calculated, and points with a number greater than MinPts are selected as core points. Directly density-reachable core points are connected to form clusters, and the centroid of each cluster serves as an attractor center. In the practical case, four attractor centers were identified in the state space, with coordinates A1(2%, 3%, 2%), A2(-3%, -4%, -3%), A3(8%, 7%, 9%), and A4(-10%, -12%, -9%), representing four drop rate states: slightly positively skewed, slightly negatively skewed, significantly positively skewed, and significantly negatively skewed, respectively.

[0092] An attraction domain boundary is formed by expanding outwards from each attractor center. The attraction domain is the range of influence of the attractor; the attraction force at the boundary decays to a specific threshold. Boundary determination employs an adaptive method based on the distribution density gradient of points surrounding the attractor center. For each attractor center, surrounding points are arranged from nearest to farthest, and the density change rate between adjacent points is calculated. When the density change rate first falls below a threshold of 0.05, the corresponding distance is determined as the attraction domain radius. For example, the attraction domain radius of attractor center A1 is 7%, A2's is 8%, A3's is 6%, and A4's is 9%. The difference in attraction domain radii among different attractors reflects the difference in drop rate stability; a larger radius indicates a larger range of state fluctuations.

[0093] Calculate the attraction strength from any point in the state space to the centers of each attractor. The attraction strength depends on the distance from the point to the attractor center and the attractor strength. The distance is calculated using Euclidean distance, and the attractor strength is determined by the number of points contained within the attractor. The formula for calculating the attraction strength is the attractor strength divided by the square of the distance from the point to the attractor center. For example, if the distances from point P(5%, 4%, 6%) in the state space to the centers of the four attractors are 4.36%, 11.53%, 4.58%, and 19.67%, respectively, and the corresponding attractor strengths are 850, 920, 480, and 350, the calculated attraction strengths for point P are 44.67, 6.92, 22.91, and 0.90.

[0094] Assign the state label corresponding to the attractor with the strongest attraction to the corresponding point. Each attractor represents a drop velocity state, labeled S1, S2, S3, and S4 respectively. For each point in the state space, calculate the attraction strength from each attractor, and assign the state label of the attractor with the strongest attraction strength to that point. For example, the strongest attraction at point P comes from attractor A1, so the state label of point P is S1. When a point is located in the boundary region of multiple attraction domains and the attraction strengths are similar, it indicates that the drop velocity is in a state transition process and requires special attention.

[0095] The drop rate deviation metric sequence is mapped to a discrete state sequence based on attractor attribution relationships. For each corresponding state space point in the original drop rate deviation metric sequence, its state label is determined according to the principle of maximum attraction, generating a discrete state sequence of equal length. For example, the drop rate deviation metric values ​​at points 500 to 510 in the original sequence correspond to state labels S2, S2, S2, S2, S2, S1, S1, S1, S1, S1, S1, indicating that the drop rate state changes from S2 to S1. The discrete state sequence simplifies the expression of drop rate changes and facilitates subsequent analysis of state transition patterns.

[0096] State transition analysis is a crucial step in anomaly early warning. By statistically analyzing the frequency and timing of transitions between different states, a normal state transition model is established. Under normal circumstances, the dripping state slowly alternates between S1 and S2, with each state lasting a relatively long time, typically 30-60 minutes. In anomalies, however, states switch frequently or rapidly jump to S3 or S4, with short durations, typically less than 5 minutes. For example, if a state sequence is detected to exhibit more than three rapid jumps from S1 to S3 or S2 to S4 within 10 minutes, it is identified as an anomaly, triggering an early warning.

[0097] In this embodiment, the continuous sequence of drip rate deviation values ​​is converted into a discrete state sequence, which simplifies the data expression while retaining key dynamic characteristics. By adaptively identifying the drip rate state characteristics of different patients through an attractor model, the trend of infusion drip rate changes is finely characterized, enabling early and accurate warning of drip rate abnormalities and improving infusion safety.

[0098] The adaptive algorithm system for early warning of abnormal infusion drip rate in this embodiment of the invention includes: The data acquisition unit is used to acquire the drip rate monitoring data stream during the infusion process. The drip rate monitoring data stream includes time-series drip rate values ​​and corresponding time identifiers. The baseline reconstruction unit is used to construct the phase space trajectory matrix of the time-series drop rate values ​​within the historical time window, and to perform singular value decomposition to obtain the singular value spectrum. The mode boundary point is determined by the bidirectional decay gradient field, and the singular value spectrum is decomposed into the dominant mode and the perturbation mode. The drop rate center baseline is reconstructed and the drop rate fluctuation boundary is determined. The deviation calculation unit is used to calculate the difference between the real-time drip rate value and the drip rate center baseline, and to perform a ratio operation with the drip rate fluctuation boundary to obtain the drip rate deviation metric value. The drip rate deviation metric values ​​within a continuous time window are organized into a drip rate deviation sequence. An anomaly identification unit is used to reconstruct the phase space of the drop rate deviation measurement value sequence and generate a recursive graph matrix. Based on the recursive graph matrix, a topological connectivity graph is constructed and recursive topological feature vectors are extracted. The isolation degree is determined by the decay rate of the recursive structure stability index. When the isolation degree exceeds the preset isolation threshold, the anomaly identification result is output. The early warning generation unit is used to calculate the state transition information transmission amount at different time intervals from the drop rate deviation measurement value sequence before the corresponding time of the anomaly identification result, determine the characteristic time scale, calculate the comprehensive anomaly index based on the frequency domain features of the characteristic time scale and the isolation degree change rate, determine the early warning level, and generate early warning information.

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

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

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

[0102] 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. An adaptive algorithm method for early warning of abnormal infusion drip rate, characterized in that, include: Acquire a drip rate monitoring data stream during the infusion process, wherein the drip rate monitoring data stream includes time-series drip rate values ​​and corresponding time identifiers; A phase space trajectory matrix is ​​constructed for the time series drop rate values ​​within the historical time window, and singular value decomposition is performed to obtain the singular value spectrum. The mode boundary point is determined by the bidirectional decay gradient field, and the singular value spectrum is decomposed into the dominant mode and the perturbation mode. The drop rate center baseline is reconstructed and the drop rate fluctuation boundary is determined. The difference between the real-time drip rate value and the drip rate center baseline is calculated and compared with the drip rate fluctuation boundary to obtain the drip rate deviation metric value. The drip rate deviation metric values ​​within a continuous time window are organized into a drip rate deviation sequence. The phase space of the drop rate deviation metric sequence is reconstructed and a recursive graph matrix is ​​generated. Based on the recursive graph matrix, a topological connectivity graph is constructed and recursive topological feature vectors are extracted. The isolation degree is determined by the decay rate of the recursive structure stability index. When the isolation degree exceeds the preset isolation threshold, the anomaly identification result is output. For the sequence of drop rate deviation measurements before the corresponding time point of the anomaly identification result, calculate the state transition information transmission amount at different time intervals, determine the characteristic time scale, calculate the comprehensive anomaly index based on the frequency domain features of the characteristic time scale and the isolation degree change rate, determine the warning level and generate warning information.

2. The method according to claim 1, characterized in that, A phase space trajectory matrix is ​​constructed for the time-series drop rate values ​​within the historical time window, and singular value decomposition is performed to obtain the singular value spectrum. The mode boundary point is determined by a two-way decaying gradient field. The singular value spectrum is decomposed into the dominant mode and the perturbation mode. The drop rate center baseline is reconstructed and the drop rate fluctuation boundary is determined, including: The mutual information function under different time delays is calculated for the time-series drop rate values ​​within the historical time window. The time delay at which the mutual information function first reaches its minimum value is identified, and the embedding delay time is determined. The pseudo-nearest neighbor ratio is calculated for the time-series drop rate values ​​under different embedding dimensions. The embedding dimension at which the pseudo-nearest neighbor ratio drops below a preset stable threshold is identified, and the phase space embedding dimension is determined. Based on the embedding delay time and the phase space embedding dimension, the time-series drop rate values ​​are expanded according to the time sliding window to generate the phase space trajectory matrix. Singular value decomposition is performed on the phase space trajectory matrix to obtain the singular value spectrum arranged in descending order of energy and the corresponding left and right singular vector sets; A bidirectional decaying gradient field is constructed on the singular value spectrum to obtain a decaying acceleration sequence. Based on the decaying acceleration sequence, a candidate mode boundary set is extracted, and the mode separation index is calculated to determine the mode boundary point. Singular values ​​and their corresponding left and right singular vectors before the mode boundary are classified as dominant modes, while singular values ​​and their corresponding left and right singular vectors after the mode boundary are classified as perturbation modes. The drop rate center baseline is reconstructed by performing matrix multiplication on the singular values ​​and left and right singular vectors in the dominant mode and projecting them back to the original time dimension. The mean and standard deviation of the singular values ​​in the perturbation mode are calculated, and confidence intervals are constructed to determine the drop rate fluctuation boundary.

3. The method according to claim 2, characterized in that, A bidirectional decaying gradient field is constructed from the singular value spectrum to obtain a decaying acceleration sequence. Based on the decaying acceleration sequence, a set of candidate mode boundaries is extracted, and the mode separation index is calculated to determine the mode boundary points, including: A two-way decay gradient field is constructed for the singular value spectrum. The forward energy decay gradient and the backward energy decay gradient are calculated at each singular value position. The difference between the forward energy decay gradient and the backward energy decay gradient is used to obtain the gradient symmetry deviation sequence. The decay acceleration sequence is obtained by stepwise recursive difference. An adaptive detection threshold is constructed based on the statistical distribution of the decay acceleration sequence. The positions of singular values ​​in the decay acceleration sequence that exceed the adaptive detection threshold are extracted to form a set of candidate mode boundaries. For each candidate modal boundary in the candidate modal boundary set, an energy segmentation topology is constructed. The energy distribution entropy difference between the front and back subspaces of the candidate modal boundary is calculated. The energy distribution entropy difference is coupled with the energy proportion variance to obtain the modal separation index. The candidate modal boundary with the largest modal separation index is selected to determine the modal boundary point.

4. The method according to claim 1, characterized in that, The drop rate deviation metric sequence is reconstructed in phase space to generate a recursive graph matrix. A topological connectivity graph is constructed based on the recursive graph matrix, and recursive topological feature vectors are extracted. The isolation degree is determined by the decay rate of the recursive structure stability index, including: The phase space of the drip rate deviation measurement sequence is reconstructed, the phase space distance between each time point in the reconstructed phase space and the historical time point is calculated, and the effective recursion points are determined based on the dual screening of adaptive recursion radius threshold and trajectory similarity, and a recursion graph matrix is ​​generated. A topological connectivity graph is constructed based on the distribution of recursive points at each time point in the recursive graph matrix. Clusters and isolated nodes are identified by graph clustering algorithm. The cohesion of clusters is obtained by calculating the connection density of clusters, and the separation degree is obtained by calculating the shortest path length between isolated nodes and the nearest cluster. A recursive topological feature vector is constructed based on the cohesion and separation degree. A temporal sliding window extraction is performed on the recursive topological feature vectors, and the temporal interval entropy and spatial distribution entropy of the recursive points within the window are calculated to construct a recursive structure stability index. The decay rate sequence on the time axis of the recursive structure stability index is calculated, the moment when the decay rate exceeds the preset mutation threshold is identified, and the product of the decay amplitude and decay rate at that moment is extracted to obtain the isolation degree.

5. The method according to claim 4, characterized in that, The drop rate deviation metric sequence is reconstructed in phase space. The phase space distance between each time point in the reconstructed phase space and the historical time points is calculated. Valid recursive points are determined based on a dual screening method of adaptive recursion radius threshold and trajectory similarity. The resulting recursion graph matrix includes: The phase space of the drop rate deviation metric sequence is reconstructed. In the reconstructed phase space, a time backtracking window is constructed for each time point. Historical time points are extracted within the time backtracking window, and the Euclidean distance between the time point and the historical time point is calculated to obtain the phase space distance. The sequence of drip rate deviation metrics is segmented and statistically analyzed according to time windows. The standard deviation of the drip rate deviation metrics within each time window is calculated, and the product of the standard deviation and the preset sensitivity coefficient is used as the adaptive recursive radius threshold. Historical time points with phase space distances less than the adaptive recursive radius threshold are selected to construct a candidate recursive point set. For each candidate recursive point, the corresponding subsequent evolution fragment is extracted and a point-to-point distance matrix is ​​formed with the subsequent evolution fragment of the current time point. Based on the distance matrix, dynamic programming is performed to calculate the cumulative distance matrix. The alignment path that minimizes the total distance is extracted from the cumulative distance matrix. The cumulative distance value of the endpoint of the alignment path is extracted to determine the dynamic time warping distance. The dynamic time warping distance is transformed by the inverse to obtain the trajectory similarity. Candidate recursive points whose trajectory similarity exceeds a preset trajectory similarity threshold are marked as valid recursive points. Based on the distribution of valid recursive points on the time axis, a recursion graph matrix is ​​generated.

6. The method according to claim 1, characterized in that, For the drop rate deviation measurement sequence before the corresponding time point of the anomaly identification result, calculate the state transition information transmission amount at different time intervals, determine the characteristic time scale, calculate the comprehensive anomaly index based on the frequency domain features of the characteristic time scale and the isolation degree change rate, determine the warning level, and generate warning information, including: Discretize the sequence of drop rate deviation metrics before the time corresponding to the anomaly identification result, calculate the transition probability between different state values ​​in the sequence of drop rate deviation metrics, calculate the information transmission amount of state transition according to different time intervals, identify the time interval with the largest information transmission amount, and determine the characteristic time scale. The isolation degree change rate sequence is obtained by calculating the difference between adjacent time points of the isolation degree sequence. The frequency domain spectrum is obtained by performing a Fourier transform on the isolation degree change rate sequence. The frequency component with the largest amplitude is identified from the frequency domain spectrum as the dominant frequency. The synchronization deviation degree is obtained by calculating the deviation ratio between the period corresponding to the dominant frequency and the characteristic time scale. The energy values ​​of the frequency components other than the dominant frequency in the frequency domain spectrum are accumulated. The ratio of the accumulated energy value to the energy value of the dominant frequency is used to obtain the spectral dispersion. The synchronization deviation and the spectral dispersion are weighted and summed to obtain the comprehensive anomaly index. The warning level is determined according to the numerical range of the comprehensive anomaly index. The warning level is encapsulated with the corresponding time of the anomaly identification result to generate warning information.

7. The method according to claim 6, characterized in that, The discretization process includes: A state space is constructed for the drop rate deviation metric sequence. Attractor centers with numerical clusters are identified in the state space. An attraction domain boundary is formed by expanding outward from each attractor center. The attraction strength from any point in the state space to each attractor center is calculated. The state label corresponding to the attractor with the largest attraction strength is assigned to the corresponding point. Based on the attractor affiliation relationship, the drop rate deviation metric sequence is mapped to a discrete state sequence.

8. An adaptive algorithm system for early warning of abnormal infusion drip rate, used to implement the method of any one of claims 1-7, characterized in that, include: The data acquisition unit is used to acquire the drip rate monitoring data stream during the infusion process. The drip rate monitoring data stream includes time-series drip rate values ​​and corresponding time identifiers. The baseline reconstruction unit is used to construct the phase space trajectory matrix of the time-series drop rate values ​​within the historical time window, and to perform singular value decomposition to obtain the singular value spectrum. The mode boundary point is determined by the bidirectional decay gradient field, and the singular value spectrum is decomposed into the dominant mode and the perturbation mode. The drop rate center baseline is reconstructed and the drop rate fluctuation boundary is determined. The deviation calculation unit is used to calculate the difference between the real-time drip rate value and the drip rate center baseline, and to perform a ratio operation with the drip rate fluctuation boundary to obtain the drip rate deviation metric value. The drip rate deviation metric values ​​within a continuous time window are organized into a drip rate deviation sequence. An anomaly identification unit is used to reconstruct the phase space of the drop rate deviation measurement value sequence and generate a recursive graph matrix. Based on the recursive graph matrix, a topological connectivity graph is constructed and recursive topological feature vectors are extracted. The isolation degree is determined by the decay rate of the recursive structure stability index. When the isolation degree exceeds the preset isolation threshold, the anomaly identification result is output. The early warning generation unit is used to calculate the state transition information transmission amount at different time intervals from the drop rate deviation measurement value sequence before the corresponding time of the anomaly identification result, determine the characteristic time scale, calculate the comprehensive anomaly index based on the frequency domain features of the characteristic time scale and the isolation degree change rate, determine the early warning level, and generate early warning information.

9. 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 7.

10. 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 7.

Images