一种基于数据建模和模型移植的通用血糖预测方法

一种基于数据建模和模型移植的通用血糖预测方法
【专利摘要】本发明公开了一种基于数据建模和模型移植的通用血糖预测方法；该方法针对人体皮下血糖测量信号、人体输注的胰岛素及摄入的饮食进行分析，提取其与未来血糖的关系，并应用到缺乏数据量的新对象上；分别采用传统最小二乘方法和基于潜变量方法建立有外源输入的自回归血糖预测模型，并进行了对比；最后基于预测性能较好的潜变量方法选取基模型实现模型移植；实际应用时针对新对象无需等待获取充分血糖测量信号后重新建模，而是可以直接将基模型移植到新对象上，根据拟合效果在线调整模型系数实现实时血糖预测，大大简化了建模工作量和复杂度。
【专利说明】一种基于数据建模和模型移植的通用血糖预测方法
【技术领域】
[0001]本发明属于血糖数据分析和预测研究领域，特别是涉及一种基于数据建模和模型移植的通用血糖预测方法。
【背景技术】[0002]人体血糖水平的一个显著特点就是时变性，即时序信号测量值存在显著的自相关关系。此外，外源输入和血糖时序数据之间亦密切相关。这两种相关关系是建立血糖预测模型的基础和关键。
[0003]针对血糖信号、外部输入胰岛素信号和饮食信号进行分析、建模，可以根据历史测量值获取未来血糖变化情况。目前，血糖预测模型的建立多采用基于数据驱动的方法。已存在的预测模型可以划分为线性(以基于最小二乘的自回归模型为其典型代表)和非线性(以神经网络方法为其典型代表)两类。线性模型由于其简单的模型结构和算法，获得了广泛应用。比较成熟的建模方法包括自回归(autoregressive, AR)、激励响应(impulse-response, IR)等。根据是否包括外源输入(胰岛素和饮食输入等)作为建模变量，AR模型又可以区分为自回归(AR)模型和具有外源输入的自回归(AR with exogenousinputs, ARX)模型。然而传统的ARX模型主要缺陷有两点:(a)直接针对测量数据利用最小二乘这种最基本的辨识方法拟合血糖预测关系，无法避免该方法本身所具有的缺陷，不能得到满意的预测精度；(b)针对不同个体的血糖动态性、外源输入不进行预先分析，直接建立不同的预测模型(即个体化模型)用于在线预测，这将耗费大量人力物力。

【发明内容】

[0004]本发明的目的在于针对现有血糖预测方法的不足，提供一种基于数据建模和模型移植的通用血糖预测方法。
[0005]本发明的目的是通过以下技术方案来实现的:一种基于数据建模和模型移植的通用血糖预测方法，该方法包括以下步骤:
[0006]步骤1:建模信号预处理:将以一定采样周期At获得的个体皮下血糖信号、输注该个体的胰岛素信号、该个体摄入的饮食信号分别组合成一维时序数据Gkxl、Ikxi, Mkxi，其中，G、1、M分别是血糖信号、胰岛素信号和饮食信号的测量值，k为以采样周期At获得的采样个数；去除血糖信号、胰岛素信号和饮食信号中的尖峰噪声。
[0007]步骤2:获取预测变量矩阵和响应矩阵:预测区间记为ph，血糖信号和两外源输入信号预测变量的长度分别记为le、1Χ> 1?,则二维预测变量矩阵表示为K十=IGmk，/_，Μ_,]，其中，二维预测变量矩阵行数n=k-ph-l+l，k为以采样周期Λ t获得的采样个数，Ph为预测区间，l=maX {1G, l^d-1, lM+d-l}，而d是两外源输入信号样本输入的延迟，二维预测变量矩阵列数?+h+ln。一维响应矩阵表示为Ynxi= [G1+ph, G1+ph+1, Gk]τO
[0008]步骤3:基于LV的血糖预测建模:通过基于LV的方法进行血糖预测建模。
[0009]步骤4:选取基模型进行模型移植:根据步骤3的预测效果，采用LV方法进行模型移植，该步骤通过以下子步骤来实现:
[0010](4.1)确定模型移植的基模型:先从实际角度考虑选取LV方法中胰岛素前的系数识别为负值，饮食前的系数识别为正值的对象的模型系数，再从预测精度考虑，选取预测RMSE最小的对象的模型系数作为基模型。
[0011](4.2)确定模型移植调节方向判断依据:根据预测值和真实测量值均值、波峰或波谷判断二者的相对大小。
[0012](4.3)确定模型移植的调节步长:权衡步长减小带来的调节时间增长和精度提高之间的矛盾，确定模型移植的调节步长。
[0013](4.4)确定模型移植的调节策略:对比先调节胰岛素前系数，再调节饮食前系数和先调节饮食前系数，再调节胰岛素前系数的拟合效果，确定调节策略为先调节胰岛素前系数，再调节饮食前系数。
[0014]步骤5:根据步骤4确定的基模型、步长和调节策略对任意个体的未来血糖值进行预测。该步骤通过以下子步骤来实现:
[0015](5.1)在线预测时，在采集到任意个体新数据 Gnewk X 1、I newk X 1、^newk X I
后(下标new代
表新样本，一般数据量较缺乏)，对新数据按步骤2所述的预测变量矩阵排列得到测试数据输入矩阵°
[0016](5.2)调用基模型进行在线修正预测:
[0017](5.2.1)调用基模型进行模型移植:
[0018]对于每个新的测试数据，由以下步骤计算直接移植基模型后对测试数据的预测值I
[0019]Y = XnmtlV;Cl)
[0020]其中，Xnew为对新数据按步骤2所述的步骤获得的预测变量矩阵，σ是PLS-CCA方法求出的基模型的回归系数向量。
[0021](5.2.2)根据预测值和测试数据的误差不断在线修正两外源输入系数，直到拟合RMSE最小，得到修正后的模型系数σ，，由以下步骤计算移植基模型在线修正后获得ph个
采样点后的预测值1:
[0022]Y, = XrC1 ;(2)
[0023]其中&为该个体采集的任意数据按步骤2所述的预测变量矩阵排列后的矩阵或向量，σ '为修正后的回归系数向量。
[0024]与现有技术相比，本发明的有益效果是:本发明所提出的通用血糖预测方法实际应用时针对新对象无需等 待获取充分血糖测量信号后重新建模，而是可以直接将基模型移植到新对象上，根据拟合效果在线调整模型系数实现实时血糖预测，大大简化了建模工作量和复杂度，将极大地降低建模花费；且预测模型采用基于潜变量建模的方法，外源输入经过平滑处理，提高了预测精度。本发明易于实施，为血糖预测建模方法的研究指明了新的方向。
【专利附图】

【附图说明】[0025]图1是本发明所述的通用血糖预测模型建模流程图。
[0026]图2是采用相同条件下的训练数据和测试数据时LS方法的预测结果图(分析对象为儿童组第十人的血糖信号);
[0027]图3是采用相同条件下的训练数据和测试数据时LV方法的预测结果图(分析对象为儿童组第十人的血糖信号);
[0028]图4是采用不同条件下的训练数据和测试数据时LS方法的预测结果图(分析对象为儿童组第十人的血糖信号);
[0029]图5是采用不同条件下的训练数据和测试数据时LV方法的预测结果图(分析对象为儿童组第十人的血糖信号);
[0030]图6是只调节胰岛素前系数，步长分别为1，0.1, 0.01, 0.001对13个训练数据三
组对象的拟合精度比较图；
[0031]图7是只调节饮食前系数，步长分别为1，0.1, 0.01, 0.001对13个训练数据三组
对象的拟合精度比较图。
[0032]图8是模型移植训练数据个数分别为13、30、50时对一天数据预测结果比较图(分析对象为青少年组第一人的血糖信号)。
[0033]图9是50个训练数据时基于LV直接建模和模型移植的预测结果比较图(分析对象为青少年组第一人的血糖信号)。
【具体实施方式】
[0034]如图1所示，本发明基于数据建模和模型移植的通用血糖预测方法，该方法包括以下步骤:
[0035]步骤1:建模信号预处理
[0036]针对以一定采样周期At获得的个体皮下血糖信号、输注该个体的胰岛素信号、该个体摄入的饮食信号(这里At=5min)，将其分别组合成一维时序数据Gkxl、Ikxl、Mkxl，其中，G、1、M分别是血糖信号、胰岛素信号和饮食信号的测量值，k为以采样周期At获得的采样个数，去除其中的尖峰噪声。本实例中，我们共有来自于三组对象的采样信号，第I组为青少年组，第2组为成人组，第3组为儿童组，每组10人，三组共30人。每个对象的采样信号包括四天的数据。
[0037]步骤2:获取预测变量矩阵和响应矩阵
[0038]该步主要用于将步骤I中采集的信号进行整合排列，得到训练数据输入输出矩阵，测试数据输入矩阵。预测区间记为ph，血糖信号和两外源输入信号预测变量的长度分别记为?。、ιχ> ιΜ，则二维预测变量矩阵表示为X,,、=[(--4;5/,<风吣]，其中，二维预测变量矩阵行数n=k-ph-l+l，k为以采样周期At获得的采样个数，Ph为预测区间，l=max{lG, l^d-1, lM+d_l}，而d是两外源输入信号样本输入的延迟，二维预测变量矩阵列数Jx=Ig+Ii+Imo 一维响应矩阵表示为 Ynxi= [G1+Ph，G1+Ph+1, GJt0
[0039]具体步骤为用一个长度为ph+le个采样点的一维滑动窗口滑过一维时序血糖信号数据Gkxi，每次移动 一个采样点，共移动η次。将每次滑动窗口中的数据作为一个新的行向量，则可以组合为一个二维数据矩阵，截取该矩阵的前IG列，则为C ;，截取该矩阵的最后一列则为Υη>α。分别用一个长度为I1和Im个米样点的一维滑动窗口从第l(；-d个米样点滑过一维时序胰岛素信号数据Ikxi和一维时序饮食信号数据Mkxi，每次移动一个采样点，共移动η次。分别将每次滑动窗口中的数据作为一个新的行向量，则组合为/,?
[0040]对个体模型的研究表明，当有7-8个历史血糖样本，1-2个外源输入样本，以及5-6个输入样本的延迟会获得比较好的预测模型。下面我们针对7个历史血糖样本，即1g=7 ;胰岛素和饮食输入均为I阶，即Ii=Im=I ；5个输入样本延迟，即d=5进行研究。这样l=max{7, 1+5-1，1+5-1} =7, jx=7+l+l=9。一般预测区间ph=6,即预测30分钟后的血糖值,这对未来血糖预测较为有意义，可以保证一定的预测精度和可靠性。
[0041]所以用一个长度为13个采样点的一维滑动窗口滑过Gkxi，每次移动一个采样点，共移动k-12次。将每次滑动窗口中的数据作为一个新的行向量，则可以组合为一个二维数据矩阵，截取该矩阵的前7列，则为G(k_12)X7，截取该矩阵的最后一列则为Y(k_12)xl。分别用一个长度为I和I个采样点的一维滑动窗口从第2个采样点开始滑过Ikxi和Mkxi，每次移动一个采样点，共移动k-12次。分别将每次滑动窗口中的数据作为一个新的行向量，则组合
为 I (k-12) X I 和 M(k-12) X I。
[0042]步骤3:基于LV的血糖预测建模:
[0043]预测模型建立有两种方法，一种是现有的基于LS (least square，最小二乘)的方法，另一种是基于LV (latent variable,潜变量)的建模方法。具体步骤如下:
[0044](3.1)根据血糖信号和两外源输入信号组成步骤2所述的预测变量矩阵和响应矩阵。
[0045]每个对象采样信号的前三天的数据用于每个对象ARX模型训练，称为训练数据，最后一天的数据用于每个对象A`RX模型性能的测试，称为测试数据。则此时训练数据有864个，测试数据有288个。则此时步骤2中训练数据的k=864，测试数据的k=288。
[0046](3.2)如果通过基于LS的方法进行血糖预测建模，则根据血糖信号和两外源输入信号采用LS方法建立ARX预测模型:
[0047]Ys = XQ +F = Ys+F.,(I)
[0048]其中，Ys为血糖真实值，X为用步骤2获取的预测变量矩阵，Θ为回归系数，F为模型误差，Fi为用LS方法获得的血糖预测值。
[0049]如果通过基于LV的方法进行血糖预测建模，根据血糖信号和两外源输入信号采用LV方法建立ARX预测模型:
[0050]基于LV 的 ARX 预测模型米用 PLS-CCA (partial least square-canonicalcorrelation analysis,偏最小二乘-典型相关性分析)的特征提取方法进行建模。先用PLS提取初步潜变量组，由于PLS方法提取的潜变量LVs不能保证其与响应变量的密切相关关系，再用CCA对其进行后处理，从而提取其中与相应变量密切相关的部分进行回归建模与预测。步骤如下:
[0051 ] (3.2.1)外源输入平滑处理
[0052]为提高预测的精度，分别使用如下两式对胰岛素和饮食两外源输入进行平滑处理:
[0053]H, =...........................r---.....................................1(2)
1?50s2+55i + l
【权利要求】
1.一种基于数据建模和模型移植的通用血糖预测方法，其特征在于，该方法包括以下步骤: 步骤1:建模信号预处理:将以一定采样周期At获得的个体皮下血糖信号、输注该个体的胰岛素信号、该个体摄入的饮食信号分别组合成一维时序数据Gkxl、Ikxl、Mkxl，其中，G、1、M分别是血糖信号、胰岛素信号和饮食信号的测量值，k为以采样周期At获得的采样个数；去除血糖信号、胰岛素信号和饮食信号中的尖峰噪声。 步骤2:获取预测变量矩阵和响应矩阵:预测区间记为ph，血糖信号和两外源输入信号预测变量的长度分别记为le、lx> 1Μ，则二维预测变量矩阵表示为HGfrxkJmttiMlmlJ,其中，二维预测变量矩阵行数n=k-ph-l+l，k为以采样周期At获得的采样个数，ph为预测区间，l=max{le，^+d-l, lM+d_l}，而d是两外源输入信号样本输入的延迟,二维预测变量矩阵列数?+W —维响应矩阵表示为Ynxi= [G1+Ph，G1+Ph+1, GJt0 步骤3:基于LV的血糖预测建模:通过基于LV的方法进行血糖预测建模。 步骤4:选取基模型进行模型移植:根据步骤3的预测效果，采用LV方法进行模型移植，该步骤通过以下子步骤来实现: (4.1)确定模型移植的基模型:先从实际角度考虑选取LV方法中胰岛素前的系数识别为负值，饮食前的系数识别为正值的对象的模型系数，再从预测精度考虑，选取预测RMSE最小的对象的模型系数作为基模型。 (4.2)确定模型移植调节方向判断依据:根据预测值和真实测量值均值、波峰或波谷判断二者的相对大小。 (4.3)确定模型移植的调节步长:权衡步长减小带来的调节时间增长和精度提高之间的矛盾，确定模型移植的调节步长。 (4.4)确定模型移植的调节策略:对比先调节胰岛素前系数，再调节饮食前系数和先调节饮食前系数，再调节胰岛素前系数的拟合效果，确定调节策略为先调节胰岛素前系数，再调节饮食前系数。 步骤5:根据步骤4确定的基模型、步长和调节策略对任意个体的未来血糖值进行预测。该步骤通过以下子步骤来实现:
(5.1)在线预测时，在采集到任意个体新数据 Gnewk X 1、I newk X 1、^newk X I后(下标new代表新样本，一般数据量较缺乏)，对新数据按步骤2所述的预测变量矩阵排列得到测试数据输入^巨夏'hm.nx j' ° (5.2)调用基模型进行在线修正预测: (5.2.1)调用基模型进行模型移植: 对于每个新的测试数据^ ?由以下步骤计算直接移植基模型后对测试数据的预测?丨.值r:? =(O 其中，Xnrat为对新数据按步骤2所述的步骤获得的预测变量矩阵，σ是PLS-CCA方法求出的基模型的回归系数向量。 (5.2.2)根据预测值和测试数据的误差不断在线修正两外源输入系数，直到拟合RMSE最小，得到修正后的模型系数O '，由以下步骤计算移植基模型在线修正后获得Ph个采样点后的预测值1: ^ -(D JL V 一* \β n\ J 其中I为该个体采集的任意数据按步骤2所述的预测变量矩阵排列后的矩阵或向量，σ/为修正后的回归系数向量。
2.根据权利要求1所述的基于数据建模和模型移植的通用血糖预测方法，，其特征在于，所述步骤2具体为:用一个长度为ph+le个采样点的一维滑动窗口滑过一维时序血糖信号数据Gkxi，每次移动一个采样点，共移动η次。将每次滑动窗口中的数据作为一个新的行向量，则可以组合为一个二维数据矩阵，截取该矩阵的前Ie列，则为，截取该矩阵的最后一列则为ΥηΧ1。分别用一个长度为I1和Im个米样点的一维滑动窗口从第l(；-d个米样点滑过一维时序胰岛素信号数据Ikxi和一维时序饮食信号数据Mkxi，每次移动一个采样点，共移动η次。分别将每次滑动窗口中的数据作为一个新的行向量，则组合为和M,,./u.
3.根据权利要求1所述基于数据建模和模型移植的通用血糖预测方法，其特征在于，所述步骤3基于LV的血糖预测建模；具体步骤如下: (3.1)根据血糖信号和两外源输入信号组成步骤2所述的预测变量矩阵和响应矩阵。 (3.2)根据血糖信号和两外源输入信号采用LV方法建立ARX预测模型:基于LV的ARX预测模型采用PLS-CCA的特征提取方法进行建模。先用PLS提取初步潜变量组，由于PLS方法提取的潜变量LVs不能保证其与响应变量的密切相关关系，再用CCA对其进行后处理，从而提取其中与相应变量密切相关的部 分进行回归建模与预测。步骤如下: (3.2.1)外源输入平滑处理 为提高预测的精度，分别使用如下两式对胰岛素和饮食两外源输入进行平滑处理:

I Y Yij/*\ Λ.Λ %a9\ J1 1750/+55s+1 H1 =-^;(4)
2450s2+55s +1 其中H1表示外源胰岛素输入平滑处理使用的传递函数，H2表示外源饮食输入平滑处理使用的传递函数。 (3.2.2)数据预处理 将所有个体的预测变量矩阵和响应矩阵组合到一起，针对组合后的预测变量矩阵和响应矩阵内任意一点的变量xu，对该变量进行减均值、除以标准差的全局标准化处理，标准化处理的计算公式如下:............-* ,χ(5)
c ’ 其中:i是组合后的矩阵任一列的均值，Sj是相应列的标准差。其计算公式如下:
I —U
X ? —- / Λ；
M ft ;[μ一....Sj =j - Xj f IM %(6) 其中:M是组合后任一列样本的总量。 (3.2.3)利用PLS提取潜变量组T:
T = XR ,；(7)
Jf = U.H ； 其中，T是由多个PLS潜变量构成的潜变量矩阵，X为用步骤2获取的预测变量矩阵，R是PLS方法的系数矩阵，W是由PLS求出的权重矩阵，P是对应于T的负载矩阵。 (3.2.4)利用CCA对PLS潜变量进行后处理，得到最终潜变量U: U=TV ；(8) 其中，V是对应于潜变量矩阵T的权重向量。因响应矩阵Y是单变量输出响应，由于CCA方法本身的特点，故最终只需提取一个潜变量。 (3.2.5)建立ARX模型: 用最小二乘法求取潜变量与响应变量间的回归系数q: Q=(UtU)-1UtY;(9) 故最终基于LV方法建立的ARX模型为:^ = (10)
O=RVq 其中，O是通过PLS与CCA两种方法结合获得的预测模型回归系数，fv为基于LV方法获得的血糖预测值。 (3.3)衡量预测精度:将获得的血糖预测结果与实际测量值进行对比。在获得一系列新的预测信号后，可以根据下式中定义的RMSE (均方根误差)性能指标计算针对新样本的血糖预测精度: RMSE = ^Σ(J7(O W) ；(11) 其中，y(i)表示原始的血糖测量值，表示采用不同方法获得的血糖预测结果，N为样本数。
【文档编号】G06F19/00GK103605878SQ201310362553
【公开日】2014年2月26日 申请日期:2013年8月19日 优先权日:2013年8月19日
【发明者】赵春晖, 喻成侠, 李文卿 申请人:浙江大学