Endocrinology multi-parameter physiological signal fusion analysis system and method
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- THE SECOND HOSPITAL OF HEBEI MEDICAL UNIV
- Filing Date
- 2026-03-06
- Publication Date
- 2026-06-26
Smart Images
Figure CN122271968A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of medical data processing technology, specifically to a multi-parameter physiological signal fusion analysis system and method for endocrinology. Background Technology
[0002] The homeostatic regulation of the human endocrine and metabolic system is a complex physiological process that relies on real-time feedback control by the autonomic nervous system. In clinical practice, the diagnosis and assessment of diabetes and related metabolic syndromes typically depend on the independent analysis of single-dimensional physiological signals, such as recording blood glucose fluctuations through continuous glucose monitoring systems or assessing autonomic nervous function by recording heart rate variability through Holter monitoring. However, this independent analysis model severs the dynamic coupling between neural regulation as a control input and blood glucose concentration as a system output, making it difficult to reveal the deep pathological mechanisms of disease from a systems theory perspective.
[0003] Although existing research attempts to jointly analyze heart rate variability and blood glucose data, it mostly employs traditional methods based on linear correlation or frequency domain analysis. Because human physiological signals inherently possess significant nonlinear and non-stationary characteristics, conventional Fourier transform or linear regression models cannot accurately capture the dynamic time-varying properties of the neuro-metabolic coupling process. Furthermore, significant technical challenges exist in the fusion and processing of multi-source heterogeneous data.
[0004] On the one hand, interstitial fluid glucose sensors and electrocardiogram (ECG) acquisition devices differ fundamentally in sampling frequency, signal transmission path, and biophysical response mechanism, resulting in a non-fixed physical lag time between the two sets of signals. Directly performing alignment analysis based on timestamps would introduce severe phase errors, leading to inverted causal relationships or distorted coupling characteristics.
[0005] On the other hand, when extracting specific metabolic regulatory components from heart rate variability signals, commonly used empirical mode decomposition algorithms are prone to mode aliasing, meaning that signals of the same physical process are dispersed into different mode components, or the same mode component contains features of different time scales. This makes it difficult to accurately separate the neural regulatory components that are truly related to blood glucose fluctuations. At the same time, existing phase space reconstruction methods lack effective signal regularization and closure mechanisms. The generated trajectories usually contain a lot of high-frequency noise and are in an open-loop state, making it impossible to calculate quantitative indicators with clear physical meaning (such as the area of hysteresis loops) using mathematical tools such as Green's theorem. Clinicians can only rely on subjective experience to qualitatively describe the shape of the graph, lacking objective and unified quantitative evaluation standards. Summary of the Invention
[0006] To address the shortcomings of existing technologies, this invention provides a multi-parameter physiological signal fusion analysis system and method for endocrinology, which solves the problems in existing neural-metabolic coupling analysis techniques, such as feature extraction distortion caused by physical lag of multi-source heterogeneous signals and nonlinear mode aliasing, and the lack of objective quantitative evaluation indicators with clear physical meaning.
[0007] To achieve the above objectives, the first aspect of the present invention provides a multi-parameter physiological signal fusion analysis system for endocrinology, which mainly consists of a data acquisition terminal, a processing server, and an output terminal.
[0008] The data acquisition terminal is configured to acquire multidimensional physiological parameters of the subject in real time, specifically including raw blood glucose sequences for monitoring interstitial fluid glucose concentration, and heart rate variability sequences for reflecting the activity of the autonomic nervous system.
[0009] The processing server communicates with the data acquisition terminal and integrates functional modules such as preprocessing calibration, mode decomposition, optimization decision, reconstruction and feature quantization.
[0010] The output terminal is used to visually display the processed phase space trajectory and corresponding quantitative evaluation indicators.
[0011] The preprocessing calibration module is configured to calculate the first-order differential cross-correlation function between the original blood glucose sequence and the heart rate variability sequence, thereby determining the physical lag time between the two sets of signals, and performing time displacement compensation and normalization processing on the original blood glucose sequence to achieve time alignment of multi-source signals.
[0012] The mode decomposition module is configured to perform ensemble empirical mode decomposition on heart rate variability sequences, generating an intrinsic mode function set that suppresses mode aliasing by repeatedly adding white noise sequences to assist the decomposition process.
[0013] The optimization decision module is configured to screen effective modal components and combine them with blood glucose sequences to construct temporary phase space trajectories. The module calculates the topological closure index, which characterizes the closeness and smoothness of the trajectory, selects the component corresponding to the maximum value of the index as the optimal modal component, and determines whether the regulation state is a steady-state closed state or an open-loop divergent state based on the comparison result of the index and the confidence threshold.
[0014] The reconstruction module is configured to construct a target phase space based on the optimal modal components and the blood glucose sequence, and introduce a virtual damping model based on the weight of the blood glucose change rate to normalize the optimal modal components. The smoothing coefficient of the virtual damping model is positively correlated with the absolute change rate of the blood glucose sequence, that is, the damping effect is reduced when the absolute change rate increases and the damping effect is enhanced when the absolute change rate decreases.
[0015] In addition, when the state is determined to be a steady-state closed state, the module performs trajectory forced closure processing by eliminating the first and last residuals.
[0016] The feature quantization module is configured to calculate quantization indicators in response to the adjustment state: Under steady-state closed-loop conditions, the area of the hysteresis loop in the phase space trajectory is calculated using the discrete Green's formula to characterize the energy loss during the regulation process. In the open-loop divergence state, the ratio of the Euclidean distance between the starting and ending points of the trajectory to the trajectory scale is calculated to obtain the phase space divergence index.
[0017] A second aspect of this invention provides a method for multi-parameter physiological signal fusion analysis in endocrinology. This method is based on the aforementioned system and includes the following steps: The raw blood glucose and heart rate variability sequences of the subjects were obtained, the first-order difference cross-correlation function of the two sets of sequences was calculated to determine the physical lag time, and time series correction and normalization were performed accordingly. The ensemble empirical mode decomposition algorithm is used to decompose the corrected heart rate variability sequence to generate an eigenmode function set; Effective modal components are screened and temporary phase space trajectories are constructed by combining blood glucose sequences. The topological closure index is calculated, and the optimal modal components are identified and the regulatory state is determined accordingly. If the state is determined to be a steady-state closed state, a virtual damping model based on the weight of blood glucose change rate is introduced to normalize the optimal mode components, and a forced closure process is performed to generate a closed hysteresis loop, and the area of the hysteresis loop is calculated. If the open-loop divergence state is determined, the forced closure process is stopped, and the normalized drift distance of the first and last points of the phase space trajectory is calculated to obtain the phase space divergence index.
[0018] This invention provides a multi-parameter physiological signal fusion analysis system and method for endocrinology. It has the following beneficial effects: 1. This invention determines and compensates for the physical lag time between the original blood glucose sequence and the heart rate variability sequence by calculating the first-order difference cross-correlation function between the two. This setting solves the problem of time synchronization issues caused by differences in the response of acquisition devices and different physiological transmission path lengths of multi-source physiological signals, ensuring that the data of the two dimensions maintain the consistency of the time reference when constructing the phase space, and effectively avoiding phase space trajectory distortion and coupling feature analysis errors caused by time sequence misalignment.
[0019] 2. This invention employs an ensemble empirical mode decomposition algorithm combined with a topological closure index optimization strategy. The ensemble empirical mode decomposition, by adding white noise to assist analysis, suppresses the mode aliasing phenomenon commonly found in the decomposition of non-stationary signals. Meanwhile, the topological closure index provides an objective evaluation standard based on geometric morphological closure, which can automatically select the mode component with the highest coupling degree to the blood glucose fluctuation cycle from the complex set of intrinsic mode functions, overcoming the uncertainty of traditional methods that rely on manual experience to select principal components.
[0020] 3. This invention is based on a virtual damping model with weighted blood glucose change rate and a hysteresis loop area quantification method. The virtual damping model dynamically adjusts the smoothing coefficient according to the rate of signal change, effectively filtering out high-frequency measurement noise while retaining the details of signal mutations that reflect pathological characteristics. The calculation of the hysteresis loop area transforms the abstract neuro-metabolic regulatory lag effect into a specific numerical indicator, which can quantify the degree of energy loss in the regulatory process, providing a comparable objective basis for clinical assessment of autonomic nerve function and insulin sensitivity. Attached Figure Description
[0021] Figure 1 This is a system architecture diagram of the present invention; Figure 2 This is a flowchart of the method of the present invention; Figure 3 A waveform diagram comparing the modal decomposition effect and the traditional EMD decomposition effect provided in the embodiments of the present invention; Figure 4 A schematic diagram comparing the phase space reconstruction trajectory based on dynamic damping with the traditional linear filtering trajectory provided in an embodiment of the present invention; Figure 5 Box plot of statistical distribution of hysteresis loop energy loss index in different populations provided in embodiments of the present invention.
[0022] Among them, 100 is the data acquisition terminal; 101 is the metabolic acquisition unit; 102 is the neural acquisition unit; 200 is the processing server; 201 is the preprocessing calibration module; 202 is the modality decomposition module; 203 is the optimization decision module; 204 is the reconstruction module; 205 is the feature quantization module; and 300 is the output terminal. Detailed Implementation
[0023] 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.
[0024] Please see the appendix Figure 1This invention provides a multi-parameter physiological signal fusion analysis system for endocrinology, which includes a data acquisition terminal 100, a processing server 200, and an output terminal 300.
[0025] The data acquisition terminal 100 is used to acquire the multidimensional physiological parameter sequence of the subject in real time.
[0026] The processing server 200 communicates with the data acquisition terminal 100 to execute data cleaning, mode decomposition, topology calculation and feature reconstruction algorithms.
[0027] The output terminal 300 is used to display the reconstructed phase space trajectory and quantitative evaluation indicators.
[0028] The data acquisition terminal 100 includes a metabolic acquisition unit 101 and a neural acquisition unit 102.
[0029] The metabolic acquisition unit 101 is configured to continuously monitor the glucose concentration in the interstitial fluid of the subject and generate a raw blood glucose sequence at a first sampling frequency.
[0030] The neural acquisition unit 102 is configured to acquire the subject's electrocardiogram signal or photoplethysmography pulse wave signal and extract the heart rate variability sequence at the second sampling frequency.
[0031] The processing server 200 includes a memory and a processor. The memory stores a computer program, and the processor executes the computer program to implement the functions of the following logical modules: The preprocessing calibration module 201, the mode decomposition module 202, the optimization decision module 203, the reconstruction module 204, and the feature quantization module 205 are included.
[0032] The preprocessing calibration module 201 is used to receive the original blood glucose sequence and heart rate variability sequence. This module resamples and normalizes the heterogeneous data, calculates the instrument transmission lag time between the two sets of signals based on the derivative cross-correlation algorithm, applies time shift compensation to the original blood glucose sequence, and generates the calibrated blood glucose sequence and heart rate variability sequence.
[0033] The mode decomposition module 202 is connected to the preprocessing calibration module 201. This module uses the ensemble empirical mode decomposition algorithm to perform time-frequency domain decomposition on the calibrated heart rate variability sequence, generating a mode set containing multiple intrinsic mode functions and residual terms.
[0034] The optimization decision module 203 is connected to the mode decomposition module 202. This module uses Hilbert transform to calculate the instantaneous frequency of each intrinsic mode function, generates an effective mode subset based on a preset effective physiological frequency band, and combines each component in the effective mode subset with the calibrated blood glucose sequence to construct a temporary phase space trajectory. It then calculates the topological closure index of each temporary phase space trajectory. This module selects the mode component corresponding to the maximum value of the topological closure index as the optimal mode component, compares the maximum value with a preset confidence threshold, and determines whether the endocrine regulation process belongs to a steady-state closed state or an open-loop divergent state.
[0035] The reconstruction module 204 is connected to the optimization decision module 203. This module constructs the target phase space based on the optimal modal components and the calibrated blood glucose sequence. During the construction process, the optimal modal components are numerically normalized using a virtual damping model based on the blood glucose change rate weight to generate the final phase space trajectory coordinates.
[0036] Feature quantization module 205 is connected to reconstruction module 204. This module executes calculation logic based on the state determination result output by optimization decision module 203. When the steady-state closed state is determined, the area of the hysteresis loop enclosed by the phase space trajectory is calculated using the discrete Green's formula. When the open-loop divergence state is determined, the ratio of the Euclidean distance between the first and last points of the phase space trajectory to the total length of the trajectory is calculated to obtain the phase space divergence index.
[0037] See attached document Figure 2 This invention provides a method for multi-parameter physiological signal fusion analysis in endocrinology, comprising the following steps: S100, the preprocessing calibration module 201 acquires the raw blood glucose sequence and heart rate variability sequence, and performs interpolation, alignment and normalization operations; Calculate the first-order difference cross-correlation function of the two sets of sequences, determine the physical lag time of the instrument, and perform time-series correction on the blood glucose sequence; S200 and mode decomposition module 202 perform ensemble empirical mode decomposition on the preprocessed heart rate variability sequence and output a set of full-band intrinsic mode functions. S300, the optimization decision module 203 removes modal components whose main frequency falls outside the preset physiological frequency band; For the remaining effective modal components, temporary trajectories are constructed by combining them with blood glucose sequences, and the corresponding topological closure index is calculated. The maximum value of the locking exponent corresponds to the modal component; S400, the optimization decision module 203 determines whether the maximum topological closure index is greater than or equal to the preset confidence threshold; If so, proceed to step S500; If not, proceed to step S600 and mark it as an open-loop divergence state; S500 and reconstructing module 204 introduce virtual damping coefficients to smooth and regularize the trajectory of the selected optimal modal components and generate a closed hysteresis loop trajectory. The hysteresis loop area is calculated by the feature quantification module 205 to quantify the endocrine regulation loss. S600 When the open-loop divergence state is determined, the forced closure reconstruction is stopped, and the trajectory phase space divergence index is calculated by the feature quantization module 205 to characterize the degree of instability of the regulation system.
[0038] The preprocessing calibration module 201 is mainly responsible for solving the problem of spatiotemporal asynchrony caused by different sampling mechanisms and differences in hardware physical characteristics of multi-source physiological signals. Since continuous blood glucose monitoring sensors are usually implanted in subcutaneous interstitial fluid, while electrocardiogram sensors or pulse wave sensors collect bioelectric or vascular volume signals, the two not only have different sampling frequencies, but also have inherent differences in physical media in the physiological transmission path. Therefore, direct data superposition cannot reflect the real endocrine-neural regulation relationship, and strict spatiotemporal alignment must be performed.
[0039] The preprocessing calibration module 201 receives the raw blood glucose sequence generated by the metabolic acquisition unit 101. and the raw heart rate variability sequence generated by the neural acquisition unit 102 Then, signal synchronization and resampling operations are performed.
[0040] Given that the metabolic acquisition unit 101 typically uses intermittent sampling (e.g., once every 3 to 5 minutes) while the neural acquisition unit 102 provides continuous high-frequency sampling, the system uses the timestamp of the original blood glucose sequence as a reference and performs sliding window averaging on the original heart rate variability sequence.
[0041] The width of the sliding window is set to 1 to 2 times the sampling interval to filter out transient high-frequency noise while preserving neural regulatory trends that match the metabolic timescale.
[0042] Alternatively, the system can use the high-frequency timestamps of the original heart rate variability sequence as a reference and employ a cubic spline interpolation algorithm to sample and reconstruct the low-frequency original blood glucose sequence.
[0043] In this embodiment, the system uniformly establishes a sampling point count of... Discrete time series This ensures that the two sets of signals have a one-to-one index relationship in the time dimension.
[0044] After timing alignment is completed, the preprocessing calibration module 201 performs instrument hysteresis compensation based on derivative cross-correlation. For the inherent physical transmission delay of the interstitial fluid glucose sensor relative to fingertip whole blood measurement, the system needs to separate this non-physiological delay from the actual physiological delay of endocrine regulation.
[0045] This physical delay only causes a shift in the overall signal along the time axis, without altering the waveform characteristics. If the original numerical sequence is used directly to calculate the correlation, the differences in baseline blood glucose levels among individuals will significantly interfere with the accuracy of the correlation coefficient calculation.
[0046] Therefore, this embodiment uses a first-order difference sequence for cross-correlation calculation. Its physical meaning is to compare the synchronicity of the rate of change of the two signals, rather than the similarity of the absolute values, thereby effectively eliminating the influence of individual baseline drift.
[0047] The specific lag time calculation process is as follows: Define the first-order difference sequence of the original blood glucose sequence as The first-order difference sequence of the original heart rate variability sequence is The calculation formula is:
[0048] in The sampling interval is defined as follows: Then, the cross-correlation function of the two difference sequences within a preset physical lag window is calculated. Preset physical lag window The range is set to .
[0049] in To maximize the allowable hysteresis time, this threshold is typically set between 15 and 30 minutes, based on the physical characteristics of existing commercial electrochemical glucose sensors. This upper limit is set to prevent the algorithm from incorrectly aligning two unrelated physiological events (such as eating and subsequent exercise).
[0050] The formula for calculating the cross-correlation function is as follows:
[0051] In the above formula, This represents the number of time displacement steps, and its value range is...
[0052] and These represent the arithmetic mean of the corresponding difference sequences within the calculation window. The preprocessing calibration module 201 traverses the physical lag window. All inside Value, find the value of the cross-correlation function. The maximum displacement is denoted as the instrument's physical lag time. :
[0053] After obtaining this parameter, the module performs time-shift compensation on the original blood glucose sequence to generate a calibrated blood glucose sequence. : ; Through the above steps, the system eliminates the pure time delay caused by the physical characteristics of the sensor, while preserving the physiological relative relationship between the two in terms of numerical response amplitude and waveform evolution.
[0054] The preprocessing calibration module 201 finally performs a normalization operation. To construct a dimensionless unified phase space coordinate system and eliminate the order-of-magnitude difference between blood glucose concentration units (mmol / L) and heart rate variability units (ms), the module uses an extreme value normalization method to normalize the calibrated blood glucose sequence. and heart rate variability sequence All values are mapped to the [0,1] interval. Unlike Z-score normalization (mean 0, variance 1), extremum normalization preserves the relative fluctuation ratio of the original signal and restricts the data to be non-negative. This is beneficial for subsequently constructing hysteresis loops in the first quadrant and calculating their geometric area for any input sequence. The normalization calculation formula is:
[0055] in and These represent the maximum and minimum values of the sequence within the current analysis window, respectively. After this processing, the normalized blood glucose sequence is output. and normalized heart rate variability sequence This serves as the data foundation for subsequent mode decomposition and phase space reconstruction.
[0056] For outliers in a data sequence, those skilled in the art can use median filtering or the Laida criterion to remove them. These are well-known techniques in the field and will not be elaborated here.
[0057] The modality decomposition module 202 is mainly responsible for adaptive separation in the time and frequency domains of the preprocessed heart rate variability sequence to resolve the different physiological regulatory components hidden in the mixed signal.
[0058] Because human heart rate variability signals are regulated by multiple feedback loops of the autonomic nervous system, including respiratory sinus arrhythmia, baroreflex, and thermoregulation and metabolism, these components superimpose each other in the time domain, and their frequencies change dynamically over time, making them typical nonstationary and nonlinear signals.
[0059] Traditional linear filtering methods are based on a fixed cutoff frequency, which can easily lead to the loss of effective signals or waveform distortion. Therefore, this embodiment adopts an ensemble empirical mode decomposition algorithm and uses noise-assisted analysis technology to solve the mode aliasing problem commonly found in traditional empirical mode decomposition. That is, it avoids the same mode component containing extremely different time scales, or signals of the same time scale being separated into different modes.
[0060] The specific processing steps performed by the mode decomposition module 202 include noise-assisted signal construction, intrinsic mode screening, and lumped average calculation.
[0061] Modality decomposition module 202 constructs an auxiliary signal set with added white noise, and sets the ensemble averaging number to 1. It is usually taken as an integer between 100 and 300 to ensure the stability of statistical properties.
[0062] For each decomposition process (in The module generates a sequence of normalized heart rate variability. Equal-length zero-mean Gaussian white noise sequence The standard deviation of the white noise sequence Set as the standard deviation of the original signal times, of which This is the noise amplitude coefficient.
[0063] In this embodiment, the value range is set to 0.1 to 0.2. The module superimposes the generated white noise sequence onto the original signal to form the first... The target signal to be decomposed :
[0064] The purpose of introducing white noise is to take advantage of the uniform distribution of the white noise spectrum to automatically fill the discontinuities in the original signal spectrum, so that signal components of different scales can be continuously decomposed by attaching to the corresponding noise background, thereby avoiding the simultaneous inclusion of extremely different time scale features in a single mode function.
[0065] Mode decomposition module 202 for each target signal When performing empirical mode decomposition, for the extreme point envelope fitting, mean calculation and screening stopping criteria in EMD decomposition, those skilled in the art can use cubic spline interpolation to construct the upper and lower envelopes and use the standard deviation as the screening termination threshold. This is a well-known technique in the field and will not be elaborated here.
[0066] After screening, the first The next decomposition will target signal It is decomposed into a set of intrinsic mode functions and a residual trend term. Let the decomposition result be the first... The intrinsic mode functions are: ,in Represents the modal order. The decomposition process satisfies the following mathematical relationship:
[0067] In the above formula, For the first The residual terms of the second decomposition. Under this decomposition structure... Representing the highest frequency fluctuation component, as the order increases... As the frequency of modal components increases, the frequency of the modal components gradually decreases. This represents the lowest frequency fluctuation component.
[0068] Finally, the mode decomposition module 202 performs a lumped average operation to eliminate the influence of auxiliary noise, due to the added... For zero-mean white noise, according to the law of large numbers, when the number of integrations... When the value is sufficiently large, the algebraic sum of the noise components in the multiple decomposition results will approach zero, and the module will... The arithmetic mean of the eigenmode functions of the corresponding order obtained from the decomposition is used to obtain the final _th_ ... eigenmode functions :
[0069] Similarly, the final residual trend term The modal decomposition module 202 outputs a set of modes strictly ordered from high to low frequency, based on the average value of each residual term after the above processing. :
[0070] Each component of the set Each component represents an independent physiological fluctuation characteristic at a specific time scale, and the components are approximately orthogonal to each other. This step provides a complete library of basis functions for the subsequent precise extraction of specific frequency band signals that are strongly coupled with glucose metabolism from complex neural signals.
[0071] The optimization decision module 203 is the core logic unit of this system. Its main function is to identify and lock onto the neural regulatory component with the strongest causal coupling relationship with blood glucose metabolism fluctuations from the intrinsic mode function set output by the mode decomposition module, using prior physiological knowledge and topological geometry evaluation indicators. This process aims to solve a core technical problem: How to distinguish between effective physiological feedback signals and background random noise in spontaneous endocrine regulation processes lacking synchronous triggering signals.
[0072] The optimization decision module 203 first executes step S301: Initial screening based on physiological spectrum masks. Although the intrinsic mode functions obtained from mode decomposition mathematically satisfy orthogonality, not all components have a clear endocrine regulatory significance.
[0073] The regulation of metabolism by the nervous system is a slow-varying process, typically exhibiting extremely low or ultra-low frequency fluctuations. Therefore, this module introduces the Hilbert transform on the set... Each modal component in Perform analytical signal construction.
[0074] For any component Its Hilbert transform Defined as and The system calculates the instantaneous phase of the analytic signal based on the convolution. Then, by differentiating the phase, the instantaneous frequency is obtained, and the average instantaneous frequency of this component over the entire time path is calculated. .
[0075] The system has a preset effective frequency band range based on endocrine physiology. In the specific application scenario of this embodiment, considering the time constant characteristics of the insulin-glucose regulation loop, the effective frequency band is set to 0.003Hz to 0.04Hz. This frequency band corresponds to the extremely low-frequency components of the autonomic nervous system and is a recognized frequency band reflecting thermoregulation, vasomotor activity, and metabolic activity.
[0076] The optimization decision module 203 will include all average instantaneous frequencies Falling Modal components outside the defined range are considered as unrelated noise or interference caused by physical activity and are therefore discarded. The remaining components that satisfy the frequency constraints constitute the effective modal subset. This process is called spectral masking, which uses physiological prior knowledge to compress the solution space and prevent the algorithm from incorrectly fitting high-frequency respiratory signals as metabolic regulatory signals.
[0077] For effective modal subset For each candidate component, the optimization decision module 203 performs the calculation of the topological closure index. This step aims to quantitatively evaluate the ability of each component to form a closed steady-state loop with the blood glucose sequence in phase space. The module traverses each component in the subset. Use it as the ordinate The corresponding normalized blood glucose sequence As the x-axis Constructing temporary phase space trajectories in a two-dimensional Cartesian coordinate system The trajectory consists of a series of discrete points It is connected in chronological order.
[0078] For each temporary trajectory The module calculates its topological closure index. The calculation formula is as follows:
[0079] The parameters in the above formula are defined as follows: The geometric area represents the convex hull formed by the set of trajectory points. The convex hull is a convex polygon that contains the set of points and has the smallest area. This area characterizes the energy range traversed by the control system in phase space; the larger the area, the more significant the control amplitude corresponding to that component. Those skilled in the art can quickly calculate the vertices and area of the convex hull using Graham's scan method or the monotonic chain algorithm.
[0080] The total path length of the trajectory represents the sum of the Euclidean distances between all adjacent time points. This parameter is used as the denominator to penalize high-frequency oscillations and detours in the trajectory. For physiological regulation, an effective regulation path should be smooth and energy-efficient. An excessively long path length usually means that it contains too much noise and jitter.
[0081] Represents the starting point of the trajectory and the finish line The Euclidean distance between them, also known as the gap distance, directly reflects the degree to which the system returns to its initial state after the observation period ends. The smaller the value, the closer the system is to a perfect closed-loop steady state.
[0082] It is a very small positive constant (e.g.) ), used to prevent when the trajectory is completely closed ( When the denominator is zero, the calculation overflows.
[0083] From the above formula, we can see that the topological closure index... It is a comprehensive indicator: it rewards trajectories with a large adjustment range, smooth evolution trend, and high initial and final closure. The optimization decision module 203 calculates the indicator values for all components in the subset and selects... The largest component is taken as the optimal modal component. The corresponding maximum value is .
[0084] Finally, the optimization decision module 203 performs a confidence decision on the adjustment state. Simple mathematical optimization can only guarantee that the relatively best result is found in this decomposition, but it cannot determine whether the result has clinical steady-state significance.
[0085] If a patient is in a state of severe metabolic breakdown or acute stress, their physiological regulatory system may be in an open-loop disordered state. In this case, even the optimal components cannot form a closed loop. Therefore, the system is set with a confidence threshold. This threshold is obtained by training on historical data from a large number of healthy people and typical diabetic patients. In this embodiment, its empirical value ranges from 2.0 to 5.0.
[0086] The module will calculate the maximum closure value. With threshold Comparison: like The system determines that the current endocrine regulation process is in a steady-state closed state, indicating that the subject's nervous system has successfully completed the negative feedback regulation of blood glucose fluctuations. The system generates a positive flag and triggers the subsequent hysteresis loop reconstruction process.
[0087] like The system determines that the adjustment process is in an open-loop divergent state, indicating that the subject's regulatory system has failed to respond effectively or has failed to complete the regression, and is in an unstable state. At this point, the system generates a negative flag and no longer forcibly performs closed reconstruction, but instead proceeds to the instability feature extraction process.
[0088] The reconstruction module 204 is activated after the optimization decision module 203 determines that the current adjustment state is a steady-state closed state. Its core task is to transform the discrete, noisy and non-strictly closed original signal mapping pair into a smooth, continuous and mathematically closed geometric hysteresis loop, so as to facilitate subsequent energy quantization calculation.
[0089] The specific implementation logic of refactoring module 204 includes the following steps: Step S501: Reconstruction module 204 performs spatial mapping of target data and construction of basic trajectory.
[0090] The system extracts the calibrated blood glucose sequence. As the input state variable of the system, the optimal modal components are extracted after screening. As a feedback response variable of the system.
[0091] By time index The one-to-one correspondence maps two one-dimensional time series to two-dimensional Euclidean space to form the original point set sequence. The original trajectory generated at this time usually has high-frequency spikes caused by measurement noise, and due to the long-term drift characteristics of biological signals, its beginning and end points usually do not coincide, showing an unclosed gap state.
[0092] S502, Reconstruction module 204 calculates the blood glucose change rate sequence as a reference benchmark for dynamic damping.
[0093] Endocrine regulation in living organisms has physical inertia, i.e., damping characteristics. This means that the actual physiological regulation trajectory should be a smooth curve. In order to retain the real fast regulation characteristics while removing noise, this system does not use a low-pass filter with fixed parameters, but introduces the concept of virtual damping.
[0094] The module first calculates the absolute rate of change of the blood glucose sequence with respect to time. This rate of change reflects the current activity level of the metabolic system: when When the value is large, it means that blood sugar is fluctuating wildly and the neural regulatory system is in an active response period. At this time, damping should be reduced to preserve signal details. when When the value is low, it means that blood sugar is in a plateau phase. At this time, the small fluctuations in the signal are mostly invalid noise, and damping should be increased to smooth the trajectory.
[0095] Step S503, Reconstruction module 204 constructs a virtual damping regularization model based on dynamic weights.
[0096] The system defines a time-varying dynamic smoothing coefficient. This coefficient is related to the rate of change in blood glucose. This exhibits a positively correlated nonlinear mapping relationship. The system can implement this mapping using a Sigmoid function, a piecewise linear function, or a hyperbolic tangent function. In this embodiment, a mapping function of the following Sigmoid form is used:
[0097] In the above formula, and These are the lower and upper limits of the smoothing coefficient, respectively, and this range determines the adjustment depth of the filter.
[0098] The method for determining key parameters is as follows: The threshold value is the center of the rate of change, and the value is taken from the entire blood glucose rate of change sequence. The median or average value ensures that the dynamic adjustment mechanism is most sensitive near the statistical center of signal changes.
[0099] To adjust the sensitivity coefficient, its value is usually set to . ,in This represents the standard deviation of the blood glucose rate of change sequence. This setting ensures that the transition band of the Sigmoid function covers the main fluctuation range of the data.
[0100] After obtaining the dynamic coefficients, the module optimizes the modal components. Perform weighted recursive smoothing to generate a normalized sequence of ordinates. :
[0101] This formula achieves time-varying inertial smoothing: when blood sugar is stable, The current value is relatively small, and it mainly depends on historical values, exhibiting large damping characteristics. During drastic changes in blood sugar It is relatively large, and the current value depends mainly on the real-time input, exhibiting a fast following characteristic.
[0102] In step S504, the reconstruction module 204 performs forced closure processing of the trajectory. In order to accurately calculate the area enclosed by the phase space trajectory using Green's formula, the trajectory must form a closed geometric loop. Even after the smoothing process in S503, there may still be small gaps at the beginning and end of the trajectory due to instrument zero drift or ultra-long period fluctuations.
[0103] The module calculates the residuals at the beginning and end of the normalized sequence. Subsequently, the module employs a linear detrending method to eliminate the gap, that is, constructing a curve that linearly increases from 0 to... The compensation sequence is subtracted from the original trajectory to ensure that the coordinates of the beginning and end of the corrected trajectory are strictly equal. The final point set sequence... This constitutes a mathematically closed hysteresis loop. In cases where the loop is determined to be open-loop divergent, the reconstruction module 204 skips the forced closure operation in step S504 and directly outputs the open trajectory.
[0104] The feature quantification module 205 is the final analysis unit of this system. Its main responsibility is to convert the geometric features output by the reconstruction module 204 into numerical indicators with clear physiological / pathological significance. Based on the regulatory state flag generated by the optimization decision module 203, this module uses different mathematical models to quantify and score the neuro-metabolic coupling characteristics of the subjects. The specific implementation logic of the feature quantification module 205 includes the following steps: Step S601: Feature quantization module 205 identifies the adjustment state flag and retrieves the corresponding dataset. When the flag is in steady-state closure, the module receives the set of hysteresis loop points processed by forced closure in S504. ,in Corresponding normalized blood glucose sequence, For the corresponding normalized optimal modal components, when the flag is open-loop divergence, the module receives the open trajectory point set that has been smoothed by S503 but is not closed.
[0105] Step S602: For the steady-state closed state, the feature quantization module 205 performs the calculation of the energy loss of the hysteresis loop. In nonlinear dynamics and control theory, the area enclosed by the phase space closed trajectory represents the energy dissipation or efficiency loss of the system during a complete cycle of excitation response.
[0106] Specifically, in the neural-metabolic coupling model of this invention, the area reflects the degree of phase lag and amplitude gain mismatch between the neural regulatory signal and the change in blood glucose concentration. In this embodiment, the discrete Green's formula is used to numerically calculate the polygon area.
[0107] The calculation formula is described as follows: For including A closed polygon with vertices, its geometric area The absolute value of half the sum of the cross products of all adjacent vertices is the module for the nth vertex in the sequence. The x-coordinate of the point and the th The product of the ordinates of the points, minus the first point The x-coordinate of the point and the th The product of the ordinates of the points, summed over the entire sequence, and the absolute value divided by 2. Mathematically, this can be expressed as:
[0108] Due to input data and All values have been normalized to the [0,1] interval, and the calculated values are... It is a dimensionless value, and must be between 0 and 1.
[0109] The physiological significance of this indicator is explained as follows: The larger the value, the wider the hysteresis loop. This means that there is a significant phase lag between the neural regulatory signal and the change in blood glucose, or that the nervous system needs to maintain a high intensity of firing for a longer period of time to pull blood glucose back to baseline. Clinically, this usually corresponds to insulin resistance or autonomic nervous system sluggishness, that is, the system has low regulatory efficiency and high energy consumption.
[0110] The smaller the value, the narrower and longer the hysteresis loop, even approaching a straight line. This means that neural regulation is highly synchronized with changes in blood glucose, and the system's response is rapid and precise. Clinically, this usually corresponds to good insulin sensitivity and sound autonomic nervous system function.
[0111] Step S603: For the open-loop divergence state, the feature quantization module 205 calculates the phase space divergence index. In this state, since the trajectory is not closed, area calculation is no longer applicable. The system needs to quantify the degree of loss of control of the subject's metabolic regulation function. The module mainly calculates the drift distance of the trajectory endpoint relative to the starting point and normalizes it in combination with the trajectory's expansion range.
[0112] Divergence Index The calculation logic is as follows: First, calculate the starting point of the trajectory. and the finish line Euclidean distance between This distance directly reflects the residual error in the control system's failure to pull the state back to the baseline.
[0113] Secondly, calculate the maximum characteristic scale of the trajectory in phase space. To avoid scale deviation in one direction, this embodiment uses the diagonal length of the bounding box of the trajectory point set as the normalized denominator, and the vertex coordinates of the bounding box are determined by the maximum / minimum horizontal and vertical coordinates in the sequence.
[0114] The final formula for calculating the divergence index is:
[0115] The index The value range is typically between [0, 1.414]. The risk threshold for this indicator is determined through subject operating characteristic curve analysis, with the value corresponding to the maximum point of the Youden index selected as the critical value. like A value >0.3 indicates a significant regulatory failure, suggesting that the subject may be experiencing acute stress, hypoglycemic rebound, or drug failure.
[0116] If 0.1 < A value ≤0.3 indicates that the subject is in a metastable state and may have mild autonomic dysfunction. Although the regulation did not completely fail, precise regression was not achieved.
[0117] Step S604, Feature Quantization Module 205 converts the raw quantization indexes obtained above into... The data is mapped to a clinically readable comprehensive score, and the system has a pre-stored percentile lookup table based on the distribution of a large sample population.
[0118] The module will calculate The value is compared with the reference distribution of healthy people of the same age and gender, and a neuro-metabolic coupling health score between 0 and 100 is output. For open-loop divergent samples, the system directly marks them as high-risk level and outputs the divergence index as an auxiliary diagnostic reference.
[0119] These quantitative results are ultimately presented to users through a display interface or stored in a database for long-term chronic disease management and tracking. For the specific percentile mapping algorithm and database construction method, those skilled in the art can implement them based on conventional medical statistical methods, which are well-known technologies in the field and will not be elaborated here.
[0120] Example: To further illustrate the implementation details and application value of this invention, the following explanation is based on a typical clinical application scenario.
[0121] 1. Subjects and Data Collection A 56-year-old male subject was selected, clinically diagnosed with early-stage type 2 diabetes mellitus, accompanied by mild autonomic dysfunction. Data was collected using this system during the subject's standardized mixed meal tolerance test.
[0122] Metabolic data: Using an implantable continuous glucose monitor (CGM), data was collected every 5 minutes for a total of 4 hours to obtain raw blood glucose sequences. ,length point.
[0123] Neural data: A single-lead ECG recorder was used, with a sampling frequency of 250 Hz. The system extracted the RR interval sequence, resampled it, and aligned it with the blood glucose sequence to obtain the raw heart rate variability sequence. .
[0124] 2. Preprocessing and Timing Compensation Effects The preprocessing calibration module 201 first normalizes the two sets of sequences. Calculations show that without time compensation, the direct Pearson correlation coefficient between the original blood glucose sequence and the heart rate variability sequence is only 0.42, and the waveforms are significantly misaligned on the time axis.
[0125] The module performs first-order difference cross-correlation calculations and searches within a preset lag window of [0, 30] minutes.
[0126] Calculation results: in lag time At 3 minutes (i.e., 3 sampling points), the cross-correlation function It has reached its peak.
[0127] Correction operation: The system performs a 15-minute time shift compensation on the blood glucose sequence.
[0128] Correction effect: After correction, the correlation between the two sets of sequences in terms of physiological trend improved to 0.89. This step ensures that subsequent analysis is based on synchronous data from the same physiological event (metabolic-neurological response caused by eating).
[0129] 3. Modal decomposition and automatic optimization The mode decomposition module 202 performs ensemble empirical mode decomposition on the corrected heart rate variability sequence and sets the noise amplitude coefficient. Number of integrations The decomposition generated seven intrinsic mode functions (IMF1-IMF7) and one residual term.
[0130] The optimization decision module 203 first performs initial screening based on the physiological frequency band (0.003Hz-0.04Hz), retaining three candidate components: IMF3, IMF4, and IMF5. Then, the topological closure index (ITC1) of each component is calculated, and the results are shown in the table below:
[0131] The module locks IMF4 as the optimal modal component, and its ITCI value (4.87) is greater than the preset confidence threshold (2.5). The system determines that the adjustment state is a steady-state closure and triggers the subsequent reconstruction process.
[0132] 4. The closure reconstruction and feature quantization reconstruction module 204 introduces a virtual damping model based on the blood glucose change rate weight to perform trajectory normalization on IMF4 and execute forced closure. Finally, the hysteresis loop area is calculated by the feature quantization module 205. .
[0133] Clinical interpretation: This value (0.42) is significantly higher than the healthy baseline value (usually <0.2), indicating that although the subject was able to complete the closed loop of blood glucose regulation, there was a significant phase lag in neural regulation relative to blood glucose changes, suggesting early insulin resistance or delayed neural conduction.
[0134] Experimental verification and effect comparison To verify the effectiveness and advancement of the method proposed in this invention, the applicant conducted comparative experiments on simulated datasets and real clinical datasets.
[0135] Experiment 1: Robustness Verification of the Cross-Correlation Timing Compensation Algorithm Experimental objective: To verify the advantages of the first-order differential cross-correlation algorithm used in this invention over the traditional direct cross-correlation algorithm in resisting baseline drift.
[0136] Experimental setup: A set of simulated coupled signals with a known delay (10 minutes) was constructed, and a significant linear baseline drift (simulating fasting blood glucose elevation) was added to the blood glucose signal.
[0137] Comparison results:
[0138] Conclusion: The preprocessing calibration scheme proposed in this invention can effectively resist baseline interference caused by individual differences and ensure the physical authenticity of timing alignment.
[0139] Experiment 2: Verification of the effectiveness of modal decomposition and optimization strategies Experimental objective: To verify the performance of ensemble empirical mode decomposition (EEMD) combined with topological closure index (ITCI) in solving mode aliasing and automatic feature extraction.
[0140] As attached Figure 3 As shown in the image, the sub-image above is from the traditional EMD method: the waveform shows severe modal aliasing, with the main metabolic regulatory components being broken down into two components, IMF3 and IMF4, accompanied by significant intermittent high-frequency noise. It is difficult for a person to determine which component represents the true neural regulation. The following figure illustrates the EEMD+ITCI method of this invention: the waveform shows that the selected optimal mode component (IMF4) has a stable oscillation period and a smooth waveform. This component accounts for more than 85% of the energy within the physiological frequency band and has the highest coupling degree with blood glucose fluctuations.
[0141] Data support: In a test on 50 samples, the consistency between the automatic optimization results of this invention and the results of expert manual annotation reached 94%, while the consistency of the traditional selection method based on energy ratio was only 72%.
[0142] Experiment 3: Improvement of Trajectory Quality by Dynamic Damping Reconstruction Experimental objective: To verify the denoising and fidelity preservation effect of the virtual damping model based on blood glucose change rate weights in phase space reconstruction.
[0143] As attached Figure 4 As shown, the dashed trajectory represents linear filtering: a low-pass filter with a fixed cutoff frequency. At the turning point where blood glucose levels change drastically (i.e., the tip of the hysteresis loop), the trajectory exhibits significant rounded corner distortion, leading to an underestimation of the hysteresis loop area and loss of peak regulation information.
[0144] The solid line trajectory represents the dynamic damping of this invention: the trajectory is smooth and burr-free in flat areas; In the transition region (where blood glucose changes rapidly), the damping automatically decreases, and the trajectory closely follows the original data points, preserving the sharp transition characteristics.
[0145] Conclusion: The dynamic damping model resolves the contradiction between denoising and preserving pathological features, making the generated hysteresis loop more realistically reflect the nonlinear regulation process.
[0146] Experiment 4: Clinical Discrimination of Hysteresis Loop Area Index Experimental objective: To verify quantitative indicators The statistical significance of (hysteresis loop area) in distinguishing different healthy populations.
[0147] As attached Figure 5 As shown, the dataset is: Healthy control group (Health): 30 cases with normal glucose tolerance.
[0148] Pre-diabetes group (Pre-DM): 30 cases with impaired fasting glucose.
[0149] Type 2 diabetes mellitus group (T2DM): 30 cases, diagnosed and with a disease duration of >5 years.
[0150] Statistical results:
[0151] Conclusion: With the deterioration of metabolic regulation function, the area of the hysteresis loop increases. It exhibits a monotonically increasing trend. This indicator can not only distinguish between healthy and diseased groups, but also effectively differentiate between prediabetes and the diagnosed stage, confirming its clinical value as a tool for quantifying endocrine regulation depletion.
[0152] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.
Claims
1. A multi-parameter physiological signal fusion analysis system for endocrinology, characterized in that, include: The data acquisition terminal is configured to acquire the subject's raw blood glucose sequence and heart rate variability sequence in real time; A processing server is communicatively connected to the data acquisition terminal; The output terminal is used to display the trajectory and quantitative indicators; The processing server includes: The preprocessing calibration module is used to perform normalization on the raw blood glucose sequence and heart rate variability sequence, and calculate the physical lag time between the two sets of sequences for time series compensation. The mode decomposition module is used to perform ensemble empirical mode decomposition on heart rate variability sequences to generate a set of eigenmode functions; The optimization decision module is used to screen effective modal components, combine them with blood glucose sequences to construct temporary phase space trajectories, select the optimal modal component based on the topological closure index, and determine whether the regulation state is a steady-state closed state or an open-loop divergent state. The reconstruction module is used to construct the target phase space based on the optimal modal components and the blood glucose sequence, and to normalize the optimal modal components using a virtual damping model based on the weight of the blood glucose change rate. The feature quantization module is used to calculate the hysteresis loop area or phase space divergence index of the trajectory in the target phase space according to the adjustment state.
2. The endocrinology multi-parameter physiological signal fusion analysis system according to claim 1, characterized in that, The preprocessing calibration module calculates the physical lag time in the following ways: Calculate the first-order difference sequence between the original blood glucose sequence and the heart rate variability sequence; Calculate the cross-correlation function of two first-order difference sequences; The time delay corresponding to the peak value of the cross-correlation function is selected as the instrument transmission lag time, and time displacement compensation is performed on the original blood glucose sequence accordingly.
3. The endocrinology multi-parameter physiological signal fusion analysis system according to claim 1, characterized in that, The specific methods by which the mode decomposition module performs set empirical mode decomposition include: A white noise sequence was added to the calibrated heart rate variability sequence, and empirical mode decomposition was performed to obtain preliminary modal components. Repeat the white noise sequence addition and decomposition steps multiple times; The initial modal components obtained from multiple decompositions are averaged to output a set of intrinsic mode functions that can suppress modal aliasing.
4. The endocrinology multi-parameter physiological signal fusion analysis system according to claim 1, characterized in that, The optimization decision module selects the optimal modal component in the following ways: The instantaneous frequency of each intrinsic mode function is calculated using Hilbert transform, and effective mode components whose dominant frequency falls within the preset physiological frequency band are retained. Calculate the topological closure index of each temporary phase space trajectory, wherein the topological closure index characterizes the proximity of the beginning and end points of the trajectory and the smoothness of the trajectory; The effective modal component corresponding to the maximum value of the topological closure index is selected as the optimal modal component.
5. The endocrinology multi-parameter physiological signal fusion analysis system according to claim 4, characterized in that, The optimization decision module determines the adjustment state in the following ways: The selected maximum topological closure index is compared with the preset confidence threshold; If the maximum topological closure index is greater than or equal to the confidence threshold, the adjustment state is determined to be a steady-state closure state. If the maximum topological closure index is less than the confidence threshold, the adjustment state is determined to be an open-loop divergent state.
6. The endocrinology multi-parameter physiological signal fusion analysis system according to claim 1, characterized in that, The virtual damping model in the reconstruction module uses a dynamic smoothing coefficient to recursively smooth the optimal modal components. The dynamic smoothing coefficient is positively correlated with the absolute rate of change of the original blood glucose sequence: when the absolute rate of change increases, the dynamic smoothing coefficient increases to reduce the damping effect. When the absolute rate of change decreases, the dynamic smoothing coefficient decreases to enhance the damping effect.
7. The endocrinology multi-parameter physiological signal fusion analysis system according to claim 6, characterized in that, The reconstruction module is also configured to perform forced trajectory closure processing: When the regulation state is a steady-state closed state, calculate the first and last residuals of the optimal modal component sequence after normalization by the virtual damping model; A linear detrending method is used to distribute the first and last residuals to each point in the sequence, so that the final generated phase space trajectory is geometrically closed.
8. The endocrinology multi-parameter physiological signal fusion analysis system according to claim 1, characterized in that, When the feature quantization module determines that the state is a steady-state closed loop, it uses the discrete Green's formula to calculate the area of the hysteresis loop: Calculate the difference of the cross products of adjacent vertices in the phase space trajectory, sum all the differences, and take half of the absolute value. The area of the hysteresis loop characterizes the energy loss and phase lag degree in the process of the nervous system regulating the metabolic system.
9. The endocrinology multi-parameter physiological signal fusion analysis system according to claim 1, characterized in that, When the feature quantization module determines that the state is open-loop divergence, it calculates the phase space divergence index: Calculate the Euclidean distance between the start and end points of the phase space trajectory; Calculate the diagonal length of the phase space trajectory bounding box; The ratio of the Euclidean distance to the diagonal length is used as the phase space divergence index to characterize the degree of instability of the control system.
10. A method for multi-parameter physiological signal fusion analysis in endocrinology, characterized in that, The method of using the multi-parameter physiological signal fusion analysis system for endocrinology according to any one of claims 1-9 includes the following steps: Obtain the raw blood glucose and heart rate variability sequences of the subjects, calculate the physical lag time between the two sets of sequences, and perform time-series correction and normalization; The ensemble empirical mode decomposition algorithm is used to decompose the corrected heart rate variability sequence to generate an eigenmode function set; Effective modal components are screened and temporary phase space trajectories are constructed by combining blood glucose sequences. The topological closure index is calculated, and the optimal modal components are identified and the regulatory state is determined accordingly. If the state is determined to be a steady-state closed state, a virtual damping model based on the weight of blood glucose change rate is introduced to normalize the optimal mode components, and a forced closure process is performed to generate a closed hysteresis loop, and the area of the hysteresis loop is calculated. If the open-loop divergence state is determined, the forced closure process is stopped, and the normalized drift distance of the first and last points of the phase space trajectory is calculated to obtain the phase space divergence index.