Biofeedback-based Diver Status Monitoring Method

Abstract

This application belongs to the technical field of diving state monitoring, and relates to a diver state monitoring method based on biofeedback. By aligning the physiological signs signal and environmental data through timestamps, and using acceleration data to extract motion features, a motion noise model is constructed to remove signal artifacts in real time; in terms of EEG signal purification, db4 wavelet basis is used for wavelet packet decomposition, features of water pressure-related frequency bands are extracted and low-frequency drift is removed. At the same time, the heart rate signal is analyzed by nonlinear dynamics, Lyapunov exponents are extracted to evaluate HRV changes, and wavelet transform is used to analyze the coupling between respiration and EEG #imgabs0# waves to further enhance the accuracy of risk assessment; in terms of environmental compensation, through the relationship between depth and SpO2 and the compensation of water temperature for respiratory rate, combined with the SLTM network to extract temporal features, and a deep causal convolution layer is used for risk trend prediction, the diver risk level and prediction curve are output in real time, and finally the comprehensive risk index is calculated to comprehensively evaluate the diver's health status and potential risks and improve safety.

Description

Technical Field

[0001] This application belongs to the technical field of diving state monitoring. More specifically, it relates to a method for monitoring the state of divers based on biofeedback. Background Art

[0002] With the popularization of underwater operations and diving activities, the health and safety issues of divers have received increasing attention. During diving, the physiological state of divers (such as heart rate, respiration, electroencephalogram activity, etc.) is affected by various factors, such as environmental parameters like water depth, pressure, temperature, as well as the physical load and psychological state of divers. To ensure the safety of divers, real-time monitoring of their physiological state and timely risk assessment have become a key research area.

[0003] Currently, biofeedback technology has been widely used in the medical and sports fields. Especially in health monitoring, some mature technologies have been formed. Through the real-time acquisition and analysis of biological signals (such as electrocardiogram, electroencephalogram, respiration signal, etc.), the physical state of divers can be accurately evaluated, and the potential risks of divers can be predicted. However, in a complex underwater environment, how to effectively remove noise interference, improve the accuracy of signals, and on this basis, conduct efficient and accurate risk assessment is still a technical problem.

[0004] For example, divers are constantly moving in the underwater environment, which causes biological signals (especially EEG and heart rate signals) to be severely interfered by motion artifacts. Existing noise suppression methods (such as traditional Kalman filtering and recursive least squares method) can remove some artifacts, but still cannot effectively remove motion artifacts in a complex environment, resulting in unstable signal quality and affecting subsequent risk assessment. Summary of the Invention

[0005] The present invention provides a method for monitoring the state of divers based on biofeedback, which solves the technical problem that the prior art cannot effectively remove noise interference in a complex environment.

[0006] The method for monitoring the state of divers based on biofeedback includes the following steps:

[0007] Data acquisition: Collect physiological signal data and environmental data, align the collected data with time stamps, and obtain multi-modal signal stream data aligned on the time axis.

[0008] Motion noise suppression: Extract motion features from the acceleration data in the physiological signal data to construct a motion noise model, and update the filter weights of the motion noise model in real time through the recursive least squares algorithm. Remove artifacts through the motion noise model to obtain blood flow signals after removing artifacts.

[0009] EEG Signal Enhancement: For the EEG signal data in the vital sign signal data, perform 6-layer wavelet packet decomposition using the db4 wavelet basis, extract features in different frequency bands from the EEG signal, identify the energy distribution in the frequency bands related to water pressure, remove the low-frequency baseline drift through thresholding, and obtain the purified EEG signal;

[0010] Time-domain Feature Generation: Perform nonlinear dynamics analysis on the denoised heart rate signal, evaluate the complexity of the signal, construct a phase space reconstruction diagram, extract the largest Lyapunov exponent to reflect the dynamic changes of the heart rate signal, and obtain the nonlinear feature vector of HRV;

[0011] Frequency-domain Feature Generation: Based on the purified EEG signal and respiratory signal, use wavelet transform to analyze the coherence between the respiratory signal and the EEG waves, and quantify the coupling degree between the respiratory signal and the EEG waves through the phase synchronization index; obtain the coupling feature between respiration and EEG;

[0012] Environmental Parameter Compensation: Based on the depth, establish a relationship to perform depth compensation on SpO2, and compensate the respiratory rate according to the water temperature;

[0013] Short-term Risk Inference: Based on the nonlinear feature vector of HRV and the coupling feature between respiration and EEG, perform temporal feature extraction through the LSTM network structure, and then map it to the short-term risk probability based on the fully connected layer to output the real-time risk level;

[0014] Medium- and Long-term Trend Prediction: Construct a deep causal convolution layer to predict the future risk evolution, and focus on the key time periods through the attention mechanism to output the risk trend prediction curve;

[0015] Comprehensive Risk Index Output: Based on the real-time risk level and the risk trend prediction curve, perform weighted fusion using the compensated environmental parameters to obtain the comprehensive risk index.

[0016] In the present invention, in the data acquisition stage, by aligning the timestamps of the vital sign signals and environmental data, the consistency of the multi-modal signals is ensured, providing accurate data support for subsequent processing; for the signal artifacts caused by the movement of divers, use the acceleration data to extract the motion features, construct a motion noise model, and update the filter weights in real time through the recursive least squares algorithm to remove the artifacts and ensure the accuracy of the blood flow signals.

[0017] In terms of EEG signal purification, the db4 wavelet basis is used for six-layer wavelet packet decomposition to extract the features in the frequency band related to water pressure, and the low-frequency baseline drift is removed through thresholding, enhancing the quality of the EEG signal. The heart rate signal is analyzed by nonlinear dynamics to extract the maximum Lyapunov exponent to evaluate the dynamic changes of HRV, providing key features for short-term risk inference; the coupling features between the respiratory signal and the EEG theta wave are analyzed by wavelet transform to further improve the accuracy of risk assessment.

[0018] In terms of environmental compensation, based on the relationship between depth and SpO2 and the compensation of water temperature for respiratory rate, the accuracy of physiological data is ensured. The time-series features are extracted through the LSTM network structure, and combined with the deep causal convolution layer for medium- and long-term risk trend prediction, and the risk level and risk prediction curve of the diver are output in real time; finally, by integrating the real-time risk level and trend prediction curve, and combining the compensation of environmental parameters, the comprehensive risk index is obtained, comprehensively and accurately evaluating the health status and potential risks of the diver, effectively improving the safety of the diver.

[0019] Preferably, the motion noise model is as follows: ; where: represents the motion noise; represents the time-varying weight parameter; represents the filtered acceleration signal;

[0020] The time-varying weight parameter is updated in real time through the recursive least squares algorithm. After obtaining the motion noise, the motion noise is subtracted from the original blood flow signal to obtain the blood flow signal after artifact removal.

[0021] Preferably, the EEG signal enhancement includes the following steps:

[0022] Wavelet packet decomposition: Six-layer decomposition is performed using the db4 wavelet basis, and the signal is decomposed into 64 sub-bands, and each band represents the energy of the signal in a specific frequency range; based on the decomposed sub-bands, the water pressure interference frequency band and the low-frequency drift frequency band are located;

[0023] Water pressure interference suppression: Energy analysis is performed in the water pressure interference frequency band, and the energy ratio of each frequency band node is calculated. If the energy ratio of the water pressure interference frequency band is higher than a predetermined threshold, the corresponding wavelet packet coefficients are attenuated through soft threshold technology, and the effective components of the EEG signal are retained;

[0024] Low-frequency baseline drift removal: For the baseline drift below 0.5 Hz, by performing additional decomposition on the wavelet packet tree, the components in the range of 0 - 1 Hz are removed from the signal;

[0025] Wavelet packet reconstruction: The processed frequency point nodes are reconstructed, and the purified EEG time-domain signal is obtained by synthesizing the wavelet packet coefficients layer by layer;

[0026] Time-frequency feature matrix generation: Extract the energy within the time window of each node from the processed signal and perform normalization to obtain the time-frequency feature matrix.

[0027] Preferably, the specific steps for generating the time-domain features are as follows:

[0028] Sample entropy calculation: Based on the denoised heart rate signal, transform the heart rate signal into a high-dimensional vector through phase space reconstruction. During the reconstruction process, the embedding dimension and time delay are adopted, where the time delay is determined by the autocorrelation function method; then calculate the matching degree between high-dimensional vectors to obtain the sample entropy;

[0029] Maximum Lyapunov exponent extraction: Based on the denoised heart rate signal, construct a phase space trajectory to obtain the trajectory, then find the nearest point for each point, calculate the divergence speed between adjacent points, and further calculate the maximum Lyapunov exponent;

[0030] Correlation dimension calculation: Based on the denoised heart rate signal, calculate the dimensional characteristics of the phase space through the correlation function to obtain the similarity between two points in the phase space, and then fit the relationship between the two points in the double logarithmic coordinate system to obtain the correlation dimension;

[0031] Quadratic sample entropy calculation: Based on the denoised heart rate signal, obtain the quadratic sample entropy value through similarity matching;

[0032] Detrended fluctuation analysis: Based on the denoised heart rate signal, calculate the cumulative deviation sequence to obtain the cumulative deviation of each time point relative to the mean value. Divide the sequence into multiple small intervals, calculate the local volatility within each interval, and finally obtain the DFA scaling exponent through double logarithmic fitting;

[0033] The results obtained based on sample entropy calculation, maximum Lyapunov exponent extraction, correlation dimension calculation, quadratic sample entropy calculation, and detrended fluctuation analysis form a five-dimensional non-linear feature vector.

[0034] Preferably, the generation of frequency-domain features includes the following steps:

[0035] EEG EEG wave component extraction: Extract the EEG wave components from the EEG time-frequency matrix based on neurophysiological characteristics wave components;

[0036] Continuous wavelet transform of respiratory signal: Use the complex Morlet wavelet to perform continuous wavelet transform on the respiratory signal to obtain the wavelet coefficients of the respiratory signal in the time-frequency domain, that is, the time-frequency representation of the respiratory signal;

[0037] Cross-wavelet spectrum calculation: Calculate the time-frequency representation of the respiratory signal and EEG The cross-wavelet spectrum of the wave is obtained by multiplying the time-frequency coefficients of the respiration signal by the complex conjugate wavelet coefficients of the EEG wave, and the normalized wavelet coherence is calculated based on the obtained cross-wavelet spectrum; The cross-wavelet spectrum is obtained by multiplying the time-frequency coefficients of the respiration signal by the complex conjugate wavelet coefficients of the EEG wave, and the normalized wavelet coherence is calculated based on the obtained cross-wavelet spectrum;

[0038] Phase synchronization index extraction: Extract the phase information from the cross-wavelet spectrum, calculate the instantaneous phase difference between the respiration signal and the EEG wave at each time-frequency point, calculate the calculated instantaneous phase difference, and calculate the global phase synchronization index within the frequency band to obtain the phase synchronization index; The global phase synchronization index within the frequency band is obtained to obtain the phase synchronization index;

[0039] Feature vector generation: Calculate the time average of the coherence within the frequency band based on the normalized wavelet coherence to obtain the average coherence intensity; identify continuous time periods with coherence exceeding the threshold, and calculate the maximum duration; finally, calculate the variance of the phase synchronization index within the sliding window, and based on this, obtain the 3D coupling characteristics of respiration and EEG.

[0040] Preferably, the environmental parameter compensation includes the following steps:

[0041] Nonlinear compensation of blood oxygen saturation by depth: A quadratic regression model is used to correct the blood oxygen saturation, where the linear term and the quadratic term of the quadratic regression model respectively compensate for the deviation and over-compensation problems caused by depth;

[0042] Linear compensation of respiration frequency by water temperature: A respiration frequency compensation model is established according to the Q10 effect, where linear compensation is achieved through the temperature sensitivity coefficient;

[0043] Multi-parameter coupling effect processing: Adjust the respiration frequency by considering the cross-action of depth and water temperature, and then use a Kalman filter to fuse the data of the pressure gauge and the sonar, and achieve accurate depth estimation by dynamically adjusting the Kalman gain.

[0044] Preferably, the comprehensive risk index output includes the following steps:

[0045] Membership function definition: Quantify the impact of environmental parameters on the comprehensive risk through the membership function, where the depth weight gradually increases with the increase of depth, and the body temperature weight is dynamically adjusted according to the degree of deviation of the body temperature from the reference value;

[0046] Weighted fusion: Perform weighted fusion according to the results of the real-time risk and the prediction trend, and output the final comprehensive risk index.

[0047] The beneficial effects of the present invention include:

[0048] In the present invention, during the data acquisition stage, the consistency of multimodal signals is ensured by aligning the timestamps of vital signs signals and environmental data, providing accurate data support for subsequent processing; in order to address the signal artifacts caused by the diver's movement, the acceleration data is used to extract motion features, construct a motion noise model, and update the filter weights in real time through a recursive least squares algorithm to remove artifacts and ensure the accuracy of the blood flow signal.

[0049] In terms of EEG signal purification, the db4 wavelet base is used to perform 6-layer wavelet packet decomposition to extract the characteristics of the water pressure-related frequency band, and the low-frequency baseline drift is removed through thresholding processing to enhance the quality of the EEG signal. The heart rate signal is analyzed through nonlinear dynamics, and the maximum Lyapunov exponent is extracted to evaluate the dynamic changes of HRV, providing key features for short-term risk reasoning. The wavelet transform is used to analyze the respiratory signal and EEG The coupled characteristics of the waves further improve the accuracy of risk assessment.

[0050] In terms of environmental compensation, the accuracy of physiological data is ensured based on the relationship between depth and SpO2 and the compensation of water temperature for respiratory rate. Time series features are extracted through the LSTM network structure, and medium- and long-term risk trend predictions are performed in combination with the deep causal convolution layer, outputting the diver's risk level and risk prediction curve in real time. Finally, by integrating the real-time risk level and trend prediction curve, and combining it with compensation for environmental parameters, a comprehensive risk index is obtained, which comprehensively and accurately assesses the diver's health status and potential risks, effectively improving the diver's safety. BRIEF DESCRIPTION OF THE DRAWINGS

[0051] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the following briefly introduces the drawings required for use in the embodiments or descriptions of the prior art. Obviously, the drawings described below are only some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without any creative work.

[0052] Figure 1 This is a flowchart of the overall steps provided by an embodiment of the present invention.

[0053] Figure 2 This is a block diagram of the EEG signal enhancement steps provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0054] In order to make the technical problems, technical solutions and beneficial effects to be solved by this application more clearly understood, this application is further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain this application and are not intended to limit this application.

[0055] See alsoFigure 1 As shown in Figure 1 , a biofeedback-based diver status monitoring method includes the following steps:

[0056] Data acquisition: Collect physiological signal data and environmental data, align the collected data with timestamps to obtain multi-modal signal stream data aligned on the time axis;

[0057] Wherein the physiological signal data includes optical blood flow signal, piezoelectric respiration signal, EEG electroencephalogram signal, and EMG electromyogram signal; the environmental parameters include depth, temperature, and motion acceleration;

[0058] Timestamp alignment uses a hardware timestamp to accurately align asynchronous signals, with a timestamp accuracy of ±1 ms; for signals with different sampling rates, resampling is performed to make their sampling rates consistent, and interpolation methods such as linear interpolation or spline interpolation are used to ensure time axis alignment; finally, the multi-modal signal data stream aligned on the time axis is output.

[0059] Motion noise suppression: Extract motion features from the acceleration data in the physiological signal data to construct a motion noise model, and update the filter weights of the motion noise model in real time through the recursive least squares algorithm. Remove artifacts through the motion noise model to obtain the blood flow signal after removing artifacts;

[0060] As a possible implementation of this embodiment, assume the original optical blood flow signal contains the true blood flow signal and motion noise ; the triaxial accelerometer signal corresponds to the acceleration signals of the X, Y, and Z axes respectively;

[0061] First, apply zero-phase Butterworth band-pass filtering to the accelerometer signal to filter out signal components above or below the 2 - 10 Hz frequency band; obtain the filtered acceleration signal;

[0062] Express the motion noise as a weighted combination of triaxial acceleration signals, so the constructed motion noise model is as follows: ;

[0063] In the formula: represents the motion noise; represents the time-varying weight parameter; represents the filtered acceleration signal;

[0064] Update the time-varying weight parameter in real time through the recursive least squares algorithm. After obtaining the motion noise, subtract the motion noise from the original blood flow signal to obtain the blood flow signal after removing artifacts;

[0065] When the time-varying weight parameters are updated in real time through the recursive least squares algorithm, the estimation error of the motion noise is minimized; an acceleration vector is constructed based on the filtered acceleration signal, then a gain vector is calculated based on the constructed acceleration vector, and then an a priori error is calculated based on the original optical blood flow signal and the acceleration vector. Finally, the current weight vector is calculated based on the calculated a priori error, the gain vector, and the weight at the previous moment; based on this, the updated time-varying weight parameters are obtained; on the basis of determining the data used and the purpose to be achieved, how to update the time-varying weight parameters through the recursive least squares algorithm is a conventional technical means in the art. Therefore, in the present invention, the implementation principle of the recursive least squares algorithm will not be described in detail.

[0066] In this embodiment, through the extraction and analysis of acceleration data by the motion noise model, the artifacts caused by the movement of the diver can be effectively removed, the purity of the blood flow signal can be significantly improved, and the filter weights of the motion noise model are updated in real time by using the recursive least squares algorithm, so that the model can dynamically adapt to the change of the diver's motion state, and the denoising effect is further improved.

[0067] EEG signal enhancement: For the EEG signal data in the vital sign signal data, 6-layer wavelet packet decomposition is performed using the db4 wavelet basis, features in different frequency bands are extracted from the EEG signal, and the energy distribution of the water pressure-related frequency band is identified. The low-frequency baseline drift is removed through thresholding to obtain the purified EEG signal;

[0068] See Figure 2 As shown, as a possible implementation manner of this embodiment, the EEG signal enhancement includes the following steps:

[0069] Wavelet packet decomposition: 6-layer decomposition is performed using the db4 wavelet basis, and the signal is decomposed into 64 sub-bands, and each band represents the energy of the signal in a specific frequency range; the water pressure interference band and the low-frequency drift band are located based on the decomposed sub-bands;

[0070] Exemplary: The original EEG signal contains underwater electromagnetic interference and low-frequency baseline drift (<0.5 Hz), and the sampling rate , that is, 256 data points are sampled per second, and the Nyquist frequency ;

[0071] The decomposition level L = 6, generating frequency band nodes, and the node numbering rule for each layer is the nth node in the lth layer;

[0072] The frequency range of each node is as follows: ;

[0073] Where: denotes the low-frequency boundary of the nth node in the lth layer; n represents the node number; L represents the decomposition layer number; denotes the high-frequency boundary of the nth node in the lth layer;

[0074] Based on the above formula, we can obtain that the water pressure interference frequency band of 14 - 18 Hz corresponds to nodes (6,7), (6,8), (6,9) (i.e., 14 - 16 Hz, 16 - 18 Hz, 18 - 20 Hz); the low-frequency drift frequency band below 0.5 Hz corresponds to node (6,0) (i.e., 0 - 2 Hz);

[0075] Water pressure interference suppression: Conduct energy analysis in the water pressure interference frequency band, calculate the energy proportion of each frequency band node. If the energy proportion in the water pressure interference frequency band is higher than the predetermined threshold, attenuate the corresponding wavelet packet coefficients through the soft threshold technology, and retain the effective components of the EEG signal;

[0076] Exemplarily: Calculate the energy of each node for the target nodes (6,7), (6,8), (6,9):

[0077] ;

[0078] In the formula: denotes the energy of node n; denotes the wavelet packet coefficient of node n in the time domain; k represents the discrete points in the time domain;

[0079] Sum up the energies calculated for each node to obtain the total energy, and then calculate the energy proportion of each node in the total energy based on the total energy; if the energy proportion is greater than the threshold, it is determined as significant interference;

[0080] Perform soft threshold processing on the interference node coefficients: ;

[0081] ;

[0082] In the formula: denotes the adaptive threshold; N represents the signal length; denotes the noise standard deviation of node n, estimated by the minimum scale coefficient; denotes the sign function of; denotes attenuating the coefficient by the threshold. If the absolute value of the coefficient is less than then set it to 0 and retain the effective components;

[0083] Low-frequency baseline drift removal: For the baseline drift below 0.5 Hz, remove the components of 0 - 1 Hz from the signal by performing additional decomposition on the wavelet packet tree;

[0084] Exemplarily, the node (6,0) contains components of 0 - 2 Hz. It is necessary to further isolate the part less than 0.5 Hz. Perform an additional 1 - layer decomposition on the node (6,0) to obtain the sub - node (7,0) of 0 - 1 Hz and the sub - node (7,1) of 1 - 2 Hz; directly discard the coefficients of the sub - node (7,0), that is ; In this way, the components below 0.5 Hz are removed;

[0085] Wavelet packet reconstruction: Reconstruct each processed frequency - point node, and obtain the purified EEG time - domain signal by synthesizing wavelet packet coefficients layer by layer;

[0086] Exemplarily: Recursively synthesize wavelet packet coefficients from the 6th layer to the 0th layer: ;

[0087] In the formula: g and h are the low - pass and high - pass filters for wavelet reconstruction (db4 wavelet basis); and respectively represent the low - frequency and high - frequency coefficients of the sub - band; k, m represent discrete points in the time domain; finally, obtain the purified EEG time - domain signal through wavelet reconstruction;

[0088] Time - frequency feature matrix generation: Extract the energy within the time window of each node from the processed signal and perform normalization processing to obtain the time - frequency feature matrix.

[0089] In this embodiment, perform 6 - layer decomposition using the db4 wavelet basis to refine the EEG signal into multiple frequency bands, which can more accurately locate and suppress water - pressure interference and low - frequency drift, effectively remove the interference of water pressure on the signal. Especially when the water pressure changes greatly in the diving environment, the quality of the EEG signal can be significantly improved; and by attenuating the wavelet packet coefficients of the water - pressure interference frequency band and retaining the effective components of the EEG signal, it is beneficial to improve the signal - to - noise ratio and enhance the reliability and accuracy of the signal; for the baseline drift signal below 0.5 Hz, by performing additional decomposition on the wavelet packet tree, the low - frequency drift can be effectively removed; in this way, the low - frequency noise that may affect the signal accuracy and stability is removed, and the time - domain accuracy of the EEG signal is improved.

[0090] Time - domain feature generation: Use the denoised heart - rate signal for non - linear dynamics analysis, evaluate the complexity of the signal, construct a phase - space reconstruction diagram, extract the maximum Lyapunov exponent to reflect the dynamic changes of the heart - rate signal, and obtain the non - linear feature vector of HRV;

[0091] As a possible implementation manner of this embodiment, the specific steps of the time - domain feature generation are as follows:

[0092] Sample entropy calculation: Based on the denoised heart rate signal, the heart rate signal is converted into a high-dimensional vector through phase space reconstruction. The embedding dimension and time delay are used in the reconstruction process, and the time delay is determined by the autocorrelation function method; then, the matching degree between high-dimensional vectors is calculated to obtain the sample entropy.

[0093] Exemplary: Map the RR interval sequence (heart rate signal) to an m-dimensional phase space to construct a state vector: ;

[0094] where: m represents the embedding dimension; represents the time delay, which is determined by the autocorrelation function method, and the delay at the first zero crossing is selected;

[0095] Similarity matching: Define the distance threshold as , where represents the standard deviation of the RR interval sequence; Calculate the number of template matches that meet the conditions:

[0096] 1.A: The number of pairs of distances of m-dimensional vectors ( );

[0097] 2.B: The number of pairs of distances between (m + 1)-dimensional vectors, represents the distance between (m + 1)-dimensional state vectors;

[0098] where the sample entropy is defined as: ;

[0099] In the formula: A and B represent the number of matches in m - 1 dimensions and (m + 1) dimensions respectively. The smaller the sample entropy, the stronger the regularity of the sequence, and the larger the sample entropy, the higher the complexity of the sequence;

[0100] Maximum Lyapunov exponent extraction: Based on the denoised heart rate signal, a phase space trajectory is constructed to obtain the trajectory, and then for each point, its nearest neighbor point is found, and the divergence speed between the neighboring points is calculated, and then the maximum Lyapunov exponent is calculated;

[0101] Exemplary, use the reconstructed m-dimensional phase space vector to describe the phase space of the heart rate signal;

[0102] For each phase point find its nearest neighbor point , satisfying ; where represents the distance between two phase space points;

[0103] Track the exponential divergence of neighboring points over time: ;

[0104] In the formula: represents the discrete time step (k = 1, 2, …, K, the maximum step length K = 50), represents the number of effective tracking points; represents the divergence rate;

[0105] Fit the linear part of the divergence rate curve: ;

[0106] In the formula: represents the maximum Lyapunov exponent, representing the sensitivity to the initial conditions; represents the function of the divergence rate changing with time; t represents the time step; c represents the fitting constant;

[0107] Calculation of correlation dimension: Based on the denoised heart rate signal, calculate the dimensional characteristics of the phase space through the correlation function, obtain the similarity between two points in the phase space, and then fit the relationship between the two points in the double logarithmic coordinates to obtain the correlation dimension;

[0108] Exemplarily, calculate the correlation integral: ;

[0109] Among them, represents the correlation integral; represents the Heaviside step function, which is used to judge whether the distance between two points is less than the threshold r; and both represent the state vectors in the phase space; N represents the total number of data points;

[0110] By calculating whether the distance between two points in the phase space is less than the threshold r and accumulating the pairs that meet the conditions;

[0111] In the double logarithmic coordinates, select the linear region for fitting: ;

[0112] In the formula: represents the correlation dimension, which is the dimensional characteristic of the phase space; represents taking the natural logarithm of the distance r; represents the constant term;

[0113] Calculation of sample entropy of order two: Based on the denoised heart rate signal, obtain the sample entropy value of order two through similarity matching;

[0114] Calculation of sample entropy of order two: Based on the denoised heart rate signal, obtain the sample entropy value of order two through similarity matching;

[0115] Exemplarily, set the embedding dimension m = 3, and other parameters are the same as the steps of sample calculation; perform phase space reconstruction and similarity matching in the same way as sample entropy calculation to obtain the sample entropy of order two , where and represent the number of matches when m = 3 obtained during the calculation process;

[0116] Detrended Fluctuation Analysis: Calculate the cumulative deviation sequence based on the denoised heart rate signal, obtain the cumulative deviation of each time point relative to the mean, divide the sequence into multiple small intervals, calculate the local volatility within each interval, and finally obtain the DFA scaling exponent through double logarithmic fitting;

[0117] Exemplarily, it includes the following steps:

[0118] Accumulate the deviation sequence: ;

[0119] In the formula: represents the mean of the RR interval sequence; represents the i-th RR interval, indicating the time interval between two R waves in the vector of the heart rate signal; represents the mean of the RR interval sequence; N represents the total length of the RR interval sequence;

[0120] Interval fitting and fluctuation calculation: Divide the sequence into intervals of length s, and calculate the local trend fitting residual: ;

[0121] In the formula: represents the local trend fitting of the v-th interval; represents the local volatility measure; s represents the interval length; represents the total number of intervals divided in the data sequence; v represents the interval index; represents the value of the cumulative deviation sequence at the k-th time point; represents the weight of each data point within the interval;

[0122] Scaling exponent extraction: ;

[0123] In the formula: represents the DFA scaling exponent, and this exponent is extracted by fitting the interval to reflect the long-term correlation; represents the local volatility of the logarithm; represents the logarithm of the interval length s; represents the interval length range;

[0124] The results obtained based on sample entropy calculation, maximum Lyapunov exponent extraction, correlation dimension calculation, quadratic sample entropy calculation, and detrended fluctuation analysis form a five-dimensional non-linear feature vector.

[0125] In this embodiment, by calculating the sample entropy and the largest Lyapunov exponent, the complexity and chaotic characteristics of the heart rate signal can be effectively evaluated. These non-linear indicators can reflect the dynamic changes in the physiological state of divers, especially their responses under stress or extreme environments, providing more in-depth information for real-time health monitoring. By calculating the correlation dimension of the heart rate signal, the physiological state of divers can be further revealed. A high-dimensional correlation dimension usually indicates a higher complexity of the signal, and this complexity may be related to the physiological states of divers such as fatigue and stress. The DFA scaling exponent can reveal the long-term dependence and local volatility of the heart rate signal, helping to evaluate the physiological adaptability of divers in the underwater environment, especially in the case of deep diving and long-term diving; it helps to monitor physiological changes such as fatigue and stress of divers.

[0126] Frequency domain feature generation: Based on the purified EEG signal and respiratory signal, wavelet transform is used to analyze the coherence between the respiratory signal and the EEG waves, and the coupling degree between the respiratory signal and the EEG waves is quantified through the phase synchronization index; the coupling characteristics of respiration and EEG are obtained;

[0127] As a possible implementation manner of this embodiment, the frequency domain feature generation includes the following steps:

[0128] EEG wave extraction, extracting the frequency components of the waves (4 - 7 Hz) from the EEG time-frequency matrix. According to the neurophysiological characteristics, the frequency range of the

[0129] waves is 4 - 7 Hz. First, select the frequency components in this frequency band from the EEG time-frequency matrix:

[0130] where: represents the wave component at time point t; represents the wavelet coefficient of the EEG time-frequency matrix at time point t and frequency k; t represents the time point; represents selecting the wave frequency band from the time-frequency matrix of the EEG signal;

[0131] Perform continuous wavelet transform on the respiratory signal, and extract the time-frequency characteristics related to the same frequency band as the waves. First, perform continuous wavelet transform (CWT) on the respiratory signal using the complex Morlet wavelet. The Morlet wavelet basis function is defined as follows: ; ;

[0132] where: f represents the analysis frequency, select the same as Wave-aligned frequency range 4 - 7 Hz; The bandwidth factor of the wavelet, set to 6, ensuring the balance of time-frequency resolution; j represents the imaginary unit, in complex form; Represents the Morlet wavelet basis function; Represents the time variable, which is the time scale in wavelet transform; Represents the Gaussian envelope of the wavelet; Represents the high-frequency component of the wavelet;

[0133] Use Morlet to perform continuous wavelet transform on the respiratory signal to obtain the wavelet coefficients of the respiratory signal in the time-frequency domain:

[0134] ;;

[0135] In the formula: Represents the wavelet coefficient of the respiratory signal at time point t and frequency f; Represents the complex conjugate of the Morlet wavelet; Represents the respiratory signal, and the input respiratory signal R at time point value;

[0136] Calculate the cross-wavelet spectrum of the respiratory signal and EEG wave, and obtain the cross-wavelet spectrum by multiplying the time-frequency coefficients of the respiratory signal with the complex conjugate wavelet coefficients of the wave signal: ;

[0137] In the formula: is the cross-wavelet spectrum, representing the coupling strength between the respiratory signal and EEG wave at time point t and frequency f; Represents the complex conjugate of the wavelet coefficient of the wave;

[0138] Calculate the normalized wavelet coherence to reflect the synchronization degree between the two signals: ;

[0139] In the formula: Represents the time-frequency smoothing operator; Represents the coherence measure; Represents the modulus value of the cross-wavelet spectrum; Represents the power spectrum of the respiratory signal; Represents EEG power spectrum of the wave;

[0140] Extract the phase information from the cross-wavelet spectrum and calculate the instantaneous phase difference between the respiratory signal and EEG wave at each time-frequency point:

[0141] ; Wherein: represents the phase angle of the respiratory signal at time point t and frequency f; represents EEG the phase angle of the wave at time t and frequency f; represents the instantaneous phase difference; represents taking the phase angle of the complex number;

[0142] Calculate the global phase synchronization index within the frequency band: ;

[0143] Wherein: represents the total number of sampling points within the time window; represents the center frequency of the wave; represents the phase synchronization index, and the closer the value is to 1, the stronger the synchronization; represents the complex exponential representation method for calculating the complex form of the phase synchronization index;

[0144] Based on the above-obtained features, construct a respiration-EEG coupling feature vector. First, calculate the average value of coherence within the frequency band:

[0145] ;

[0146] Wherein: represents the average value of coherence within the frequency band;

[0147] The maximum coherence duration : Identify continuous time periods with coherence exceeding the threshold (0.6) and calculate its maximum duration;

[0148] Phase synchronization stability: Calculate the variance of PSI within the sliding window, with the window length set to 30 seconds and the step size to 10 seconds: ;

[0149] Wherein: represents the number of sliding windows; represents the phase synchronization index within the i-th window; represents the average value of PSI within the sliding window;

[0150] Based on this, output a three-dimensional respiration-EEG coupling feature vector: .

[0151] In this embodiment, by extracting the theta wave component from the EEG signal, the brain activity state of the diver can be accurately captured, especially the changes in mental states such as stress, concentration, and relaxation, providing key electroencephalogram physiological characteristics for diver state monitoring; and by performing continuous wavelet transform on the respiratory signal and calculating its cross-wavelet spectrum with the EEG theta wave, the time-frequency feature coupling between the two can be obtained, revealing the interaction between respiration and brain activity, especially the impact of the respiratory pattern on the diver's brain activity; secondly, by calculating the phase difference between the EEG theta wave and the respiratory signal and extracting the phase synchronization index, the synchronization between respiration and brain electrical activity can be revealed; high synchronization is usually related to characteristics such as the diver's concentration state and mental stability, which is of great significance for monitoring the mental and physiological states of divers.

[0152] Environmental parameter compensation: Based on the depth, establish a relationship to perform depth compensation on SpO2, and compensate the respiratory rate according to the water temperature;

[0153] As a possible implementation manner of this embodiment, the environmental parameter compensation includes the following steps:

[0154] Establish the compensation model: ;

[0155] In the formula: represents the blood oxygen saturation after depth compensation; represents the original blood oxygen saturation measurement value; D represents the underwater depth, with the unit of meter; 0.12 represents the linear compensation coefficient; -0.003 represents the quadratic compensation coefficient, which is used to suppress over-compensation when the depth is greater than 20 meters;

[0156] Verify through an underwater pressure chamber experiment and calculate the model error: ;

[0157] In the formula: represents the measurement value of the arterial blood gas analyzer; N represents the number of measurement samples; it is required that the error after compensation ;

[0158] Linear compensation of the respiratory rate by water temperature:

[0159] Establish the compensation formula for the respiratory rate: ;

[0160] In the formula: represents the compensated respiratory rate; represents the originally measured respiratory rate; represents the current water temperature; 25 represents the reference water temperature; 0.023 represents the temperature sensitivity coefficient;

[0161] When the water temperature is lower than 10°C, the compensation coefficient of the breathing rate is adjusted to 0.015 to avoid overcompensation in the low-temperature zone;

[0162] Multi-parameter coupling effect processing: When the diver is in a deep hypothermic environment (D > 30 m, ), the combined influence of depth and water temperature is more complex. Therefore, a cross-term compensation needs to be introduced:

[0163] ;

[0164] In the formula: represents the cross-compensation coefficient of depth and water temperature;

[0165] To improve the depth measurement accuracy, the depth data of the pressure gauge and sonar sensor are fused using the Kalman filter algorithm: ;

[0166] In the formula: represents the filtered depth value; represents the depth value measured by the pressure gauge; represents the depth value measured by the sonar; represents the Kalman gain.

[0167] In this embodiment, by using a quadratic regression model to perform depth compensation on blood oxygen saturation, the non-linear influence of depth change on blood oxygen saturation can be effectively eliminated, ensuring the accuracy of blood oxygen monitoring, helping the diver to maintain accurate physiological data at different depths, and thus improving the safety monitoring ability of the diver; secondly, a compensation model is established according to the influence of water temperature on breathing rate, so that the breathing rate at different water temperatures can be linearly compensated, thus more accurately reflecting the physiological reactions and needs of the diver; by processing the depth and temperature cross-effect with a Kalman filter, the accuracy of depth estimation is improved, ensuring that the physiological state of the diver in a complex environment can be accurately estimated.

[0168] Short-term risk inference: Based on the non-linear feature vector of HRV and the coupling feature of breathing and EEG, time series features are extracted through the LSTM network structure, and then mapped to short-term risk probabilities based on the fully connected layer to output the real-time risk level;

[0169] Exemplarily, the input data includes the non-linear feature vector of HRV (5 dimensions); the breathing-EEG coupling feature (3 dimensions);

[0170] The time series structure is to construct a time window at a sampling rate of 10 Hz, with a total of 100 time points;

[0171] Z-score normalization is performed on each feature dimension to eliminate the sensor dimension difference;

[0172] Input the standardized data into the LSTM, calculate the hidden state at each time step \(t\), and pass the sequence of hidden states for 100 time steps to the next layer, retaining only the hidden state at the last time step as the time series feature;

[0173] Based on the time series feature output by the LSTM, use a fully connected layer to map it to a risk probability; divide the risk levels based on the output probability value; that is, by setting a threshold, it can be determined which risk level the output probability value belongs to.

[0174] Medium- and long-term trend prediction: Construct a deep causal convolutional layer to predict the future risk evolution, and use the attention mechanism to focus on key time periods to output the risk trend prediction curve;

[0175] In medium- and long-term trend prediction, the input data is a historical feature sequence, which contains 30 time steps, and each time step corresponds to features within 10 seconds; the feature dimension in the historical feature sequence is the same as that in short-term risk inference;

[0176] Construction of the causal convolutional network: ;

[0177] In the formula: \(t\) represents the current time step; represents the feature at the \(t\) -th time step of the output of the \(l\) -th layer; \(K\) represents the size of the convolutional kernel; \(d\) represents the dilation factor; represents the weight matrix of the \(k\) -th convolutional kernel of the \(l\) -th layer; represents the bias term of the \(l\) -th layer;

[0178] Among them, causal convolution ensures that the prediction at each time step only depends on the current and previous time steps, avoiding future information leakage, and the dilation convolution \(d = 2\) to improve the long-range dependence modeling ability of the network by expanding the receptive field;

[0179] Use the output of the last convolutional network to predict the risk trend in the next 3 minutes: ;

[0180] In the formula: represents the output of the last causal convolutional network; represents the projection matrix, which is used to map the output of the convolutional network to 18 time steps (each time step is 10 seconds); represents the bias term; represents the normalization function; through function, the predicted value at each time step is normalized to a probability value, representing the risk probability every ten seconds within the next 3 minutes;

[0181] On this basis, use the attention mechanism for optimization:

[0182] Calculate the attention score: ;

[0183] Wherein: represents a learnable query vector; represents the feature vector at the t-th time step (from the last convolutional network);

[0184] Weighted context based on the calculated attention score: ; Modify the prediction output based on the weighted context: ;

[0185] Wherein: represents the attention weight matrix;

[0186] The final output is the risk probability every ten seconds within the next 3 minutes , where each represents the risk probability at the i-th time step.

[0187] Comprehensive risk index output: Based on the real-time risk level and the risk trend prediction curve, weighted fusion is performed using the compensated environmental parameters to obtain the comprehensive risk index.

[0188] As a possible implementation of this embodiment:

[0189] Deep weight calculation: ;

[0190] Wherein: D represents the diving depth; represents the weight of the depth, which increases with the increase of the depth;

[0191] Body temperature weight calculation: ;

[0192] Wherein: T represents the temperature; represents the weight of the body temperature. When the body temperature is higher than 28°C, it indicates an overheating risk, and when it is lower than 28°C, it indicates a hypothermia risk;

[0193] Evaluate the reliability of the prediction trend and avoid overfitting: ;

[0194] Wherein: represents the prediction mean, which represents the mean of the risk probabilities at all time steps; represents the i-th predicted value;

[0195] Weighted fusion: According to the real-time risk and the prediction trend, fuse and output the final comprehensive risk index: ;

[0196] Wherein: represents the real-time risk level; represents the maximum predicted risk in the next 3 minutes.

[0197] In the present invention, in the data acquisition stage, by aligning the timestamps of the physiological signs signals and the environmental data, the consistency of the multimodal signals is ensured, providing accurate data support for subsequent processing. For the signal artifacts caused by the movement of the diver, the movement features are extracted using the acceleration data, a movement noise model is constructed, and the filter weights are updated in real time through the recursive least squares algorithm, removing the artifacts and ensuring the accuracy of the blood flow signals.

[0198] In the purification of EEG signals, the db4 wavelet basis is used for 6-layer wavelet packet decomposition, the features in the frequency band related to water pressure are extracted, and the low-frequency baseline drift is removed through thresholding processing, enhancing the quality of the EEG signals. The heart rate signals are analyzed through nonlinear dynamics, and the maximum Lyapunov exponent is extracted to evaluate the dynamic changes of HRV, providing key features for short-term risk inference. The coupling features between the respiratory signals and the EEG waves are analyzed using wavelet transform, further improving the accuracy of risk assessment. In the environmental compensation aspect, based on the relationship between depth and SpO2 and the compensation of water temperature for the respiratory rate, the accuracy of the physiological data is ensured. The temporal features are extracted through the LSTM network structure, and combined with the deep causal convolution layer for medium- and long-term risk trend prediction, and the risk level and risk prediction curve of the diver are output in real time. Finally, by integrating the real-time risk level and the trend prediction curve, and combining with the compensation of environmental parameters, a comprehensive risk index is obtained, comprehensively and accurately evaluating the health status and potential risks of the diver, effectively improving the safety of the diver.

[0199] The above are only the preferred embodiments of the present application and are not intended to limit the present application. Any modifications, equivalent replacements, and improvements made within the spirit and principle of the present application shall be included within the protection scope of the present application.

[0200] The above are only the preferred embodiments of the present application and are not intended to limit the present application. Any modifications, equivalent replacements, and improvements made within the spirit and principle of the present application shall be included within the protection scope of the present application.

Claims

1. A method for monitoring the status of divers based on biofeedback, characterized in that, It includes the following steps: Data acquisition: Collect physiological signal data and environmental data, align the collected data by timestamp, and obtain multi-modal signal flow data after time-axis alignment; Motion noise suppression: Extract motion features from the acceleration data in the physiological signal data to construct a motion noise model, and update the filter weights of the motion noise model in real time through the recursive least squares algorithm. Remove artifacts through the motion noise model to obtain blood flow signals after artifact removal; EEG signal enhancement: Perform 6-layer wavelet packet decomposition on the EEG signal data in the physiological signal data using the db4 wavelet basis, extract features in different frequency bands from the EEG signal, and identify the energy distribution in the frequency bands related to water pressure. Remove low-frequency baseline drift through thresholding to obtain purified EEG signals; Time-domain feature generation: Perform non-linear dynamics analysis on the denoised heart rate signal, evaluate the complexity of the signal, construct a phase space reconstruction diagram, extract the maximum Lyapunov exponent to reflect the dynamic changes of the heart rate signal, and obtain the non-linear feature vector of HRV; The specific steps of the above-mentioned time-domain feature generation are as follows: Sample entropy calculation: Based on the denoised heart rate signal, convert the heart rate signal into a high-dimensional vector through phase space reconstruction. The embedding dimension and time delay are used in the reconstruction process, and the time delay is determined by the autocorrelation function method; then calculate the matching degree between high-dimensional vectors to obtain the sample entropy; Maximum Lyapunov exponent extraction: Based on the denoised heart rate signal, construct a phase space trajectory to obtain the trajectory, then find the nearest point for each point, calculate the divergence speed between adjacent points, and further calculate the maximum Lyapunov exponent; Correlation dimension calculation: Based on the denoised heart rate signal, calculate the dimensional characteristics of the phase space through the correlation function to obtain the similarity between two points in the phase space, and then fit the relationship between the two points in the double logarithmic coordinate system to obtain the correlation dimension; Quadratic sample entropy calculation: Based on the denoised heart rate signal, obtain the quadratic sample entropy value through similarity matching; Detrended fluctuation analysis: Based on the denoised heart rate signal, calculate the cumulative deviation sequence to obtain the cumulative deviation of each time point relative to the mean value. Divide the sequence into multiple small intervals, calculate the local volatility within each interval, and finally obtain the DFA scaling exponent through double logarithmic fitting; The results obtained from sample entropy calculation, maximum Lyapunov exponent extraction, correlation dimension calculation, quadratic sample entropy calculation, and detrended fluctuation analysis form a five-dimensional non-linear feature vector; Frequency domain feature generation: Based on the purified EEG signal and respiratory signal, wavelet transform is used to analyze the coherence between the respiratory signal and the EEG wave, and the coupling degree between the respiratory signal and the EEG wave is quantified by the phase synchronization index; the coupling feature between respiration and EEG is obtained; Environmental parameter compensation: Perform depth compensation on SpO2 based on the established relationship by depth, and compensate the respiratory rate according to the water temperature; [[ID=?]]Short-term risk inference: Based on the non-linear feature vector of HRV and the coupling features of respiration and EEG, extract time-series features through the LSTM network structure, and then map them to short-term risk probabilities based on the fully connected layer to output the real-time risk level; Medium- and long-term trend prediction: Construct a deep causal convolution layer to predict the future risk evolution, and focus on key periods through the attention mechanism to output the risk trend prediction curve; Comprehensive risk index output: Based on the real-time risk level and the risk trend prediction curve, the compensated environmental parameters are used for weighted fusion to obtain the comprehensive risk index.

2. The method for monitoring the status of divers based on biofeedback according to claim 1, characterized in that, The motion noise model is as follows: ; Where: represents the motion noise; represents the time-varying weight parameter; represents the filtered acceleration signal; When the time-varying weight parameters are updated in real time through the recursive least squares algorithm, after obtaining the motion noise, the motion noise is subtracted from the original blood flow signal to obtain the blood flow signal after artifact removal.

3. The method for monitoring the status of divers based on biofeedback according to claim 1, characterized in that, The EEG signal enhancement includes the following steps: Wavelet packet decomposition: Using the db4 wavelet basis for 6-layer decomposition, the signal is decomposed into 64 sub-bands, and each band represents the energy of the signal in a specific frequency range; Based on the decomposed sub-bands, the water pressure interference band and the low-frequency drift band are located; Water pressure interference suppression: Energy analysis is performed in the water pressure interference band, and the energy ratio of each band node is calculated. If the energy ratio of the water pressure interference band is higher than a predetermined threshold, the corresponding wavelet packet coefficients are attenuated through soft threshold technology, and the effective components of the EEG signal are retained; Removal of low-frequency baseline drift: For the baseline drift below 0.5 Hz, by performing additional decomposition on the wavelet packet tree, the components in the range of 0-1 Hz are removed from the signal; Wavelet packet reconstruction: Reconstruct each processed frequency point node, and the purified EEG time-domain signal is obtained by synthesizing the wavelet packet coefficients layer by layer; Time-frequency feature matrix generation: Extract the energy within the time window of each node from the processed signal and perform normalization processing to obtain the time-frequency feature matrix.

4. The method for monitoring the status of a diver based on biofeedback according to claim 1, characterized in that The frequency domain feature generation includes the following steps: EEG EEG wave component extraction: Extract EEG wave components from the EEG time-frequency matrix based on neurophysiological characteristics wave components; Continuous wavelet transform of the respiration signal: Use the complex Morlet wavelet to perform continuous wavelet transform on the respiration signal to obtain the wavelet coefficients of the respiration signal in the time-frequency domain, that is, the time-frequency representation of the respiration signal; Cross-wavelet spectrum calculation: Calculate the time-frequency representation of the respiratory signal and the cross-wavelet spectrum of the EEG wave, and obtain the cross-wavelet spectrum by multiplying the time-frequency coefficients of the respiratory signal with the complex conjugate wavelet coefficients of the EEG wave, and calculate the normalized wavelet coherence based on the obtained cross-wavelet spectrum; Phase synchronization index extraction: Extract phase information from the cross-wavelet spectrum, and calculate the instantaneous phase difference between the respiration signal and the EEG waves at each time-frequency point, and calculate the calculated instantaneous phase difference The global phase synchronization index within the frequency band is obtained to get the phase synchronization index; Feature vector generation: Based on the calculation of normalized wavelet coherence The time average of the coherence within the frequency band is obtained to get the average coherence intensity; the continuous time periods with coherence exceeding the threshold are identified, and the maximum duration is calculated; finally, the variance of the phase synchronization index within the sliding window is calculated, and based on this, the 3D coupling features of respiration and EEG are obtained.

5. The method for monitoring the status of a diver based on biofeedback according to claim 1, characterized in that The environmental parameter compensation includes the following steps: Nonlinear compensation of depth on blood oxygen saturation: Use a quadratic regression model to correct the blood oxygen saturation, where the linear term and the quadratic term of the quadratic regression model compensate for the deviation and over-compensation problems caused by depth respectively; Linear compensation of water temperature on respiration frequency: Establish a respiration frequency compensation model according to the Q10 effect, where linear compensation is achieved through the temperature sensitivity coefficient; Processing of multi-parameter coupling effect: Adjust the respiration frequency by considering the cross-action of depth and water temperature, and then use a Kalman filter to fuse the data of the pressure gauge and the sonar, and achieve accurate depth estimation by dynamically adjusting the Kalman gain.

6. The method for monitoring the status of a diver based on biofeedback according to claim 1, characterized in that, The comprehensive risk index output includes the following steps: Definition of membership function: Quantify the impact of environmental parameters on the comprehensive risk through the membership function, where the depth weight gradually increases with the increase of depth, and the body temperature weight is dynamically adjusted according to the degree of deviation of the body temperature from the reference value; Weighted fusion: Perform weighted fusion according to the results of real-time risk and prediction trend, and output the final comprehensive risk index.

Images