Noninvasive continuous blood pressure prediction method based on multi-dimensional data driving

Abstract

The invention discloses a non-invasive continuous blood pressure prediction method based on multi-dimensional data driving, which comprises the following steps: reading ECG and PPG original signals, and preprocessing the ECG and PPG original signals; wherein the preprocessing comprises filtering, polarity correction, time alignment and normalization preprocessing; time sequence features, amplitude features and waveform morphological features are obtained through calculation based on the preprocessed ECG and PPG signals, and a feature set is constructed in combination with physiological parameters; selecting a deep learning framework, constructing a Transform hybrid model, and designing a loss function; and selecting an optimization function amount learning rate attenuation strategy, setting hyper-parameters including the number of iterations and the batch size, inputting the feature set into the network for training, and testing the model. The electrocardiogram signals, the photoplethysmography signals and the physiological parameters are fused, blood pressure prediction is achieved through an improved Transform architecture, and the adaptability and stability of the model are enhanced while the prediction accuracy is improved.

Description

Technical Field

[0001] This invention relates to the field of machine learning technology, specifically to a non-invasive continuous blood pressure prediction method based on multidimensional data-driven approaches. Background Technology

[0002] Hypertension is one of the leading risk factors for serious health problems such as cardiovascular disease and stroke worldwide. Traditional blood pressure measurement relies on cuff-based methods (such as auscultation and oscillometric methods), which require inflatable cuffs and can only provide intermittent readings, making continuous monitoring difficult. Although methods such as arterial pressure volume clamping can provide continuous monitoring, the equipment is bulky and inconvenient to use, hindering widespread application. Non-invasive continuous blood pressure monitoring technology, due to its non-intrusive and continuous measurement characteristics, holds the promise of revolutionizing traditional methods. In particular, research on estimating blood pressure using vital signs such as electrocardiogram (ECG) and photoplethysmography (PPG) signals has attracted attention. Among these, pulse transit time (PTT), the time difference between cardiac electrical activity and peripheral pulse waves, is negatively correlated with blood pressure and has been proven by numerous studies to be an important indicator for cuffless blood pressure estimation.

[0003] However, relying solely on PTT (Pulse Test Tolerance) for blood pressure estimation has limited accuracy and is susceptible to individual differences and environmental factors, requiring calibration compensation. To improve model accuracy, recent studies have introduced other physiological features as aids. Ding et al. proposed the PPG intensity ratio (PIR) as a new indicator, combining it with PTT to significantly improve the accuracy of non-invasive blood pressure estimation. Studies have shown that introducing additional indicators besides PTT helps improve the accuracy of cuffless blood pressure measurement and facilitates better tracking of dynamic changes in blood pressure. Meanwhile, machine learning and deep learning methods have also been applied to non-invasive blood pressure prediction. For example, some studies have comprehensively extracted physiological features and signal statistical features, using regression algorithms to estimate blood pressure, achieving accuracy conforming to the AAMI standard on large-scale data. Furthermore, deep models such as convolutional neural networks and recurrent neural networks have also been used to directly learn features from PPG or ECG+PPG signals for blood pressure prediction, achieving progressively improving performance.

[0004] Despite the success of the above methods, challenges remain: the models have poor generalization ability across different populations and states, insufficient robustness to motion disturbances, and a lack of interpretability. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a non-invasive continuous blood pressure prediction method based on multidimensional data-driven approaches. This method integrates electrocardiogram (ECG) signals, photoplethysmography (PPG) signals, and physiological parameters, and achieves blood pressure prediction through an improved Transformer architecture. This method enhances the adaptability and stability of the model while improving prediction accuracy.

[0006] This invention discloses a non-invasive continuous blood pressure prediction method based on multidimensional data, comprising: Step 1: Read the raw ECG and PPG signals and perform preprocessing; wherein, the preprocessing includes filtering, polarity correction, time alignment and normalization preprocessing; Step 2: Calculate the timing features, amplitude features, and waveform morphology features based on the preprocessed ECG and PPG signals, and construct a feature set by combining physiological parameters; Step 3: Select a deep learning framework, build a Transformer hybrid model, and design a loss function; Step 4: Select the optimization function and learning rate decay strategy, set hyperparameters including the number of iterations and batch size, input the feature set into the network for training and test the model.

[0007] As a further improvement of the present invention, step 1 specifically includes: Step 1.1: Perform second-order Butterworth high-pass filtering and zero-phase band-pass filtering on the original ECG and PPG signals to obtain signals with noise and baseline drift removed for relocation; Step 1.2: After filtering, ensure the peak and valley directions of the PPG signal are consistent through waveform polarity self-check; Step 1.3: After polarity correction, the global micro-delay between ECG and PPG is estimated using the cross-correlation method, and the position of the R peak is corrected.

[0008] As a further improvement of the present invention, step 1.3 specifically includes: Step 1.3.1, sampling rate (For an ECG signal with a maximum frequency of 45Hz, theoretically only a 90Hz sampling rate is needed. 125Hz is much higher than twice the maximum frequency, providing a safety margin for the anti-aliasing filter and the high-frequency components of the signal.) Let the sequence of the ECG and PPG signals be... and n=0,1,...,N-1; the sequence after bandpass filtering is and The existing R-peak index set is The goal is to estimate the lag of the entire channel. And obtain the refined R peak ; Step 1.3.2: Select a time period near the midpoint of the record. Second window (Within a 30-second window, there are approximately 30 to 50 heartbeats. This provides a large number of data points for cross-correlation calculation, which is sufficient to obtain a smooth and stable cross-correlation function curve, thereby reliably estimating the global delay. The window parameters are shown in equation (1).) (1) In the formula, Sampling rate, a and b are the start and end indices of the window used for cross-correlation calculation, M is the number of valid samples in the window, N is the total length of the signal, and L is the window length; Step 1.3.3, construct using equation (2) and Equal length; (2) In the formula, The first-order difference signal (approximate derivative) of the PPG signal is used for cross-correlation calculation with the ECG signal; Step 1.3.4: Using equations (3) and (4) in the window obtained in step 1.3.2 Top and Perform z-score standardization; (3) (4) In the formula, , These are the sample means of the ECG signal and the PPG differential signal within the window, respectively. Let be the sample standard deviations of the ECG signal and the PPG differential signal within the window, respectively. The ECG and PPG differential signals are Z-Score standardized to have a mean of 0 and a standard deviation of 1. Step 1.3.5: Let the discrete time delay index be... Search scope To ensure that the summation does not exceed the bounds, define... The aligned sequence is obtained through equations (5) and (6). Perform discrete cross-correlation and lag estimation; (5) (6) In the formula, This is a time delay index, representing the relative sample point offset between signals. , For a given delay Effective summation range indexes when calculating cross-correlation to prevent array out-of-bounds errors. Let be the cross-correlation value at time delay l. To make the cross-correlation value To achieve the maximum optimal latency index. For optimal latency; Step 1.3.6: Arrange the R peak as a whole according to formula (7). Translate and crop to the legal area; (7) In the formula, To be based on global latency The index of the i-th R-peak after preliminary translation correction; Step 1.3.7: Locally relock the R peak on the ECG wave using equation (8); (8).

[0009] In the formula, To the position of the R peak in the initial correction The radius of the window for performing a fine-grained local search in the vicinity. For the start and end indices of the local fine-grained search window, For the index of the i-th R-peak after final precise location, the window [ Internal ECG filter signal The maximum value point.

[0010] As a further improvement of the present invention, step 2 specifically includes: Step 2.1: Relocate the R peak of ECG using the method described in Step 1.3. Within a certain time window after the R peak, perform local upsampling on the rising segment of PPG to obtain the upsampled sequence. ;exist Find the point of maximum slope Its slope is denoted as ; A baseline of "valley value priority, anomaly backoff" is adopted: the amplitude of the previous valley value within the window is taken first, and if there is no valley value or an anomaly, the median of the preset upsampling start time is used; then let The original sampling index of the left endpoint of the local window is used as the basis for obtaining the PPG footpoint using equation (9). ,in The original sampling index of the R peak is used, and the PTT is finally calculated by equation (10); (9) (10) In the formula, This is the start time index for the local search window used to locate the foot point. This is the factor for local upsampling; in the actual code, it's 8. Let y(t) be the index of the point with the maximum slope in the upsampled signal y(t). This represents the maximum slope value of the rising segment of PPG. To determine the baseline amplitude of the foot point, To calculate the PPG foot time, The time index is the R peak, and PTT is the pulse wave transit time, which is the time difference from the R peak of the ECG to the PPG foot point. Step 2.2: Calculate the time from the PPG wave base to the main peak using equation (11). ,in This is the index of the PPG contraction peak in the original sampling coordinate system; (11) In the formula, For the time index of the main PPG peak (P1 peak), The rise time of the PPG wave; Step 2.3: Calculate the maximum upward slope of PPG using equation (12). That is, the maximum first derivative value of the rising segment of PPG, where PPG after bandpass and smoothing; (12) In the formula, This is the maximum first derivative of the rising segment of the PPG, multiplied by the sampling rate. To correct for amplitude scaling caused by the difference; Step 2.4: Calculate the ratio PIR of the AC component to the DC component of the pulse wave using equation (13), which is the ratio of the AC component to the DC component of the pulse wave; where Peak amplitude, It is the range of its previous trough value; For pulse wave alternating component, , is the DC component of the pulse wave; (13) In the formula, This represents the amplitude of the main peak of PPG. The amplitude value of the preceding trough paired with the main peak. The alternating component of the pulse wave. PIR represents the DC component of the pulse wave, and PIR represents the pulse wave intensity ratio. Step 2.5: Calculate the heart rate using formula (14), which is the time interval between two consecutive R waves; (14) In the formula, The i-th RR interval is the time interval between two consecutive R peaks. Heart rate; Step 2.6: Detect the p2 peak and retrieve the original sampling index. ,time In the fragment Above, let the smoothed amplitude be denoted as . ,in To be the main peak, The valley before the main peak. The next valley after the main peak; in the descending segment after the PPG main peak, the global minimum point is first located using equation (15). Then, all local maxima are found behind the notch using equation (16). As candidate P2 peaks, the candidate set is: Then, index all candidate indices. Convert to time And as in Equation (17), the P2 peak is screened by curvature verification, and finally the optimal candidate point is located by parabolic interpolation using Equation (18). (15) (16) (17) (18) In the formula, These are sampling indices representing the main peak, the valley before the main peak, and the valley after the main peak, respectively. This is the index of the notch location found in the PPG descent segment between the main peak and the next trough, i.e., the global minimum point of the signal in that segment. The set of all local maxima found after the notch is taken as the candidate P2 peak. The second derivative of the PPG signal. , This represents the time corresponding to the notch and the candidate peak of P2. , This represents the original index of the optimal P2 candidate peak and the signal amplitude at that point. This is the sub-pixel offset obtained through parabolic interpolation. The index of the P2 peak after subpixel-level precise positioning; Step 2.7: Calculate the reflectance index RI using equation (19), where It is the main peak of this auction. Index of the previous valley value of the pair Conversion time; (19) In the formula, The reflection index reflects the intensity of the reflected wave of the waveform. The amplitude of the main peak (P1), The amplitude of the reflection peak (P2) The time of the previous trough value paired with the main peak of this auction; Step 2.8: Calculate the systolic-diastolic period using equation (20). ; (20) In the formula, The systolic-diastolic phase; Step 2.9: Combine the obtained feature set with the corresponding subject's age and weight to jointly construct a physiological parameter feature set.

[0011] As a further improvement of the present invention, step 3 specifically includes: Step 3.1: Select PyTorch as the deep learning framework and construct a fusion network. The network consists of a feature extraction module, the "Multi-Source Temporal Encoder Module" (MSEM), and a fusion-decoding module, the "Dynamic Conditional Fusion-Decoder Layer" (DCFD). The MSEM consists of an encoder and a middle layer. The encoder uses a 6-layer Transformer Encoder (128 hidden dimensions, 8 attention heads) to integrate ECG and PPG multi-channel temporal information and enhance long-range dependencies and noise resistance. The middle layer contains two types of sub-blocks: a branch generator and a Feature-wise Linear Modulation (FiLM) parameter generator. The branch generator incorporates three types of lightweight quantum heads: the Main head (a linear regression head after the encoder, with small parameters and stable convergence), the KAN head (Kolmogorov–Arnold network, with analytical approximation capabilities, suitable for nonlinear explicit features), and the eXtreme Gradient Boosting (XGB) head (gradient boosting tree ensemble, adept at tabular and statistical features, robust to small samples and outliers), and extends the residual correction head for error compensation and drift correction. The FiLM parameter generator is a two-layer perceptron structure, outputting the scale and offset coefficients per token for conditional recalibration, achieving individualization and domain adaptation. The Dynamic Conditional Fusion-Decoder Layer (DCFD) consists of a FiLM conditional modulation unit, a branch attention projection unit, a multi-head differential attention (MHDA) unit, a Transformer decoder, a token weighted aggregation unit, and a baseline-gated residual merging unit. The FiLM unit provides lightweight affine modulation to suppress inter-domain differences, while the projection unit unifies the representations of each branch into a comparable subspace. MHDA uses differential enhancement of dual-path similarity to highlight branches and reduce redundancy. The Transformer Decoder uses cross-attention to align branch representations with the main temporal context to improve temporal consistency. Token-weighted aggregation uses learnable weights to achieve interpretable fusion. Baseline gating uses the XGB output as a stable baseline, combined with sigmoid gating and a training-phase ramp strategy to achieve residual superposition and error correction, ultimately forming a dual-channel output of SBP and DBP. Step 3.2, set the input window as... Their window mean / variance are respectively Construct two types of flat features: the original expanded ⊕ statistics for input to XGB and residual tree. ) and the standardized expansion ⊕statistic for FiLM and KAN inputs ( ),in ; Step 3.3, will After embedding and position encoding, the input to the Transformer encoder yields the temporal implicit representation. Then, the intermediate branch generator produces three tokens in parallel. The Main branch takes... Terminal correlation step, via linear have to The KAN branch passes through get The XGB branch is passed through get The residual branch is obtained through the residual tree in the k-th round during the residual training process. Meanwhile, the FiLM parameter generator is based on Output the scale and offset parameters that match the current token count S. ; In the formula, For batch size, The time step length of the input sequence. The number of channels for the input feature. For the input tensor of the model, The mean and standard deviation of the input window X; Step 3.4: Dynamically pre-scale the weights of the three basic branches Main / KAN / XGB using equation (21), and reset the weights of the remaining residual tokens to 1. (twenty one) In the formula, The learnable scalar weights represent the importance of the three branches: Main, KAN, and XGB. To The branch weights after Softmax normalization of the three weights. The original token generated for the i-th branch. The token is pre-scaled with weights. Step 3.5: Conditional modulation is performed using FiLM. An affine transformation is applied to each token using equation (22) to achieve token-by-token feature reweighting; where, This represents the three basic branches and the original residual token after prescaling according to equation (21); (twenty two) In the formula, The original scale and offset parameters generated by MLP, , These are the scale and offset parameters that ultimately apply to the token after transformation. The token is FiLM modulated; Step 3.6: Each branch enters the multi-head differential attention module, first linearly increasing the dimensionality. Stacked Then, the similarity between the two paths is calculated at each attention head using equation (23), and differential reinforcement is performed. (twenty three) In the formula, V represents the Query, Key, and Value matrices, respectively, for two different attention paths. , The attention weight matrix calculated for the two paths. This is a weighting coefficient used to control the influence of the second attention path A2 in the difference process. This is the attention weight matrix after differential enhancement; Step 3.7: Each branch enters the Decoder and is associated across time-series decoders using equation (24); (twenty four) Step 3.8: Weight and aggregate each Token using equation (25) and map it to a 2D dimension, then obtain the result. Copy and align to the predicted step size get ; (25) Step 3.9: Integrate the residual round summation using equations (26) and (27) and merge it with the baseline gated residual, where For the first The learnable weights of the wheel residuals, For learnable gating, This refers to the warm-up coefficient during training. (26) (27) In the formula, The learnable ensemble weights for the residual tokens in the k-th round. For the residual token generated by the residual tree in the kth round, The signal is a learnable gated signal, constrained to the range (0,1) by the Sigmoid function. This is a warm-up coefficient during training, typically increasing linearly from 0 to 1, used to stabilize the initial stages of training. This is the final output of the model; Step 3.10, record the first One sample in the channel The truth value of is Predicted as The error is denoted as Let U be the total number of elements after batch and time step expansion, and let the k-th channel be... Statistics are Total Normalized weight definition Each channel also learns a log-variance parameter. l c Used for uncertainty estimation; the total loss function is designed based on the above, as follows: The total loss function is shown in Equation (28), which includes Gaussian negative log-likelihood, focal regression loss, quantile loss, Huber loss, KAN regularization, log-variance regularization and residual consistency. (28) The Gaussian negative log-likelihood loss is shown in Equation (29) and is used to learn the mean and variance simultaneously. (29) The focus regression loss is shown in Equation (30), emphasizing large errors and difficult samples; (30) The quantile loss is shown in Equation (31), which directly constrains the upper and lower quantiles, allowing the model to learn to predict "blood pressure range" rather than a single value, and is insensitive to abnormal systolic or diastolic readings, thus improving interval prediction and robustness. (31) The Huber loss is shown in Equation (32), which is quadratic for small errors and linear for large errors, taking into account both robustness and differentiability. (32) KAN regularization and logarithmic variance are shown in equations (33) and (34), respectively, which suppress overfitting and stabilize variance; (33) (34) The residual consistency is shown in Equation (35), which makes the residuals close to the true residuals when there is a baseline or residual branch; (35).

[0012] In formulas 28-35, For sample index, For output channel index, This represents the true blood pressure value of the u-th sample in the c-th channel. Let u be the predicted blood pressure value of the u-th sample in the c-th channel. For prediction error, The total number of elements in a batch, including samples and time periods. The log-variance parameter learned for the c-th channel is used for uncertainty estimation. Let be the L1 loss statistic for the k-th channel.

[0013] The normalized weight for the k-th channel. , , , , , , The weight hyperparameters for various loss terms, The target quantile for quantile loss is... The threshold parameter for Huber loss determines the switching point between secondary and primary losses.

[0014] As a further improvement of the present invention, step 4 specifically includes: Choose the Lookahead-packaged RAdam optimizer, select a learning rate scheduling strategy of 10 warm-up rounds (Linear Learning Rate, LinearLR) and cosine annealing, set hyperparameters including the number of iterations and batch size, input the feature set into the network for training, save the model for testing, and evaluate the model performance using the BHS standard and AAMI standard.

[0015] Compared with the prior art, the beneficial effects of the present invention are as follows: This invention integrates electrocardiogram (ECG) signals, photoplethysmography (PPG) signals, and physiological parameters to achieve blood pressure prediction through an improved Transformer architecture, which enhances the model's adaptability and stability while improving prediction accuracy. Attached Figure Description

[0016] Figure 1 Bland-Altman consistency analysis plot for DBP (diastolic blood pressure); Figure 2 Bland-Altman concordance analysis plot for SBP (systolic blood pressure); Figure 3 The error distribution histogram and kernel density curve for DBP (diastolic blood pressure); Figure 4 The kernel density curve of the error distribution histogram of SBP (systolic blood pressure); Figure 5 A scatter plot of DBP (diastolic blood pressure); Figure 6 A scatter plot of SBP (systolic blood pressure); Figure 7 Time-series line graphs of DBP (diastolic blood pressure) and SBP (systolic blood pressure); Figure 8 A comparison table for model evaluation; Figure 9 This is a flowchart of the method of the present invention; Figure 10 This is a diagram of the deep learning framework of the present invention. Detailed Implementation

[0017] 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, 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.

[0018] The present invention will now be described in further detail with reference to the accompanying drawings: like Figure 9 As shown, this invention discloses a non-invasive continuous blood pressure prediction method based on multidimensional data-driven methods, including: Step 1: Read the raw ECG and PPG signals and perform preprocessing; wherein, the preprocessing includes filtering, polarity correction, time alignment and normalization preprocessing; Specifically, it includes: Step 1.1: Perform second-order Butterworth high-pass filtering and zero-phase band-pass filtering on the original ECG and PPG signals to obtain signals with noise and baseline drift removed for relocation; Step 1.2: After filtering, ensure the peak and valley directions of the PPG signal are consistent through waveform polarity self-check; Step 1.3: After polarity correction, the global micro-delay between ECG and PPG is estimated using the cross-correlation method, and the position of the R peak is corrected.

[0019] Specifically, step 1.3 includes: Step 1.3.1, sampling rate (For an ECG signal with a maximum frequency of 45Hz, theoretically only a 90Hz sampling rate is needed. 125Hz is much higher than twice the maximum frequency, providing a safety margin for the anti-aliasing filter and the high-frequency components of the signal.) Let the sequence of the ECG and PPG signals be... and n=0,1,...,N-1; the sequence after bandpass filtering is and The existing R-peak index set is The goal is to estimate the lag of the entire channel. And obtain the refined R peak ; Step 1.3.2: Select a time period near the midpoint of the record. Second window (Within a 30-second window, there are approximately 30 to 50 heartbeats. This provides a large number of data points for cross-correlation calculation, which is sufficient to obtain a smooth and stable cross-correlation function curve, thereby reliably estimating the global delay.) The window parameters are shown in Equation (1). (1) In the formula, Sampling rate, a and b are the start and end indices of the window used for cross-correlation calculation, M is the number of valid samples in the window, N is the total length of the signal, and L is the window length; Step 1.3.3, construct using equation (2) and Equal length; (2) In the formula, The first-order difference signal (approximate derivative) of the PPG signal is used for cross-correlation calculation with the ECG signal; Step 1.3.4: Using equations (3) and (4) in the window obtained in step 1.3.2 Top and Perform z-score standardization; (3) (4) In the formula, , These are the sample means of the ECG signal and the PPG differential signal within the window, respectively. Let be the sample standard deviations of the ECG signal and the PPG differential signal within the window, respectively. The ECG and PPG differential signals are Z-Score standardized to have a mean of 0 and a standard deviation of 1. Step 1.3.5: Let the discrete time delay index be... Search scope To ensure that the summation does not exceed the bounds, define... The aligned sequence is obtained through equations (5) and (6). Perform discrete cross-correlation and lag estimation; (5) (6) In the formula, This is a time delay index, representing the relative sample point offset between signals. , For a given delay Effective summation range indexes when calculating cross-correlation to prevent array out-of-bounds errors. Let be the cross-correlation value at time delay l. To make the cross-correlation value To achieve the maximum optimal latency index. For optimal latency; Step 1.3.6: Arrange the R peak as a whole according to formula (7). Translate and crop to the legal area; (7) In the formula, To be based on global latency The index of the i-th R-peak after preliminary translation correction; Step 1.3.7: Locally relock the R peak on the ECG wave using equation (8); (8).

[0020] In the formula, To the position of the R peak in the initial correction The radius of the window for performing a fine-grained local search in the vicinity. For the start and end indices of the local fine-grained search window, For the index of the i-th R-peak after final precise location, the window [ Internal ECG filter signal The maximum value point.

[0021] Step 2: The pulse transit time (PTT) determined by the R peak and PPG foot point of ECG is used as the core temporal feature; morphological features are extracted from the rising segment of PPG and the reflected wave, and intensity features (PIR) are calculated; at the same time, rhythmic features such as HR are calculated based on the RR interval, and a set of physiological parameter features is constructed by combining the subject's age and weight. Specifically, it includes: Step 2.1: Relocate the R peak of ECG using the method described in Step 1.3. Within a certain time window after the R peak, perform local upsampling on the rising segment of PPG to obtain the upsampled sequence. ;exist Find the point of maximum slope Its slope is denoted as ; A baseline of "valley value priority, anomaly backoff" is adopted: the amplitude of the previous valley value within the window is taken first, and if there is no valley value or an anomaly, the median of the preset upsampling start time is used; then let The original sampling index of the left endpoint of the local window is used as the basis for obtaining the PPG footpoint using equation (9). ,in The original sampling index of the R peak is used, and the PTT is finally calculated by equation (10); (9) (10) In the formula, This is the start time index for the local search window used to locate the foot point. This is the factor for local upsampling; in the actual code, it's 8. Let y(t) be the index of the point with the maximum slope in the upsampled signal y(t). This represents the maximum slope value of the rising segment of PPG. To determine the baseline amplitude of the foot point, To calculate the PPG foot time, The time index is the R peak, and PTT is the pulse wave transit time, which is the time difference from the R peak of the ECG to the PPG foot point. Step 2.2: Calculate the time from the PPG wave base to the main peak using equation (11). ,in This is the index of the PPG contraction peak in the original sampling coordinate system; (11) In the formula, For the time index of the main PPG peak (P1 peak), The rise time of the PPG wave; Step 2.3: Calculate the maximum upward slope of PPG using equation (12). That is, the maximum first derivative value of the rising segment of PPG, where PPG after bandpass and smoothing; (12) In the formula, This is the maximum first derivative of the rising segment of the PPG, multiplied by the sampling rate. To correct for amplitude scaling caused by the difference; Step 2.4: Calculate the ratio PIR of the AC component to the DC component of the pulse wave using equation (13), which is the ratio of the AC component to the DC component of the pulse wave; where Peak amplitude, It is the range of its previous trough value; For pulse wave alternating component, , is the DC component of the pulse wave; (13) In the formula, This represents the amplitude of the main peak of PPG. The amplitude value of the preceding trough paired with the main peak. The alternating component of the pulse wave. PIR represents the DC component of the pulse wave, and PIR represents the pulse wave intensity ratio. Step 2.5: Calculate the heart rate using formula (14), which is the time interval between two consecutive R waves; (14) In the formula, The i-th RR interval is the time interval between two consecutive R peaks. Heart rate; Step 2.6: Detect the p2 peak and retrieve the original sampling index. ,time In the fragment Above, let the smoothed amplitude be denoted as . ,in To be the main peak, The valley before the main peak. The next valley after the main peak; in the descending segment after the PPG main peak, the global minimum point is first located using equation (15). Then, all local maxima are found behind the notch using equation (16). As candidate P2 peaks, the candidate set is: Then, index all candidate indices. Convert to time And as in Equation (17), the P2 peak is screened by curvature verification, and finally the optimal candidate point is located by parabolic interpolation using Equation (18). (15) (16) (17) (18) In the formula, These are sampling indices representing the main peak, the valley before the main peak, and the valley after the main peak, respectively. This is the index of the notch location found in the PPG descent segment between the main peak and the next trough, i.e., the global minimum point of the signal in that segment. The set of all local maxima found after the notch is taken as the candidate P2 peak. The second derivative of the PPG signal. , This represents the time corresponding to the notch and the candidate peak of P2. , This represents the original index of the optimal P2 candidate peak and the signal amplitude at that point. This is the sub-pixel offset obtained through parabolic interpolation. The index of the P2 peak after subpixel-level precise positioning; Step 2.7: Calculate the reflectance index RI using equation (19), where It is the main peak of this auction. Index of the previous valley value of the pair Conversion time; (19) In the formula, The reflection index reflects the intensity of the reflected wave of the waveform. The amplitude of the main peak (P1), The amplitude of the reflection peak (P2) The time of the previous trough value paired with the main peak of this auction; Step 2.8: Calculate the systolic-diastolic period using equation (20). ; (20) In the formula, The systolic-diastolic phase; Step 2.9: Combine the obtained feature set with the corresponding subject's age and weight to jointly construct a physiological parameter feature set.

[0022] Step 3: Select a deep learning framework, build a Transformer hybrid model, and design a loss function; Specifically, it includes: Step 3.1: Select PyTorch as the deep learning framework and construct a fusion network. The network consists of a feature extraction module, the "Multi-Source Temporal Encoder Module" (MSEM), and a fusion-decoding module, the "Dynamic Conditional Fusion-Decoder Layer" (DCFD). The MSEM consists of an encoder and a middle layer. The encoder uses a 6-layer Transformer Encoder (128 hidden dimensions, 8 attention heads) to integrate ECG and PPG multi-channel temporal information and enhance long-range dependencies and noise resistance. The middle layer contains two types of sub-blocks: a branch generator and a Feature-wise Linear Modulation (FiLM) parameter generator. The branch generator incorporates three types of lightweight quantum heads: the Main head (a linear regression head after the encoder, with small parameters and stable convergence), the KAN head (Kolmogorov–Arnold network, with analytical approximation capabilities, suitable for nonlinear explicit features), and the eXtreme Gradient Boosting (XGB) head (gradient boosting tree ensemble, adept at tabular and statistical features, robust to small samples and outliers), and extends the residual correction head for error compensation and drift correction. The FiLM parameter generator is a two-layer perceptron structure, outputting the scale and offset coefficients per token for conditional recalibration, achieving individualization and domain adaptation. The Dynamic Conditional Fusion-Decoder Layer (DCFD) consists of a FiLM conditional modulation unit, a branch attention projection unit, a multi-head differential attention (MHDA) unit, a Transformer decoder, a token weighted aggregation unit, and a baseline-gated residual merging unit. The FiLM unit provides lightweight affine modulation to suppress inter-domain differences, while the projection unit unifies the representations of each branch into a comparable subspace. MHDA uses differential enhancement of dual-path similarity to highlight branches and reduce redundancy. The Transformer Decoder uses cross-attention to align branch representations with the main temporal context to improve temporal consistency. Token-weighted aggregation uses learnable weights to achieve interpretable fusion. Baseline gating uses the XGB output as a stable baseline, combined with sigmoid gating and a training-phase ramp strategy to achieve residual superposition and error correction, ultimately forming a dual-channel output of SBP and DBP. Step 3.2, set the input window as... Their window mean / variance are respectively Construct two types of flat features: the original expanded ⊕ statistics for input to XGB and residual tree. ) and the standardized expansion ⊕statistic for FiLM and KAN inputs ( ),in ; Step 3.3, will After embedding and position encoding, the input to the Transformer encoder yields the temporal implicit representation. Then, the intermediate branch generator produces three tokens in parallel. The Main branch takes... Terminal correlation step, via linear have to The KAN branch passes through get The XGB branch is passed through get The residual branch is obtained through the residual tree in the k-th round during the residual training process. Meanwhile, the FiLM parameter generator is based on Output the scale and offset parameters that match the current token count S. ; In the formula, For batch size, The time step length of the input sequence. The number of channels for the input feature. For the input tensor of the model, The mean and standard deviation of the input window X; Step 3.4: Dynamically pre-scale the weights of the three basic branches Main / KAN / XGB using equation (21), and reset the weights of the remaining residual tokens to 1. (twenty one) In the formula, The learnable scalar weights represent the importance of the three branches: Main, KAN, and XGB. To The branch weights after Softmax normalization of the three weights. The original token generated for the i-th branch. The token is pre-scaled with weights. Step 3.5: Conditional modulation is performed using FiLM. An affine transformation is applied to each token using equation (22) to achieve token-by-token feature reweighting; where, This represents the three basic branches and the original residual token after prescaling according to equation (21); (twenty two) In the formula, The original scale and offset parameters generated by MLP, , These are the scale and offset parameters that ultimately apply to the token after transformation. The token is FiLM modulated; Step 3.6: Each branch enters the multi-head differential attention module, first linearly increasing the dimensionality. Stacked Then, the similarity between the two paths is calculated at each attention head using equation (23), and differential reinforcement is performed. (twenty three) In the formula, V represents the Query, Key, and Value matrices, respectively, for two different attention paths. , The attention weight matrix calculated for the two paths. This is a weighting coefficient used to control the influence of the second attention path A2 in the difference process. This is the attention weight matrix after differential enhancement; Step 3.7: Each branch enters the Decoder and is associated across time-series decoders using equation (24); (twenty four) Step 3.8: Weight and aggregate each Token using equation (25) and map it to a 2D dimension, then obtain the result. Copy and align to the predicted step size get ; (25) Step 3.9: Integrate the residual round summation using equations (26) and (27) and merge it with the baseline gated residual, where For the first The learnable weights of the wheel residuals, For learnable gating, This refers to the warm-up coefficient during training. (26) (27) In the formula, The learnable ensemble weights for the residual tokens in the k-th round. For the residual token generated by the residual tree in the kth round, The signal is a learnable gated signal, constrained to the range (0,1) by the Sigmoid function. This is a warm-up coefficient during training, typically increasing linearly from 0 to 1, used to stabilize the initial stages of training. This is the final output of the model; Step 3.10, record the first One sample in the channel The truth value of is Predicted as The error is denoted as Let U be the total number of elements after batch and time step expansion, and let the k-th channel be... Statistics are Total Normalized weight definition Each channel also learns a log-variance parameter. l c Used for uncertainty estimation; the total loss function is designed based on the above, as follows: The total loss function is shown in Equation (28), which includes Gaussian negative log-likelihood, focal regression loss, quantile loss, Huber loss, KAN regularization, log-variance regularization and residual consistency. (28) The Gaussian negative log-likelihood loss is shown in Equation (29) and is used to learn the mean and variance simultaneously. (29) The focus regression loss is shown in Equation (30), emphasizing large errors and difficult samples; (30) The quantile loss, as shown in Equation (31), directly constrains the upper and lower quantiles, allowing the model to learn to predict "blood pressure range" rather than a single value, and is insensitive to abnormal systolic or diastolic readings, thus improving interval prediction and robustness. (31) The Huber loss is shown in Equation (32), which is quadratic for small errors and linear for large errors, taking into account both robustness and differentiability. (32) KAN regularization and logarithmic variance are shown in equations (33) and (34), respectively, which suppress overfitting and stabilize variance; (33) (34) The residual consistency is shown in Equation (35), which makes the residuals close to the true residuals when there is a baseline or residual branch; (35).

[0023] In formulas 28-35, For sample index, For output channel index, This represents the true blood pressure value of the u-th sample in the c-th channel. Let u be the predicted blood pressure value of the u-th sample in the c-th channel. For prediction error, The total number of elements in a batch, including samples and time periods. The log-variance parameter learned for the c-th channel is used for uncertainty estimation. Let be the L1 loss statistic for the k-th channel.

[0024] The normalized weight for the k-th channel. , , , , , , The weight hyperparameters for various loss terms, The target quantile for quantile loss is... The threshold parameter for Huber loss determines the switching point between secondary and primary losses.

[0025] Step 4: Select the learning rate decay strategy of the optimization function, set hyperparameters including the number of iterations and batch size, input the feature set into the network for training and test the model; Specifically, it includes: Choose the Lookahead-packaged RAdam optimizer, select the 10-round warm-up (Linear Learning Ratescheduler, LinearLR) and cosine annealing learning rate scheduling strategy, set hyperparameters including the number of iterations and batch size, input the feature set into the network for training, save the model for testing, and evaluate the model performance using the BHS standard and AAMI standard.

[0026] verify: like Figures 1-8 As shown, the detection accuracy of this invention for systolic blood pressure (SBP) and diastolic blood pressure (DBP) is 0.44±3.23 mmHg and 0.97±5.37 mmHg, respectively; the percentages of samples with absolute errors ≤5 mmHg are 71.8% and 92.7%, respectively.

[0027] The advantages of this invention are: 1. Multi-source signal fusion: Simultaneously utilizes ECG, PPG and physiological explicit features, resulting in high information utilization and greater robustness to individual differences and scene transfer; 2. Accurate data acquisition and positioning: Algorithms such as ECG R peak relocation, PPG foot point sampling-level positioning, and cross-modal time delay alignment can significantly reduce time stamp errors. At the same time, it has anomaly self-checking and polarity correction processing, and strong stability in data preprocessing and feature extraction. 3. Accurate and stable blood pressure estimation: Compared with traditional machine learning methods, this method performs better in blood pressure prediction; the evaluation results meet the AAMI standard and BHS Class A standard, and are suitable for continuous blood pressure monitoring, with certain medical application value.

[0028] The above are merely preferred embodiments of the present invention and are not intended to limit the present invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A non-invasive continuous blood pressure prediction method based on multidimensional data-driven approach, characterized in that, include: Step 1: Read the raw ECG and PPG signals and perform preprocessing; wherein, the preprocessing includes filtering, polarity correction, time alignment and normalization preprocessing; Step 2: Calculate the timing features, amplitude features, and waveform morphology features based on the preprocessed ECG and PPG signals, and construct a feature set by combining physiological parameters; Step 3: Select a deep learning framework, build a Transformer hybrid model, and design a loss function; Step 4: Select the optimization function and learning rate decay strategy, set hyperparameters including the number of iterations and batch size, input the feature set into the network for training and test the model.

2. The non-invasive continuous blood pressure prediction method as described in claim 1, characterized in that, Step 1 specifically includes: Step 1.1: Perform second-order Butterworth high-pass filtering and zero-phase band-pass filtering on the original ECG and PPG signals to obtain signals with noise and baseline drift removed for relocation; Step 1.2: After filtering, ensure the peak and valley directions of the PPG signal are consistent through waveform polarity self-check; Step 1.3: After polarity correction, the global micro-delay between ECG and PPG is estimated using the cross-correlation method, and the position of the R peak is corrected.

3. The non-invasive continuous blood pressure prediction method as described in claim 2, characterized in that, Step 1.3 specifically includes: Step 1.3.1: Let the sequences of ECG and PPG signals be... and n=0,1,...,N-1; the sequence after bandpass filtering is and The existing R-peak index set is The goal is to estimate the lag of the entire channel. And obtain the refined R peak ; Step 1.3.2: Select a window near the midpoint of the record. The window parameters are shown in equation (1); （1） In the formula, Sampling rate, a and b are the start and end indices of the window used for cross-correlation calculation, M is the number of valid samples in the window, N is the total length of the signal, and L is the window length; Step 1.3.3, construct using equation (2) and Equal length; （2） In the formula, It is the first-order differential signal of the PPG signal, used for cross-correlation calculation with the ECG signal; Step 1.3.4: Using equations (3) and (4) in the window obtained in step 1.3.2 Top and Perform z-score standardization; （3） （4） In the formula, , These are the sample means of the ECG signal and the PPG differential signal within the window, respectively. Let be the sample standard deviations of the ECG signal and the PPG differential signal within the window, respectively. The ECG and PPG differential signals are Z-Score standardized to have a mean of 0 and a standard deviation of 1. Step 1.3.5: Let the discrete time delay index be... Search scope To ensure that the summation does not exceed the bounds, define... The aligned sequence is obtained through equations (5) and (6). Perform discrete cross-correlation and lag estimation; （5） （6） In the formula, This is a time delay index, representing the relative sample point offset between signals. , For a given delay Effective summation range indexes when calculating cross-correlation to prevent array out-of-bounds errors. Let be the cross-correlation value at time delay l. To make the cross-correlation value To achieve the maximum optimal latency index. For optimal latency; Step 1.3.6: Arrange the R peak as a whole according to formula (7). Translate and crop to the legal area; （7） In the formula, To be based on global latency The index of the i-th R-peak after preliminary translation correction; Step 1.3.7: Locally relock the R peak on the ECG wave using equation (8); （8） In the formula, To the position of the R peak in the initial correction The radius of the window for performing a fine-grained local search in the vicinity. For the start and end indices of the local fine-grained search window, For the index of the i-th R-peak after final precise location, the window [ Internal ECG filter signal The maximum value point.

4. The non-invasive continuous blood pressure prediction method as described in claim 3, characterized in that, Step 2 specifically includes: Step 2.1: Relocate the R peak of ECG using the method described in Step 1.

3. Within a certain time window after the R peak, perform local upsampling on the rising segment of PPG to obtain the upsampled sequence. ;exist Find the point of maximum slope Its slope is denoted as ; A baseline of "valley value priority, anomaly backoff" is adopted: the amplitude of the previous valley value within the window is taken first, and the median of the preset upsampling start time is used when there is no valley value or an anomaly; then let The original sampling index of the left endpoint of the local window is used as the basis for obtaining the PPG footpoint using equation (9). ,in The original sampling index of the R peak is used, and the PTT is finally calculated by equation (10); （9） （10） In the formula, This is the start time index for the local search window used to locate the foot point. This is the factor for local upsampling; in the actual code, it's 8. Let y(t) be the index of the point with the maximum slope in the upsampled signal y(t). This represents the maximum slope value of the rising segment of PPG. To determine the baseline amplitude of the foot point, To calculate the PPG foot time, The time index is the R peak, and PTT is the pulse wave transit time, which is the time difference from the R peak of the ECG to the PPG foot point. Step 2.2: Calculate the time from the PPG wave base to the main peak using equation (11). ,in This is the index of the PPG contraction peak in the original sampling coordinate system; （11） In the formula, For the time index of the main PPG peak (P1 peak), The rise time of the PPG wave; Step 2.3: Calculate the maximum upward slope of PPG using equation (12). That is, the maximum first derivative value of the rising segment of PPG, where PPG after bandpass and smoothing; （12） In the formula, This is the maximum first derivative of the rising segment of the PPG, multiplied by the sampling rate. To correct for amplitude scaling caused by the difference; Step 2.4: Calculate the ratio PIR of the AC component to the DC component of the pulse wave using equation (13), which is the ratio of the AC component to the DC component of the pulse wave; where Peak amplitude, It is the range of its previous trough value; For pulse wave alternating component, , is the DC component of the pulse wave; （13） In the formula, This represents the amplitude of the main peak of PPG. The amplitude value of the preceding trough paired with the main peak. The alternating component of the pulse wave. PIR represents the DC component of the pulse wave, and PIR represents the pulse wave intensity ratio. Step 2.5: Calculate the heart rate using formula (14), which is the time interval between two consecutive R waves; （14） In the formula, The i-th RR interval is the time interval between two consecutive R peaks. Heart rate; Step 2.6: Detect the p2 peak and retrieve the original sampling index. ,time In the fragment Above, let the smoothed amplitude be denoted as . ,in To be the main peak, The valley before the main peak. The next valley after the main peak; in the descending segment after the PPG main peak, the global minimum point is first located using equation (15). Then, all local maxima are found behind the notch using equation (16). As candidate P2 peaks, the candidate set is: Then, index all candidate indices. Convert to time And as in Equation (17), the P2 peak is screened by curvature verification, and finally the optimal candidate point is located by parabolic interpolation using Equation (18). （15） （16） （17） （18） In the formula, These are sampling indices representing the main peak, the valley before the main peak, and the valley after the main peak, respectively. This is the index of the notch location found in the PPG descent segment between the main peak and the next trough, i.e., the global minimum point of the signal in that segment. The set of all local maxima found after the notch is taken as the candidate P2 peak. The second derivative of the PPG signal. , This represents the time corresponding to the notch and the candidate peak of P2. , This represents the original index of the optimal P2 candidate peak and the signal amplitude at that point. This is the sub-pixel offset obtained through parabolic interpolation. The index of the P2 peak after subpixel-level precise positioning; Step 2.7: Calculate the reflectance index RI using equation (19), where It is the main peak of this auction. Index of the previous valley value of the pair Conversion time; （19） In the formula, The reflection index reflects the intensity of the reflected wave of the waveform. The amplitude of the main peak (P1), The amplitude of the reflection peak (P2) The time of the previous trough value paired with the main peak of this auction; Step 2.8: Calculate the systolic-diastolic period using equation (20). ; （20） In the formula, The systolic-diastolic phase; Step 2.9: Combine the obtained feature set with the corresponding subject's age and weight to jointly construct a physiological parameter feature set.

5. The non-invasive continuous blood pressure prediction method as described in claim 4, characterized in that, Step 3 specifically includes: Step 3.1: Select PyTorch as the deep learning framework; Step 3.2, set the input window as... Their window mean / variance are respectively Construct two types of flat features: the original expanded ⊕ statistics for input to XGB and residual tree. ) and the standardized expansion ⊕statistic for FiLM and KAN inputs ( ),in ; Step 3.3, will After embedding and position encoding, the input to the Transformer encoder yields the temporal implicit representation. Then, the intermediate branch generator generates three tokens in parallel; the Main branch takes... Terminal correlation step, via linear have to The KAN branch passes through get The XGB branch is passed through get The residual branch is obtained through the residual tree in the k-th round during the residual training process. Meanwhile, the FiLM parameter generator is based on Output the scale and offset parameters that match the current token count S. ; In the formula, For batch size, The time step length of the input sequence. The number of channels for the input feature. For the input tensor of the model, The mean and standard deviation of the input window X; Step 3.4: Dynamically pre-scale the weights of the three basic branches Main / KAN / XGB using equation (21), and reset the weights of the remaining residual tokens to 1. （21） In the formula, The learnable scalar weights represent the importance of the three branches: Main, KAN, and XGB. To The branch weights after Softmax normalization of the three weights. The original token generated for the i-th branch. The token is pre-scaled with weights. Step 3.5: Conditional modulation is performed using FiLM. An affine transformation is applied to each token using equation (22) to achieve token-by-token feature reweighting; where, This represents the three basic branches and the original residual token after prescaling according to equation (21); （22） In the formula, The original scale and offset parameters generated by MLP, , These are the scale and offset parameters that ultimately apply to the token after transformation. The token is FiLM modulated; Step 3.6: Each branch enters the multi-head differential attention module, first linearly increasing the dimensionality. Stacked Then, the similarity between the two paths is calculated at each attention head using equation (23), and differential reinforcement is performed. （23） In the formula, V represents the Query, Key, and Value matrices, respectively, for two different attention paths. , The attention weight matrix calculated for the two paths. This is a weighting coefficient used to control the influence of the second attention path A2 in the difference process. This is the attention weight matrix after differential enhancement; Step 3.7: Each branch enters the Decoder and is associated across time-series decoders using equation (24); （24） Step 3.8: Weight and aggregate each Token using equation (25) and map it to a 2D dimension, then obtain the result. Copy and align to the predicted step size get ; （25） Step 3.9: Integrate the residual round summation using equations (26) and (27) and merge it with the baseline gated residual, where For the first The learnable weights of the wheel residuals, For learnable gating, This refers to the warm-up coefficient during training. （26） （27） In the formula, The learnable ensemble weights for the residual tokens in the k-th round. For the residual token generated by the residual tree in the kth round, The signal is a learnable gated signal, constrained to the range (0,1) by the Sigmoid function. This is a warm-up coefficient during training, typically increasing linearly from 0 to 1, used to stabilize the initial stages of training. This is the final output of the model; Step 3.10, record the first One sample in the channel The truth value of is Predicted as The error is denoted as Let U be the total number of elements after batch and time step expansion, and let the k-th channel be... Statistics are Total Normalized weight definition Each channel also learns a log-variance parameter. l c Used for uncertainty estimation; the total loss function is designed based on the above, as follows: The total loss function is shown in Equation (28), which includes Gaussian negative log-likelihood, focal regression loss, quantile loss, Huber loss, KAN regularization, log-variance regularization and residual consistency. （28） The Gaussian negative log-likelihood loss is shown in Equation (29) and is used to learn the mean and variance simultaneously. （29） The focus regression loss is shown in Equation (30), emphasizing large errors and difficult samples; （30） The quantile loss is shown in equation (31), which directly constrains the upper and lower quantiles; （31） The Huber loss is shown in Equation (32), which is quadratic for small errors and linear for large errors, taking into account both robustness and differentiability. （32） KAN regularization and logarithmic variance are shown in equations (33) and (34), respectively, which suppress overfitting and stabilize variance; （33） （34） The residual consistency is shown in Equation (35), which makes the residuals close to the true residuals when there is a baseline or residual branch; （35） In formulas 28-35, For sample index, For output channel index, This represents the true blood pressure value of the u-th sample in the c-th channel. Let u be the predicted blood pressure value of the u-th sample in the c-th channel. For prediction error, The total number of elements in a batch, including samples and time periods. The log-variance parameter learned for the c-th channel is used for uncertainty estimation. Let L1 loss statistic be the L1 loss statistic for the k-th channel; The normalized weight for the k-th channel. , , , , , , The weight hyperparameters for various loss terms, The target quantile for quantile loss is... The threshold parameter for Huber loss determines the switching point between secondary and primary losses.

6. The non-invasive continuous blood pressure prediction method as described in claim 5, characterized in that, Step 4 specifically includes: Choose the Lookahead-packaged RAdam optimizer, select a learning rate scheduling strategy of 10 rounds of warm-up and cosine annealing, set hyperparameters including the number of iterations and batch size, input the feature set into the network for training, save the model for testing, and evaluate the model performance using the BHS standard and AAMI standard.

Images