







技术领域technical field
本发明涉及数据处理技术领域,尤其涉及一种呼吸数据处理方法。The invention relates to the technical field of data processing, in particular to a method for processing respiratory data.
背景技术Background technique
人体呼吸波形作为一种重要的生理信号,在众多医疗器械或解决方案中使用,如肺部结节穿刺、肺部及其周围脏器的肿瘤放疗、肺部4D CT重建、肺部MRI成像等领域。对于患者类周期性、不稳定呼吸曲线的准确显示或特征提取,为上述领域更高端发展提供契机。As an important physiological signal, the human respiratory waveform is used in many medical devices or solutions, such as lung nodule puncture, tumor radiotherapy of the lung and its surrounding organs, lung 4D CT reconstruction, lung MRI imaging, etc. field. The accurate display or feature extraction of patient-like periodic and unstable breathing curves provides an opportunity for higher-end development in the above fields.
从传感器直接采集的呼吸数据存在零点基准漂移、杂波干扰以及呼吸幅值和周期不稳定情况,无法直接供后续处理使用。The breathing data directly collected from the sensor has zero reference drift, clutter interference, and unstable breathing amplitude and period, which cannot be directly used for subsequent processing.
发明内容SUMMARY OF THE INVENTION
发明目的:针对上述不足,本发明提供一种呼吸数据处理方法,处理过的数据能够对肺部穿刺手术提供呼吸引导功能,可以找到相对准确的呼吸周期内的呼吸曲线,为医生提供穿刺指导。Purpose of the invention: In view of the above deficiencies, the present invention provides a breathing data processing method. The processed data can provide a breathing guidance function for pulmonary puncture surgery, and can find a relatively accurate breathing curve within the breathing cycle to provide puncture guidance for doctors.
技术方案:Technical solutions:
一种呼吸数据处理方法,包括:A breathing data processing method, comprising:
获取患者呼吸数据;Obtain patient breathing data;
遍历患者呼吸数据识别得到潜在波峰和波谷;Identify potential peaks and valleys by traversing patient breathing data;
对比相邻波峰、波谷的间隔和幅值剔除伪波峰和伪波谷,得到患者呼吸数据的真实的波峰波谷;Compare the intervals and amplitudes of adjacent peaks and troughs to eliminate false peaks and troughs, and obtain the real peaks and troughs of the patient's breathing data;
在真实波谷处识别呼气末停留阶段。The end-expiratory dwell phase is identified at the true trough.
在识别潜在波峰和波谷之前包括异常值检测替换步骤:Include an outlier detection replacement step before identifying potential peaks and troughs:
检测的异常值Xout范围如下:The detected outliers Xout range are as follows:
xi>Q3+c*(IQR)∪xi<Q1-c*(IQR)xi >Q3+c*(IQR)∪xi <Q1-c*(IQR)
其中,xi表示采集的患者呼吸数据中第i个采样点的呼吸数据,i=1,2,…,n;n表示采集的患者呼吸数据总的采样点数量;Q1、Q3分别表示患者呼吸数据中的幅值通过四分位差得到的第一和第三个四分位数;c为阈值因数,c>0,IQR=Q3-Q1;Among them, xi represents the respiration data of the ith sampling point in the collected patient respiration data, i=1,2,...,n; n represents the total number of sampling points of the collected patient respiration data; Q1 and Q3 represent the patient's respiration respectively The first and third quartiles of the amplitude in the data obtained by the interquartile difference; c is the threshold factor, c>0, IQR=Q3-Q1;
对异常值进行插值替换,异常值范围的两端点以近邻点替换,异常值范围内的点以线性插值替换,如下式:The outliers are replaced by interpolation, the two ends of the outlier range are replaced by the nearest neighbors, and the points within the outlier range are replaced by linear interpolation, as follows:
其中,k表示第k个采样点,m代表中间异常点个数。Among them, k represents the kth sampling point, and m represents the number of intermediate abnormal points.
在所述异常值检测替换步骤之前包括滤波步骤。A filtering step is included before the outlier detection replacement step.
在识别潜在波峰和波谷之前还包括基准校准步骤:A reference calibration step is also included before identifying potential peaks and valleys:
通过下式拟合一阶多项式;Fit a first-order polynomial by:
得到线性趋势:Get a linear trend:
f(i)=w*i+bf(i)=w*i+b
其中,w、b分别表示一阶多项式的系数;Among them, w and b respectively represent the coefficients of the first-order polynomial;
对呼吸数据进行基准校准:xi′=xi-f(i)。The respiration data is benchmarked: xi '=xi -f(i).
遍历患者呼吸数据识别得到潜在波峰和波谷具体为:The potential peaks and troughs obtained by traversing the patient's breathing data are as follows:
根据呼吸数据采样点数设置若干滑窗和滑窗起始位置,采用每个滑窗依次遍历患者呼吸数据,若其中某一采样点的呼吸数据中的呼吸幅值满足波峰和波谷的接受阈值,则认为其为波峰或波谷,从而得到每个滑窗遍历完呼吸数据之后的波峰数目和波谷数目;Set a number of sliding windows and the starting position of the sliding window according to the number of respiratory data sampling points, and use each sliding window to traverse the patient's respiratory data in turn. Consider it as a peak or a trough, so as to obtain the number of peaks and troughs after each sliding window traverses the respiratory data;
依次计算每个滑窗相对于上一个滑窗遍历完呼吸数据得到波峰数目和波谷数目的增量,并计算得到增量对应的多个波峰和多个波谷对应的呼吸幅值的平均值,并以两平均值的均值作为分割阈值;Calculate the increments of the number of peaks and troughs obtained by traversing the respiratory data of each sliding window relative to the previous sliding window in turn, and calculate the average value of the respiratory amplitudes corresponding to multiple peaks and troughs corresponding to the increments, and Take the mean of the two averages as the segmentation threshold;
以呼吸数据中呼吸幅值大于该分割阈值的波峰作为潜在波峰,以呼吸数据中呼吸幅值小于该分割阈值的波谷作为潜在波谷。The peaks in the respiration data whose respiration amplitude is greater than the segmentation threshold are used as potential peaks, and the troughs in the respiration data whose respiration amplitudes are smaller than the segmentation threshold are used as potential troughs.
所述接受阈值计算如下:The acceptance threshold is calculated as follows:
计算呼吸数据的均值M和均方差S,设置波峰接受阈值为T1>M+S/2、波谷接受阈值为T2<M-S/2。Calculate the mean M and mean square error S of the respiratory data, and set the peak acceptance threshold as T1 >M+S/2 and the trough acceptance threshold as T2 <MS/2.
按照斐波那契数列设置各个滑窗尺寸,且至少存在一个尺寸大于呼吸周期的滑窗。Set the size of each sliding window according to the Fibonacci sequence, and there is at least one sliding window whose size is greater than the breathing period.
还包括对呼吸数据进行镜像扩展的步骤:Also includes steps to mirror the respiration data:
镜像扩展的尺寸为呼吸数据的采样尺寸和该滑窗遍历至呼吸数据末端时用于补足该滑窗尺寸所需要的尺寸中的较小值。The size of the mirror extension is the smaller of the sampling size of the respiration data and the size required to make up the size of the sliding window when the sliding window traverses to the end of the respiration data.
剔除伪波峰和伪波谷具体为:Removing false peaks and false valleys is as follows:
计算相邻波峰和波谷之间的间隔和相邻波峰和波峰之间的间隔,若后者大于前者,则认为当前波峰为真波峰,否则认为存在伪波峰;若当前波峰对应的呼吸幅值小于下一波峰对应的呼吸幅值,则认为当前波峰为伪波峰;否则认为下一波峰为伪波峰;将伪波峰直接删除;Calculate the interval between adjacent crests and troughs and the interval between adjacent crests and crests. If the latter is greater than the former, the current crest is considered to be a true crest, otherwise, a false crest exists; if the respiratory amplitude corresponding to the current crest is less than If the respiratory amplitude corresponding to the next wave peak, the current wave peak is regarded as a pseudo wave peak; otherwise, the next wave peak is regarded as a pseudo wave peak; the pseudo wave peak is directly deleted;
计算相邻波谷和波峰之间的间隔和相邻波谷和波谷之间的间隔,若后者大于前者,则认为当前对应的波谷为真波谷,否则认为存在伪波谷;若当前波谷对应的呼吸幅值大于下一波谷对应的呼吸幅值,则认为当前波谷为伪波谷;否则认为下一波谷为伪波谷;将伪波谷直接删除。Calculate the interval between adjacent troughs and peaks and the interval between adjacent troughs and troughs. If the latter is greater than the former, the current corresponding trough is considered to be a true trough; otherwise, a false trough exists; if the respiratory amplitude corresponding to the current trough is If the value is greater than the breathing amplitude corresponding to the next trough, the current trough is considered as a pseudo-trough; otherwise, the next trough is considered as a pseudo-trough; the pseudo-trough is directly deleted.
所述在真实波谷处识别呼气末停留阶段具体为:The identification of the end-expiratory dwell stage at the real trough is specifically:
截取相邻两波峰中间部分区段的呼吸数据,以呼吸幅值为纵坐标,以采样点频数为横坐标进行直方图统计;Intercept the respiratory data in the middle part of the two adjacent peaks, take the respiratory amplitude as the ordinate and the frequency of the sampling points as the abscissa for histogram statistics;
某一直方图每个数组bin的范围内的采样点频数为Nbin和其呼吸幅值的中心值为Vbin,计算直方图内采样点频数最大的数组的采样点频数max(Nbin)和其呼吸幅值的中心值Vmaxbin,以Vmaxbin<facmag*range(Vbin)+min(Vbin)及max(Nbin)>facbin1*T(Nbin)作为存在呼气末停留阶段的筛选条件;其中,range(Vbin)表示所有截取的部分呼吸数据的幅值范围,facmag表示直方图中幅值的缩放因子;min(Vbin)表示直方图所有数组的中心值的最小值;T(Nbin)为根据采样频率设置的阈值;facbin1表示直方图数组bin数量的缩放因子;The sampling point frequency within the range of each array bin of a certain histogram is Nbin and the center value of its respiration amplitude is Vbin , and the sampling point frequency max(Nbin ) and The central value Vmaxbin of its respiration amplitude takes Vmaxbin <facmag *range(Vbin )+min(Vbin ) and max(Nbin )>facbin1 *T(Nbin ) as the end-expiratory dwell stage where, range(Vbin ) represents the amplitude range of all intercepted partial respiratory data, facmag represents the scaling factor of the amplitude in the histogram; min(Vbin ) represents the minimum value of the center value of all the arrays in the histogram value; T(Nbin ) is the threshold set according to the sampling frequency; facbin1 represents the scaling factor of the number of bins in the histogram array;
若不满足以上条件,则认为不存在所需要识别的呼气末停留阶段;If the above conditions are not met, it is considered that there is no end-expiratory dwell phase that needs to be identified;
若满足以上条件,则对某一数组bin的采样点数量Nbin在直方图中向上检索,得到满足Nbin>facbin2*Nmaxbin所在数组bin对应的中心值Tbin作为对呼吸数据分割的阈值,facbin2表示直方图数组bin数量的缩放因子;凡是呼吸幅值小于Tbin的呼吸数据认定为呼气末所在阶段点,从而得到呼气末停留阶段。If the above conditions are met, the number of sampling points Nbin of a certain array bin is searched upward in the histogram, and the center value Tbin corresponding to the array bin where Nbin >facbin2 *Nmaxbin is obtained is obtained as the threshold for dividing the respiratory data. , facbin2 represents the scaling factor of the number of bins in the histogram array; any breathing data whose breathing amplitude is less than Tbin is identified as the stage point at the end of expiration, thus obtaining the end-expiratory dwell stage.
有益效果:本发明通过对原始1维数据的滤波降噪、消除线性趋势归零、峰谷搜索识别及呼气末短暂暂停点识别得到有效信号特征,并据此得到一个相对准确的呼吸周期内的呼吸曲线,以此供医生进行肺部穿刺手术指导,提高手术精准度。Beneficial effects: the present invention obtains effective signal characteristics through filtering and noise reduction of the original 1-dimensional data, eliminating the linear trend and returning to zero, searching and identifying peaks and valleys, and identifying the short pause point at the end of expiration, and accordingly obtains a relatively accurate breathing cycle. The breathing curve can be used for doctors to guide the lung puncture operation and improve the accuracy of the operation.
附图说明Description of drawings
图1为呼吸数据处理方法的流程图;Fig. 1 is the flow chart of respiratory data processing method;
图2为原始呼吸数据的示意图;Fig. 2 is the schematic diagram of raw respiration data;
图3为滤波后呼吸数据的示意图;Fig. 3 is the schematic diagram of respiratory data after filtering;
图4为消除线性趋势的呼吸数据的示意图;Fig. 4 is the schematic diagram of the respiration data that eliminates linear trend;
图5为波峰与波谷识别的流程图;Fig. 5 is the flow chart of wave crest and wave trough identification;
图6为搜索得到的呼吸数据中波峰与波谷的示意图;Fig. 6 is the schematic diagram of wave crest and wave trough in the respiration data that search obtains;
图7为呼气末停留阶段识别的示意图;Fig. 7 is the schematic diagram of end-expiratory dwell phase identification;
图8为搜索呼气末端停留阶段的示意图。FIG. 8 is a schematic diagram of searching for a dwell phase at the end of exhalation.
具体实施方式Detailed ways
下面结合附图和具体实施例,进一步阐明本发明。The present invention will be further illustrated below in conjunction with the accompanying drawings and specific embodiments.
本发明的呼吸数据处理方法如图1所示,包括如下步骤:The respiratory data processing method of the present invention, as shown in Figure 1, includes the following steps:
(1)获取患者原始呼吸数据;(1) Obtain the original respiratory data of the patient;
患者的呼吸波形数据可以来自各种呼吸数据采集设备,从人体呼吸端来讲,如面罩式呼吸流量阀、鼻腔温度场等,从呼吸引起的体表变化来讲,如空气压力式腹带、ECG心电信号、短波雷达UWB、光学导航、磁导航等。The patient's breathing waveform data can come from various respiratory data acquisition devices. From the perspective of human breathing, such as mask-type breathing flow valve, nasal temperature field, etc., from the surface changes caused by breathing, such as air pressure abdominal belt, ECG electrocardiogram signal, short-wave radar UWB, optical navigation, magnetic navigation, etc.
(2)对步骤(1)得到的患者原始呼吸数据的异常值进行处理;(2) processing the abnormal value of the patient's original respiratory data obtained in step (1);
从呼吸数据采集设备得到的原始呼吸数据为序列数据,如图2所示,因人状态而异,如咳嗽、紧张、抖动等会对数据造成干扰,其存在局部极值、杂波干扰等,需要对呼吸数据执行滤波及异常值检测替换处理,然后提取到特征信号才能更好的应用数据,如呼吸门控、模式识别等,具体处理如下:The original respiratory data obtained from the respiratory data acquisition equipment is sequence data, as shown in Figure 2, which varies with the state of the person, such as coughing, nervousness, shaking, etc., which will interfere with the data, and there are local extreme values, clutter interference, etc. It is necessary to perform filtering and outlier detection and replacement processing on the respiratory data, and then extract the characteristic signal to better apply the data, such as respiratory gating, pattern recognition, etc. The specific processing is as follows:
(21)数据平滑滤波;(21) data smoothing filtering;
鉴于频率域滤波存在相位偏移问题,本发明在异常值检测替换之前采用Savitzky-Golay滤波,如图3所示;In view of the phase offset problem in frequency domain filtering, the present invention adopts Savitzky-Golay filtering before outlier detection and replacement, as shown in FIG. 3 ;
(22)异常值检测替换;(22) Outlier detection and replacement;
本发明检测的异常值Xout为如下范围:The abnormal value Xout detected by the present invention is the following range:
xi>Q3+c*(IQR)∪xi<Q1-c*(IQR) (1)xi >Q3+c*(IQR)∪xi <Q1-c*(IQR) (1)
其中,xi表示采集的序列数据中第i个采样点的呼吸数据,i=1,2,…,n;n表示采集的序列数据总的采样点数量;Q1、Q3分别表示序列数据中的幅值通过四分位差得到的第一和第三个四分位数;c为阈值因数,c>0,一般取1.5;IQR=Q3-Q1;Among them, xi represents the respiration data of the i-th sampling point in the collected sequence data, i=1,2,...,n; n represents the total number of sampling points of the collected sequence data; Q1 and Q3 respectively represent the number of sampling points in the sequence data. The first and third quartiles of the amplitude obtained by the interquartile difference; c is the threshold factor, c>0, generally 1.5; IQR=Q3-Q1;
通过该步骤可以对步骤(21)滤波后得到的呼吸数据中急剧变化点进行检测替换,其中式(1)检测得到了滤波后呼吸数据异常值,采用式(2)对检测得到的异常值进行插值替换,其中,插值替换具体为:异常值范围的两端点以近邻点替换,异常值范围内的点以线性插值替换,通过式(2)表示:Through this step, the sharp change points in the respiratory data obtained after the filtering in step (21) can be detected and replaced, wherein the abnormal value of the filtered respiratory data is detected by the formula (1), and the detected abnormal value is obtained by using the formula (2). Interpolation replacement, where the interpolation replacement is specifically: the two ends of the outlier range are replaced by the nearest neighbor points, and the points in the outlier range are replaced by linear interpolation, which is expressed by formula (2):
其中,k表示第k个采样点,m代表中间异常点个数;Among them, k represents the kth sampling point, and m represents the number of intermediate abnormal points;
通过该步骤可以对步骤(1)得到的呼吸数据中进行平稳性检查,并进行异常值数据替换;Through this step, the respiration data obtained in step (1) can be checked for stationarity, and abnormal value data can be replaced;
(3)进行基准校准,消除线性趋势;(3) Perform benchmark calibration to eliminate linear trends;
呼吸采样设备长时间采集数据过程中,人体移动等导致基准逐步发生偏移,因此呼吸数据中存在因体位变动、低频干扰带来的趋势变化,因此需要消除这部分趋势变化对呼吸数据的影响,需要将数据转换到零值附近,使得数据更准确。During the long-term data collection process by the respiratory sampling device, the movement of the human body causes the reference to gradually shift. Therefore, there are trend changes in the respiratory data caused by body position changes and low-frequency interference. Therefore, it is necessary to eliminate the impact of these trend changes on the respiratory data. The data needs to be converted to around zero to make the data more accurate.
如图4所示,本发明先采用一阶多项式拟合计算数据的趋势线,之后对降噪的数据减去趋势线得到基准为零的无零飘的数据信号,具体如下:As shown in FIG. 4 , the present invention first uses a first-order polynomial to fit the trend line of the calculated data, and then subtracts the trend line from the noise-reduced data to obtain a zero-zero data signal with zero drift, as follows:
通过下式拟合一阶多项式;Fit a first-order polynomial by:
得到线性趋势:Get a linear trend:
f(i)=w*i+bf(i)=w*i+b
其中,w、b分别表示一阶多项式的系数;Among them, w and b respectively represent the coefficients of the first-order polynomial;
从而可以通过下式对经步骤(2)处理后的呼吸数据消除线性趋势得到准确的呼吸数据:Thus, accurate respiratory data can be obtained by eliminating the linear trend of the respiratory data processed by step (2) by the following formula:
xi′=xi-f(i)xi ′=xi -f(i)
(4)识别潜在波峰和波谷;(4) Identify potential peaks and valleys;
搜索经步骤(3)处理的呼吸数据的波峰和波谷,得到对应的呼气末和吸气末,从而明确呼吸数据中具体的呼气阶段(同一周期内波谷至波峰阶段)和吸气阶段(同一周期内波峰至波谷阶段);通常的波峰波谷二阶差分识别算法,即将局部极值点识别为呼气末或吸气末,容易出现误识情况,本发明基于滑窗算法识别波峰波谷,可以大大提高识别准确率,如图5所示,具体如下:Search the peaks and troughs of the respiratory data processed in step (3) to obtain the corresponding end-expiration and end-inspiration, so as to clarify the specific expiratory phase (trough to peak phase in the same cycle) and inspiratory phase ( Peak to trough stage in the same cycle); the usual second-order differential identification algorithm of peak and trough, that is to identify the local extreme point as the end of expiration or the end of inspiration, is prone to misrecognition, and the present invention identifies the peak and trough based on the sliding window algorithm, The recognition accuracy can be greatly improved, as shown in Figure 5, as follows:
(41)计算呼吸数据的均值M和均方差S,设置波峰接受阈值T1>M+S/2,波谷接受阈值T2<M-S/2;(41) Calculate the mean value M and the mean square error S of the respiratory data, set the peak acceptance threshold T1 >M+S/2, and the trough acceptance threshold T2 <MS/2;
(42)根据呼吸数据采样点数设置若干滑窗和滑窗起始位置,其中,可按照斐波那契数列设置各个滑窗尺寸得到序列滑窗,且必须存在大于呼吸周期的尺寸的滑窗,其中,滑窗尺寸即呼吸数据的采样点数,也即呼吸数据的采样尺寸的某一段尺寸;通过每个滑窗尺寸对呼吸数据进行遍历时,可根据实际需求从不同的滑窗起始位置多次遍历呼吸数据,滑窗起始位置可根据滑窗尺寸设置;保证每个滑窗尽可能多的覆盖数据因周期、幅值不稳定导致的局部波峰波谷。(42) Set several sliding windows and sliding window starting positions according to the number of respiratory data sampling points, wherein, the size of each sliding window can be set according to the Fibonacci sequence to obtain a sequence sliding window, and there must be a sliding window larger than the size of the breathing period, Among them, the size of the sliding window is the number of sampling points of the breathing data, that is, the size of a certain segment of the sampling size of the breathing data; when traversing the breathing data through each sliding window size, the starting position of different sliding windows can be selected according to actual needs. After traversing the respiratory data for the second time, the starting position of the sliding window can be set according to the size of the sliding window; ensure that each sliding window covers as much as possible the local peaks and valleys of the data due to unstable periods and amplitudes.
其中,滑窗的尺寸与采样频率正相关,如波峰-波谷对的数据个数为50,可将滑窗尺寸分别设置10、20、30、50、80、130、210,保证波峰-波谷对在滑窗的中心位置;分别设置滑窗起始位置为0%、10%、30%等;Among them, the size of the sliding window is positively related to the sampling frequency. For example, the number of data of the peak-valley pair is 50, and the sliding window size can be set to 10, 20, 30, 50, 80, 130, and 210 respectively to ensure the peak-valley pair. At the center of the sliding window; set the starting position of the sliding window as 0%, 10%, 30%, etc.;
为防止滑窗遍历发生数据不足导致滑窗无法完全遍历呼吸数据的情况,必要时需对呼吸数据进行镜像扩展,扩展的呼吸数据为从呼吸数据采样点的0点开始的呼吸数据,扩展的尺寸h为呼吸数据的采样尺寸和该滑窗遍历至呼吸数据末端时用于补足该滑窗尺寸以完整遍历呼吸数据所需要的尺寸中的较小值;则将呼吸数据{xi}扩展成扩展呼吸数据{Xj},j=1,2,…,j,j=n+h;In order to prevent the sliding window from traversing the respiratory data due to insufficient data in the sliding window traversal, the respiratory data needs to be mirrored and expanded if necessary. The expanded respiratory data is the respiratory data starting from the 0 point of the respiratory data sampling point. h is the smaller value of the sampling size of the respiration data and the size required to supplement the sliding window size to complete the traversal of the respiration data when the sliding window traverses to the end of the respiration data; then the respiration data {xi } is expanded into extended respiration data {Xj }, j=1,2,...,j,j=n+h;
(43)采用每个滑窗依次遍历经步骤(42)扩展后的扩展呼吸数据,若其中某一采样点的呼吸数据中的呼吸幅值满足阈值T1或T2,则认为其为波峰或波谷,并增加对应的索引数目,具体的:某一采样点的呼吸数据大于阈值T1,则认为其为波峰,并将波峰数目增1,某一采样点的呼吸数据小于阈值T2,则认为其为波谷,并将波谷数目增1,从而得到每个滑窗遍历完扩展呼吸数据之后的波峰数目Nap和波谷数目Nat;(43) Using each sliding window to traverse the expanded respiratory data expanded in step (42) in turn, if the respiratory amplitude in the respiratory data at a certain sampling point satisfies the threshold T1 or T2 , it is considered to be a peak or trough, and increase the corresponding index number, specifically: if the respiratory data of a certain sampling point is greater than the threshold T1 , it is considered as a peak, and the number of peaks is increased by 1, and the respiratory data of a certain sampling point is less than the threshold T2 , then Consider it a trough, and increase the number of troughs by 1, so as to obtain the number of peaks Nap and the number of troughs Nat after each sliding window traverses the extended respiratory data;
(44)依次计算每个滑窗相对于上一个滑窗遍历完扩展呼吸数据得到波峰数目和波谷数目的增量,如有N个滑窗,则各有N-1个潜在波峰数目和波谷数目的增量;从而计算得到增量对应的多个波峰和多个波谷对应的呼吸幅值的平均值,并以两平均值的均值作为分割阈值T,即:(44) Calculate the increments of the number of peaks and troughs obtained by traversing the extended breathing data of each sliding window relative to the previous sliding window in turn. If there are N sliding windows, each has N-1 potential number of peaks and troughs. Therefore, the average value of the respiratory amplitudes corresponding to multiple peaks and multiple troughs corresponding to the increment is calculated, and the average of the two averages is used as the segmentation threshold T, namely:
T={median[diff(Nap)]+median[diff(Nat)]}*0.5T={median[diff(Nap)]+median[diff(Nat)]}*0.5
其中,diff()表示求取两相邻滑窗中波峰或波谷数目的增量的函数,median()表示求取其内增量对应的呼吸幅值的平均值;Wherein, diff() represents a function to obtain the increment of the number of peaks or troughs in two adjacent sliding windows, and median() represents the average value of the respiratory amplitude corresponding to the increment within it;
(45)根据步骤(44)得到的分割阈值,以步骤(43)得到的呼吸数据中呼吸幅值大于该分割阈值的波峰作为潜在波峰,以步骤(43)得到的呼吸数据中呼吸幅值小于该分割阈值的波谷作为潜在波谷;(45) According to the segmentation threshold obtained in step (44), take the peak whose breathing amplitude is greater than the segmentation threshold in the breathing data obtained in step (43) as a potential peak, and take the breathing data obtained in step (43) whose breathing amplitude is less than The trough of the segmentation threshold is regarded as a potential trough;
若潜在波峰中第一个波峰对应的呼吸幅值小于潜在波谷中任一个波谷对应的呼吸幅值,则删除潜在波峰中的第一个波峰以保证波峰和波谷对,并初始化波峰和波谷的索引pt=0和tt=0。If the respiration amplitude corresponding to the first peak in the potential peak is smaller than the respiration amplitude corresponding to any trough in the potential trough, delete the first peak in the potential peak to ensure the peak and trough pair, and initialize the index of the peak and trough pt =0 and tt =0.
(5)剔除潜在波峰和波谷中的伪波峰和伪波谷;(5) Eliminate false peaks and troughs in potential peaks and troughs;
(51)识别伪波峰:(51) Identify false peaks:
计算相同索引的波峰和波谷(即相邻波峰和波谷)之间的间隔d1=vt(tt)-vt(pt)和相邻索引的波峰和波峰(即相邻波峰)之间的间隔d2=vt(pt+1)-vt(pt),其中vt()表示某一索引的采样点时刻,vp()表示某一索引的原始数据采样点的呼吸幅值;Calculate the interval between peaks and troughs of the same index (ie adjacent peaks and troughs) d1 =vt(tt )-vt(pt ) and the interval between peaks and peaks of adjacent indices (ie adjacent peaks) Interval d2 =vt(pt+1 )-vt(pt ), wherein vt( ) represents the sampling point moment of a certain index, and vp( ) represents the respiration amplitude of the original data sampling point of a certain index;
若d2>d1,则认为当前波峰为真波峰,执行下一步,否则认为其为伪波峰并对其分类处理,具体为:若vp(pt)<vp(pt+1),则认为索引pt对应的波峰为伪波峰,直接跳过;否则认为索引pt+1对应的波峰为伪波峰,直接删除索引pt+1对应的波峰的呼吸数据;If d2 >d1 , the current peak is considered to be a true peak, and the next step is performed; otherwise, it is considered to be a false peak and classified and processed, specifically: if vp(pt )<vp(pt+1 ), then The peak corresponding to the indexpt is considered to be a pseudo peak, and skipped directly; otherwise, the peak corresponding to the index pt+1 is considered to be a pseudo peak, and the respiration data of the peak corresponding to the index pt+1 is directly deleted;
重复本步骤,直至对所有潜在波峰遍历完毕;Repeat this step until all potential peaks are traversed;
(52)识别伪波谷;(52) Identifying false troughs;
计算相邻索引的波谷和波峰(同样即相邻波峰和波谷)之间的间隔d3=vt(tt+1)-vt(tt)和相邻索引的波谷和波谷(即相邻波谷)之间的间隔d4=vt(tt+1)-vt(tt),若d4>d3,则认为索引tt对应的波谷为真波谷,执行下一步,否则认为其为伪波谷并对其分类处理,具体为:Calculate the interval d3 =vt(tt+1 )-vt(tt ) between adjacent indexed troughs and peaks (again, adjacent peaks and troughs) and adjacent indexed troughs and troughs (i.e. adjacent troughs) ) between d4 =vt(tt+1 )-vt(tt ), if d4 >d3 , the trough corresponding to the index tt is considered to be a true trough, and the next step is performed, otherwise it is considered a false one The trough is classified and processed, specifically:
若vp(tt)>vp(tt+1),则认为索引tt对应的波谷为伪波谷,直接跳过;否则认为索引tt+1对应的波谷为伪波谷,直接删除索引tt+1对应的波谷的呼吸数据;If vp(tt )>vp(tt+1 ), then the trough corresponding to the index tt is considered to be a false trough and skipped directly; otherwise, the trough corresponding to the index tt+1 is considered to be a false trough, and the index tt is directly deleted+1 for the respiration data corresponding to the trough;
重复本步骤,直至对所有潜在波峰和潜在波谷遍历完毕;Repeat this step until all potential peaks and valleys are traversed;
(53)将经过步骤(51)和(52)剔除伪波峰和伪波谷后剩余的波峰和波谷作为真实波峰和真实波谷;(53) will pass through steps (51) and (52) to remove the remaining peaks and valleys after false peaks and false valleys as real peaks and true valleys;
遍历完潜在的波峰波谷后,保留未经镜像扩展的呼吸数据的索引对应的真实波峰波谷即为实际的波峰波谷值,如图6所示。After traversing the potential peaks and troughs, the real peaks and troughs corresponding to the indexes of the respiration data without mirror expansion are retained as the actual peaks and troughs, as shown in FIG. 6 .
(6)搜索呼气末停留阶段;(6) Search for the end-expiratory dwell stage;
人体吸气时间比呼气时间占比小,正常呼吸存在呼气停留阶段,也可能因为呼吸节奏快跳过这个过程,因此需要判断是否存在一个相对较长的呼气停留阶段并进行识别,具体的:截取经步骤(5)得到的呼吸数据的相邻两波峰中间部分区段的部分呼吸数据,对该部分数据做直方图统计,如图7所示,以呼吸幅值为纵坐标、以采样点频数为横坐标进行直方图统计,直方图的数组的数量根据实际需求确定;Human inhalation time is smaller than exhalation time. Normal breathing has an expiratory stop phase, and this process may be skipped because of the fast breathing rhythm. Therefore, it is necessary to determine whether there is a relatively long expiratory stop phase and identify it. : intercepting part of the respiratory data in the middle part of the adjacent two peaks of the respiratory data obtained in step (5), and making histogram statistics on this part of the data, as shown in Figure 7, with the respiratory amplitude as the ordinate and the The frequency of sampling points is the abscissa for histogram statistics, and the number of arrays in the histogram is determined according to actual needs;
直方图中每个数组bin对应的采样点频数为Nbin、呼吸幅值的中心值为Vbin,计算直方图内采样点频数最大的数组的采样点频数Nmaxbin和呼吸幅值的中心值Vmaxbin,以Vmaxbin<facmag*range(Vbin)+min(Vbin)及Nmaxbin>facbin1*T(Nbin)作为存在呼气末停留阶段的筛选条件,若不满足以上条件,则认为不存在所需要识别的呼气末停留阶段,其中,range(Vbin)表示所截取的部分呼吸数据的幅值范围,facmag表示直方图中幅值的缩放因子,推荐取0.3;min(Vbin)表示直方图所有数组的中心值的最小值;T(Nbin)为根据采样频率设置的阈值,采样频率越高,此阈值越大,本申请中,采样频率为30赫兹,T(Nbin)设为5;facbin1表示直方图数组bin数量的缩放因子,推荐取0.9;The frequency of sampling points corresponding to each array bin in the histogram is Nbin , and the center value of respiration amplitude is Vbin . Calculate the frequency of sampling points Nmaxbin of the array with the largest sampling point frequency in the histogram and the center value V of respiration amplitude.maxbin , take Vmaxbin <facmag *range(Vbin )+min(Vbin ) and Nmaxbin >facbin1 *T(Nbin ) as the screening conditions for the existence of the end-expiratory dwell stage, if the above conditions are not met, then It is considered that there is no end-expiratory stay stage that needs to be identified, where range(Vbin ) represents the amplitude range of the intercepted partial respiratory data, facmag represents the scaling factor of the amplitude in the histogram, and it is recommended to take 0.3; min( Vbin ) represents the minimum value of the center value of all the arrays in the histogram; T(Nbin ) is the threshold set according to the sampling frequency. The higher the sampling frequency, the larger the threshold. In this application, the sampling frequency is 30 Hz, and T ( Nbin ) is set to 5; facbin1 represents the scaling factor of the number of bins in the histogram array, and it is recommended to take 0.9;
当满足以上条件(即存在需要识别的呼气停留阶段)时,对某一数组bin的采样点数量Nbin在直方图中向上检索,得到满足Nbin>facbin2*Nmaxbin所在数组bin对应的中心值Tbin作为对呼吸数据分割的阈值,facbin2表示直方图数组bin数量的缩放因子,此处推荐取0.75;凡是呼吸幅值小于Tbin的呼吸数据均认定为呼气末所在阶段点,从而得到呼气末停留阶段,如图8所示;When the above conditions are met (that is, there is an expiratory stop stage that needs to be identified), the number of sampling points Nbin of a certain array bin is searched upward in the histogram, and the corresponding bin of the array bin where Nbin > facbin2 *Nmaxbin is obtained is obtained. The central value Tbin is used as the threshold for dividing the respiratory data, facbin2 represents the scaling factor of the number of bins in the histogram array, and it is recommended to take 0.75 here; any respiratory data whose respiratory amplitude is less than Tbin is regarded as the end-expiratory stage point, Thus, the end-expiratory dwell phase is obtained, as shown in Figure 8;
(7)指导穿刺;(7) guide puncture;
通过以上步骤对患者呼吸数据进行处理之后,能够得到干扰较小的患者呼吸波形,以该波形通过步骤(6)搜索患者的呼气末停留阶段并据此指导穿刺,能够实现较为精准的穿刺手术。After the patient's breathing data is processed through the above steps, the patient's breathing waveform with less interference can be obtained, and the waveform can be used to search for the patient's end-expiratory stay stage through step (6) and guide the puncture accordingly, which can realize a more accurate puncture operation. .
本发明中,在获取得到患者原始呼吸数据后,可以无需对齐进行异常值处理和基准校准,直接识别其中潜在波峰和波谷;当然也可以进行异常值处理或基准校准,然后再识别其中潜在波峰和波谷,虽然相比前一方法效果较差,但也可以取得本发明的技术效果,解决本发明所要解决的技术问题。In the present invention, after obtaining the original respiratory data of the patient, outlier processing and benchmark calibration can be performed without alignment, and potential peaks and troughs can be directly identified; of course, outlier processing or benchmark calibration can also be performed, and then potential peaks and troughs can be identified. Although the trough is less effective than the previous method, the technical effect of the present invention can also be obtained, and the technical problem to be solved by the present invention can be solved.
本发明通过对原始1维数据的滤波降噪、消除线性趋势归零、峰谷搜索识别及呼气末短暂暂停点识别得到有效信号特征,并据此得到一个相对准确的呼吸周期内的呼吸曲线,以此供医生进行肺部穿刺手术指导,提高手术精准度。The present invention obtains effective signal characteristics through filtering and noise reduction of the original 1-dimensional data, eliminating linear trends and returning to zero, peak and valley search and identification, and identification of short pause points at the end of exhalation, and accordingly obtains a relatively accurate breathing curve in the breathing cycle. , so as to provide guidance for doctors to perform lung puncture surgery and improve the accuracy of surgery.
以上详细描述了本发明的优选实施方式,但是本发明并不限于上述实施方式中的具体细节,在本发明的技术构思范围内,可以对本发明的技术方案进行多种等同变换(如数量、形状、位置等),这些等同变换均属于本发明的保护范围。The preferred embodiments of the present invention are described in detail above, but the present invention is not limited to the specific details of the above-mentioned embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations (such as quantity, shape, etc.) can be performed on the technical solutions of the present invention. , position, etc.), these equivalent transformations all belong to the protection scope of the present invention.
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210610264.3ACN114938951B (en) | 2022-05-31 | 2022-05-31 | Respiration data processing method |
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210610264.3ACN114938951B (en) | 2022-05-31 | 2022-05-31 | Respiration data processing method |
| Publication Number | Publication Date |
|---|---|
| CN114938951Atrue CN114938951A (en) | 2022-08-26 |
| CN114938951B CN114938951B (en) | 2025-07-25 |
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210610264.3AActiveCN114938951B (en) | 2022-05-31 | 2022-05-31 | Respiration data processing method |
| Country | Link |
|---|---|
| CN (1) | CN114938951B (en) |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116474220A (en)* | 2022-12-30 | 2023-07-25 | 深圳融昕医疗科技有限公司 | Method for detecting ineffective ventilation event of breathing machine and breathing machine equipment |
| WO2024067030A1 (en)* | 2022-09-29 | 2024-04-04 | 荣耀终端有限公司 | Signal denoising method and electronic device |
| CN118370531A (en)* | 2024-06-19 | 2024-07-23 | 宁波星巡智能科技有限公司 | Multi-mode human respiration monitoring method, device, equipment and medium |
| EP4427672A1 (en)* | 2023-03-10 | 2024-09-11 | Lynx Health Science GmbH | Method, device and system for detecting distinct repetitive events from signals |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6099481A (en)* | 1997-11-03 | 2000-08-08 | Ntc Technology, Inc. | Respiratory profile parameter determination method and apparatus |
| JP2003000561A (en)* | 2001-06-18 | 2003-01-07 | Canon Inc | R wave recognition method, RR interval measuring method, heart rate measuring method, RR interval measuring device, and heart rate measuring device |
| US20050119586A1 (en)* | 2003-04-10 | 2005-06-02 | Vivometrics, Inc. | Systems and methods for respiratory event detection |
| US20070135726A1 (en)* | 2005-12-08 | 2007-06-14 | Shenzhen Mindray Bio-Medical Electronics Co., Ltd. | Method for improving recognition rate of respiratory wave |
| US20080082018A1 (en)* | 2003-04-10 | 2008-04-03 | Sackner Marvin A | Systems and methods for respiratory event detection |
| AU2011203044A1 (en)* | 2003-04-10 | 2011-07-14 | Adidas Ag | Systems and methods for respiratory event detection |
| US20150059755A1 (en)* | 2012-04-13 | 2015-03-05 | Resmed Limited | Apparatus and methods for ventilatory treatment |
| JP2016116746A (en)* | 2014-12-22 | 2016-06-30 | 国立大学法人三重大学 | Respiration stability evaluation method, and device for the same |
| WO2017118127A1 (en)* | 2016-01-05 | 2017-07-13 | 深圳和而泰智能控制股份有限公司 | Heartbeat signal processing method, device and system |
| US20170209074A1 (en)* | 2014-07-28 | 2017-07-27 | S & V Siu Associates, Llc | Method and apparatus for assessing respiratory distress |
| KR101800739B1 (en)* | 2016-10-13 | 2017-11-24 | 한국한의학연구원 | Apparatus and method for detecting respiration rate |
| US20190117097A1 (en)* | 2017-10-19 | 2019-04-25 | Hill-Rom Services Pte. Ltd. | Respiration rate estimation from a photoplethysmography signal |
| CN114266812A (en)* | 2021-12-10 | 2022-04-01 | 南京佗道医疗科技有限公司 | A Phase Registration Method of Lesion Feature Motion Curve and Respiratory Fluctuation Curve |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6099481A (en)* | 1997-11-03 | 2000-08-08 | Ntc Technology, Inc. | Respiratory profile parameter determination method and apparatus |
| JP2003000561A (en)* | 2001-06-18 | 2003-01-07 | Canon Inc | R wave recognition method, RR interval measuring method, heart rate measuring method, RR interval measuring device, and heart rate measuring device |
| US20050119586A1 (en)* | 2003-04-10 | 2005-06-02 | Vivometrics, Inc. | Systems and methods for respiratory event detection |
| US20080082018A1 (en)* | 2003-04-10 | 2008-04-03 | Sackner Marvin A | Systems and methods for respiratory event detection |
| AU2011203044A1 (en)* | 2003-04-10 | 2011-07-14 | Adidas Ag | Systems and methods for respiratory event detection |
| US20070135726A1 (en)* | 2005-12-08 | 2007-06-14 | Shenzhen Mindray Bio-Medical Electronics Co., Ltd. | Method for improving recognition rate of respiratory wave |
| US20150059755A1 (en)* | 2012-04-13 | 2015-03-05 | Resmed Limited | Apparatus and methods for ventilatory treatment |
| US20170209074A1 (en)* | 2014-07-28 | 2017-07-27 | S & V Siu Associates, Llc | Method and apparatus for assessing respiratory distress |
| JP2016116746A (en)* | 2014-12-22 | 2016-06-30 | 国立大学法人三重大学 | Respiration stability evaluation method, and device for the same |
| WO2017118127A1 (en)* | 2016-01-05 | 2017-07-13 | 深圳和而泰智能控制股份有限公司 | Heartbeat signal processing method, device and system |
| KR101800739B1 (en)* | 2016-10-13 | 2017-11-24 | 한국한의학연구원 | Apparatus and method for detecting respiration rate |
| US20190117097A1 (en)* | 2017-10-19 | 2019-04-25 | Hill-Rom Services Pte. Ltd. | Respiration rate estimation from a photoplethysmography signal |
| CN114266812A (en)* | 2021-12-10 | 2022-04-01 | 南京佗道医疗科技有限公司 | A Phase Registration Method of Lesion Feature Motion Curve and Respiratory Fluctuation Curve |
| Title |
|---|
| JOHN A. BERKEBILE 等: "Towards Estimation of Tidal Volume and Res- piratory Timings via Wearable-Patch-Based Im- pedance Pneumography in Ambulatory Settings", EMB, 31 December 2021 (2021-12-31)* |
| 李德旺;仇原鹰;盛英;: "呼吸力学参数测量方法的研究", 生物医学工程学杂志, no. 04, 30 August 2006 (2006-08-30)* |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2024067030A1 (en)* | 2022-09-29 | 2024-04-04 | 荣耀终端有限公司 | Signal denoising method and electronic device |
| CN116474220A (en)* | 2022-12-30 | 2023-07-25 | 深圳融昕医疗科技有限公司 | Method for detecting ineffective ventilation event of breathing machine and breathing machine equipment |
| EP4427672A1 (en)* | 2023-03-10 | 2024-09-11 | Lynx Health Science GmbH | Method, device and system for detecting distinct repetitive events from signals |
| WO2024188827A1 (en)* | 2023-03-10 | 2024-09-19 | Lynx Health Science Gmbh | Method, device and system for detecting distinct repetitive events from signals |
| CN118370531A (en)* | 2024-06-19 | 2024-07-23 | 宁波星巡智能科技有限公司 | Multi-mode human respiration monitoring method, device, equipment and medium |
| Publication number | Publication date |
|---|---|
| CN114938951B (en) | 2025-07-25 |
| Publication | Publication Date | Title |
|---|---|---|
| CN114938951A (en) | Respiration data processing method | |
| CN108388912B (en) | Sleep staging method based on multi-sensor feature optimization algorithm | |
| Li | Wavelets for electrocardiogram: overview and taxonomy | |
| CN101815465B (en) | Apnea/hypopnea index derived from ECG | |
| Casaseca-de-la-Higuera et al. | Weaning from mechanical ventilation: a retrospective analysis leading to a multimodal perspective | |
| CN113558584B (en) | Pulse wave preprocessing method based on signal quality evaluation | |
| Janbakhshi et al. | ECG-derived respiration estimation from single-lead ECG using gaussian process and phase space reconstruction methods | |
| WO2006054306A2 (en) | Sleep staging based on cardio-respiratory signals | |
| Houssein et al. | Estimation of respiratory variables from thoracoabdominal breathing distance: a review of different techniques and calibration methods | |
| CN111563451A (en) | Identification method of ineffective inspiratory effort in mechanical ventilation based on multi-scale wavelet features | |
| CN112363139A (en) | Human body breathing time length detection method and device based on amplitude characteristics and storage medium | |
| EP3612087A1 (en) | Methods and system for detecting inhalations and extracting measures of neural respiratory drive from an emg signal | |
| Boucheham et al. | Piecewise linear correction of ECG baseline wander: a curve simplification approach | |
| Korten et al. | Respiratory waveform pattern recognition using digital techniques | |
| CN114266812B (en) | A phase registration method for lesion characteristic motion curve and respiratory fluctuation curve | |
| CN114668384A (en) | An improved method and device for estimating quasi-static lung compliance under pressure-controlled mechanical ventilation | |
| CN116108345B (en) | A second heart sound width split detection method based on parameter estimation | |
| CN113520368A (en) | Cough monitoring method and system and storage device | |
| Rolink et al. | Improving sleep/wake classification with recurrence quantification analysis features | |
| Kumar et al. | Robust multiresolution wavelet analysis and window search based approach for electrocardiogram features delineation | |
| CN117281490A (en) | Method and system for telemetering vital signs of non-anesthetized animals | |
| Yang et al. | Algorithmic detection of sleep-disordered breathing using respiratory signals: a systematic review | |
| Koley et al. | Classification of sleep apnea using cross wavelet transform | |
| Waser et al. | A blind source-based method for automated artifact-correction in standard sleep EEG | |
| Avcı et al. | Comparison of the ANN based classification accuracy for real time sleep apnea detection methods |
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| CB02 | Change of applicant information | ||
| CB02 | Change of applicant information | Address after:210000 building 3, No. 34, Dazhou Road, Yuhuatai District, Nanjing, Jiangsu Province Applicant after:Tuodao Medical Technology Co.,Ltd. Address before:210000 building 3, No. 34, Dazhou Road, Yuhuatai District, Nanjing, Jiangsu Province Applicant before:Nanjing Tuodao Medical Technology Co.,Ltd. | |
| GR01 | Patent grant | ||
| GR01 | Patent grant |