CN106093795A - 一种脑起搏器的电量监测方法以及电量监测系统 - Google Patents

Abstract

本发明涉及一种脑起搏器的电量监测方法，该方法包括：获取脑起搏器的刺激参数；通过最小二乘算法功耗预测模型计算该脑起搏器的电池相关参数获得第一组电池功耗参数；通过支持向量机功耗预测模型计算该脑起搏器的电池相关参数获得第二组电池功耗参数；以及将该第一组电池功耗参数和第二组电池功耗参数中的至少一组显示给用户。本发明还涉及一种采用上述电量监测方法的脑起搏器电量监测系统。

Description

一种脑起搏器的电量监测方法以及电量监测系统

技术领域

本发明涉及医疗器械相关技术领域，尤其，涉及一种植入式医疗器械(Implantable Medical Device,IMD)。

背景技术

脑起搏器，又称脑深部电刺激(Deep Brain Stimulation,DBS)，是目前治疗晚期和药物难治性运动障碍和精神障碍疾病的有效技术。特别是在帕金森病的治疗上疗效显著。近期研究成果表明脑起搏器也可用于治疗阿兹海默症以及强迫症、抑郁症等精神类疾病。参见图1，脑深部电刺激的实施器件即为脑起搏器，它是一套可植入式微电子装置，包括脉冲发生器(IPG)、电极和延伸导线。脑起搏器开刺激后，脉冲发生器将根据体外程控仪设定的刺激参数产生持续电脉冲，电流经由延伸导线通过电极触点刺激丘脑底核或苍白球内侧核的神经核团。由于这些部件均需植入体内，因此脑起搏器具有功耗低、体积小、安全性强等特点，同时对其材料特性、温度特性和外形设计都有较高的要求。通过体外程控仪就可以调整刺激幅度、频率和脉宽在内的所有刺激参数。脑起搏器采用的电池包括可充电电池和不可充电电池。

在临床上需要知道脑起搏器植入者的DBS的不可充电电池的电量水平和剩余寿命。不准确的电池寿命预测将会导致过早地进行DBS更换手术，造成电池的浪费和患者财产损失；或者过晚进行更换手术造成意外停机，对治疗效果不利。对于DBS适应症患者而言，术后的刺激参数(幅度、频率、脉宽)调节是保证疗效最为重要的手段，几乎每位DBS植入者都会经历调参的过程，目前的寿命预测方案只适用于典型(常用)刺激参数组合，并且刺激参数的每一次调整都将带来误差的累积，这使得预测结果很不可靠以致无法应用于临床。但从实际的角度看，目前只能采用这种近似估计的方法，因为在不考虑阻抗的前提下，幅度、频率、脉宽三个刺激参数的组合有一百余万种，其中电压模式下就有五十余万种组合，将每种参数组合下的功耗都测出来显然是不切实际的。

因此，需要寻求高精度的预测方法对脑起搏器进行电量预测，并且要求该预测方案可以实现与临床特异性相结合的针对每一种刺激参数组合和患者个体情况的(如阻抗)DBS电量监测和寿命评估，并且能够给出对应的误差范围和分布特点，具备一定的精度和可靠性，最终可应用于临床。

发明内容

本发明提供一种脑起搏器的电量监测方法以及电量监测系统。

一种脑起搏器的电量监测方法，该方法包括：获取脑起搏器的刺激参数；通过最小二乘算法功耗预测模型计算该脑起搏器的电池相关参数获得第一组电池功耗参数；通过支持向量机功耗预测模型计算该脑起搏器的电池相关参数获得第二组电池功耗参数；以及将该第一组电池功耗参数和第二组电池功耗参数中的至少一组显示给用户。

根据上述脑起搏器的电量监测方法，其中，所述刺激参数包括：幅度amp，脉宽pw，阻抗im，以及频率cf；所述通过最小二乘算法功耗预测模型计算该脑起搏器的电池相关参数的方法为：通过公式(1)-(2)计算脑起搏器的功耗Cd：

η(amp)＝-0.0007×amp⁴+0.0178×amp³-0.1590×amp²+0.6172×amp-0.0473 (2)

其中，η(amp)表示任意幅度下的效率。

根据上述脑起搏器的电量监测方法，其中，所述公式(1)中的空载功耗通过公式(3)-(5)计算，所述公式(1)中的负载理论功耗通过公式(6)计算：

Ncd(30)＝-0.0091×amp³+0.6963×amp²+1.0049×amp+9.7751 (3)

Ncd(150)＝0.0903×amp³-0.6880×amp²+15.059×amp+14.909 (4)

Ncd(amp,cf)＝Ncd(30)+(Ncd(150)-Ncd(30))/(150-30)·(cf-30) (5)

负载理论功耗＝1000×n×amp²·pw/(2.8×im/cf) (6)

其中，Ncd(30)和Ncd(150)分别表示频率为30Hz和150Hz的空载功耗，Ncd(amp,cf)表示任意幅度和频率下的空载功耗，公式(6)中当脑起搏器为双通道时n＝2，当脑起搏器为单通道时n＝1。

根据上述脑起搏器的电量监测方法，其中，所述通过最小二乘算法功耗预测模型计算该脑起搏器的电池相关参数的方法进一步包括：根据公式(7)计算脑起搏器电量消耗D；根据公式(8)计算脑起搏器开机前剩余电量F；根据公式(9)计算脑起搏器n次调节参数后的剩余电量Q；以及根据公式(10)计算脑起搏器电池寿命P：

$D=\begin{array}{cc}{\∫}_0^{t}Cd& dt\end{array}---\left(7\right)$

F＝BC-BC×1％×t₁/12-10×24×30t₂/1000 (8)

$Q=F-\underset{i=1}{\overset{n}{\Σ}}{D}_i/1000=BC-BC\×1\%\×t_1/12-10\×24\×30\×t_2/1000-\underset{i=1}{\overset{n}{\Σ}}\left(\begin{array}{cc}{\∫}_0^{t_i}{\mathrm{Cd}}_i& dt\end{array}\right)/1000---\left(9\right)$

$P=\frac{BC}{\begin{array}{cc}{\∫}_0^{t}Cd& dt\end{array}}=\frac{7200\×1000}{Cd\×24\×365}---\left(10\right),$

其中，公式(8)中，BC为电池初始容量，单位：mAh；t1为电池从购买到生产加工结束的时间，单位：月；t2为从脑起搏器出厂到开机前时间，单位：月。

根据上述脑起搏器的电量监测方法，其中，所述支持向量机功耗预测模型通过以下方法建立：步骤S10，确定所选的核函数，进入步骤S11；步骤S11，初始化核函数的参数，进入步骤S12；步骤S12，获取实验样本数据，进入步骤S13；步骤S13，将样本数据用于学习或预测对比，如果是学习，进入步骤S14，如果是预测对比，则进入步骤S18；步骤S14，对样本数据进行预处理，进入步骤S15；步骤S15，判断所选参数特征是否具有代表性，如果是，进入步骤S16，如果否，则返回步骤S14；步骤S16，算法训练建立模型，进入步骤S17；步骤S17，采用网络搜索方法优化参数，进入步骤S18；步骤S18，将预测数据与实验数据进行对比，进入步骤S19；步骤S19，判断误差是否在允许范围内，如果是，进入步骤S20，如果否，则返回步骤S16；以及步骤S20，作为预测模型。

根据上述脑起搏器的电量监测方法，其中，所述步骤S10中，所选的核函数为：径向基核函数K(x_i,x_j)＝exp(-γ||x_i-x_j||²),γ>0；所述步骤S12中，获取实验样本数据的方法为：利用自动测试系统采集三个维度的数据，幅度、频率和脉宽；所述步骤S14中，对样本数据进行预处理的方法为归一化或取序列号。

根据上述脑起搏器的电量监测方法，其中，所述归一化公式(1)为：X₁＝2((X-X_min)/(X_max-X_min))-1(1)，其中，X为原数据；X₁为X经归一化处理后的值；X_max和X_min分别为原数据X所在数据组中的最大值和最小值；所述取序列号的方法为：幅度序列号＝幅度×20；脉宽序列号＝(脉宽-30)/10；频率序列号按照自然数以1为增量进行编号。

根据上述脑起搏器的电量监测方法，其中，所述将该第一组电池功耗参数和第二组电池功耗参数中的至少一组显示给用户的方法包括：比较该第一组电池功耗参数和第二组电池功耗参数的误差；以及将误差较小的一组参数显示给用户。

根据上述脑起搏器的电量监测方法，其中，所述将该第一组电池功耗参数和第二组电池功耗参数中的至少一组显示给用户的方法为将该第一组电池功耗参数和第二组电池功耗参数一起显示给用户。

一种脑起搏器电量监测系统，其包括：一控制模块，所述控制模块控制控制整个脑起搏器电量监测系统的工作；一计算模块，所述计算模块用于通过上述任意一种方法进行计算；一显示模块，所述显示模块用于将所述计算模块计算得到的结果显示给用户；一数据输入模块，所述数据输入模块用于输入脑起搏器的数据信息；以及一存储模块，所述存储模块用于存储数据信息。

根据上述脑起搏器电量监测系统，进一步包括一数据获取模块和一通讯模块；所述数据获取模块用于通过该通讯模块与所述脑起搏器的体外程控仪进行通讯，从而获取脑起搏器的刺激参数。

根据上述脑起搏器电量监测系统，进一步包括一比较模块和一提示模块；所述比较模块用于将所述计算模块计算得到的结果与一安全阈值进行比较，从而在危险时通过该提示模块进行提示；且所述示模块通过一图形用户界面将所述计算模块计算得到的结果显示给用户。

相较于现有技术，本发明提供的脑起搏器的电量监测方法以及电量监测系统可以实现对对脑起搏器的电量进行高精度的预测，并且该预测方案可以实现与临床特异性相结合的针对每一种刺激参数组合和患者个体情况的DBS电量监测和寿命评估，具备一定的可靠性，最终可应用于临床。

附图说明

图1为本发明实施例采用的脑起搏器的工作原理图。

图2为本发明实施例采用的脑起搏器与功耗相关的简易电路原理图。

图3为本发明实施例采用的脑起搏器的效率-幅度关系曲线图。

图4为本发明实施例采用的脑起搏器不同幅度和供电电压下的效率误差图。

图5为本发明实施例采用的脑起搏器幅度4.5V不同频率、脉宽和阻抗下的效率误差。

图6为本发明实施例采用的脑起搏器幅度4.5V不同频率(不含2Hz)、脉宽和阻抗下的效率误差。

图7为本发明实施例采用的脑起搏器效率误差随频率和脉宽的分布等值线图。

图8为本发明实施例采用的脑起搏器效率误差随频率和阻抗的分布等值线图。

图9为本发明实施例采用的脑起搏器效率误差随阻抗和脉宽的分布等值线图。

图10为本发明实施例采用的脑起搏器幅度0.9V不同频率、脉宽和阻抗下的效率误差。

图11为本发明实施例采用的脑起搏器幅度0.9V不同频率(不含2Hz)、脉宽和阻抗下的效率误差。

图12为本发明实施例采用的脑起搏器空载功耗误差随频率变化趋势。

图13为本发明实施例采用的脑起搏器空载功耗误差整体趋势。

图14为本发明模拟脑起搏器电量实际监测过程(电池容量变化)。

图15为本发明模拟脑起搏器电量实际监测过程(电池容量百分比变化)。

图16为结构风险最小化示意图。

图17为支持向量机的特征空间示意图。

图18为核函数原理图。

图19为支持向量机预测模型建模流程图。

图20为SVM不对训练和测试样本进行预处理的预测绝对误差与曲线拟合预测绝对误差对比图。

图21为SVM不对训练和测试样本进行预处理的预测相对误差与曲线拟合预测相对误差对比图。

图22为SVM对训练和预测样本进行取序列号预处理的预测绝对误差与曲线拟合预测绝对误差对比图。

图23为SVM对训练和预测样本进行取序列号预处理的预测相对误差与曲线拟合预测相对误差对比图。

图24为SVM样本归一化和取序列号处理的预测绝对误差对比图。

图25为SVM样本归一化和取序列号预测相对误差对比图。

图26为SVM参数优选mse等值线图。

图27为本发明提供的脑起搏器电量监测系统的结构示意图。

图28为本发明提供的脑起搏器电量监测系统的图形用户界面示意图。

主要元件符号说明

脑起搏器电量监测系统	10
		控制模块	100
计算模块	101
		数据获取模块	102
通讯模块	103
		显示模块	104
数据输入模块	105
		存储模块	106
比较模块	107
		提示模块	108

如下具体实施例将结合上述附图进一步说明本发明。

具体实施方式

本发明提供了一种脑起搏器的电量监测方法以及采用该方法的脑起搏器电量监测系统。该方法可以预测不可充电电池的电量和可充电电池的电量消耗情况。对不可充电电池尤为重要。

对脑起搏器进行电量监测和寿命评估的前提是知道电池容量的消耗情况，而电池容量的消耗除了与其内部化学反应相关的自放电和脑起搏器在关刺激时由于电路内部的时钟和内存管理仍在运行而产生的静态功耗以外，最重要的一部分就是不同刺激参数(幅度、频率、脉宽)和阻抗影响下的由脉冲发生器的刺激电流产生的功耗。而刺激参数的调整将伴随每个DBS植入者，对其症状改善有着重要意义。由于自放电和静态功耗可以认为是不变的，用于治疗的刺激电流所产生的功耗却受刺激参数和阻抗的影响很大，甚至可以认为幅度、频率、脉宽和阻抗决定了刺激电流所产生功耗的大小。并且在开刺激时这部分功耗就是总的功耗，自放电和静态功耗只有在关刺激时才有意义。因此，构建电量监测和寿命管理系统的基础是对受刺激参数和阻抗影响的刺激电流产生功耗的准确预测，本发明分别采用将传统的最小二乘法与脑起搏器电路原理相结合建立基于电路原理的最小二乘回归和专门适用于多维、小样本、非线性问题的支持向量机预测方法构建功耗预测模型。

本发明提供一种脑起搏器的电量监测方法，该方法包括：获取脑起搏器的刺激参数；通过最小二乘算法功耗预测模型计算该脑起搏器的电池相关参数获得第一组电池功耗参数；通过支持向量机功耗预测模型计算该脑起搏器的电池相关参数获得第二组电池功耗参数；以及将该第一组电池功耗参数和第二组电池功耗参数中的至少一组显示给用户。

本发明首先介绍基于电路原理的最小二乘回归模型。基于电路原理的最小二乘回归模型是根据简易电路原理图，在不同的刺激参数组合和负载下通过分别计算负载理论功耗、效率和空载功耗来获得对应的功耗预测值。其中效率和空载功耗是利用最小二乘法拟合曲线得到的计算公式。

请参见图2，为与功耗相关的脑起搏器简易电路原理图。从能量损耗的角度来讲，假若没有DC-DC转换电路，从电池流出的能量应该和负载端消耗的能量相等。这里负载理论功耗就是在不考虑DC-DC转换电路的情况下负载端消耗的能量，由刺激参数幅度amp、频率cf、脉宽pw、阻抗im(图2中的R)以及电池供电电压U决定。然而，实际上在脑起搏器电路中有一个很重要的部分，即DC-DC转换电路。这一部分电路中包含有储能电容，可以将较低的电池供电电压以倍数的形式抬升到治疗所需的刺激电压(刺激幅度)。这势必会产生能量的损耗，而在图2中从电池流出的能量和负载端消耗的能量都是可以相对准确地计算出来的，但是DC-DC转换电路所消耗的能量是未知的，就如同一个黑盒子一样，并且会受刺激参数设定、人体阻抗、电池供电电压等影响，所以这一部分能量的损耗步不但是未知的而且是不确定的，这给功耗预测和电量监测带来了很大的难题。

基于电路原理的最小二乘回归模型在这个问题上的处理想法是，通过计算DC-DC电路的转换效率来将这一部分的能量损耗考虑进来。效率即为负载理论功耗与负载实际功耗相除。在特定的频率(150Hz)、脉宽(90μs)、阻抗(1000Ω)、供电电压(2.8V)下拟合一条效率-幅度关系曲线，并利用最小二乘法由该曲线获得效率和幅度之间的关系式，作为在不同刺激参数组合、阻抗和供电电压下的效率计算公式。这里相当于认为DC-DC电路的转换效率只与设定的刺激幅度有关，而没有考虑其他刺激参数、阻抗以及电池供电电压的影响。这种做法显然会存在误差，会在后面对误差来源和影响范围进行全面分析。不过这种想法也有其合理性，这就要追溯到在脑起搏器电路中加入DC-DC转换电路的原因。对于大部分的DBS适应症而言，脉冲发生器典型的刺激幅度设定为2-5V，这也是临床常用的刺激幅度，特殊情况下某些治疗需要刺激幅度高达10V。通常脑起搏器电池在其寿命初期可以提供3V左右的电压，而在植入人体后的长期工作中实际提供的电压会小于这个值，一般认为是2.8V左右。这就造成了一种不平衡，刺激所需电压为2-5V，甚至是10V，而电池实际可提供电压却只有2.8V。为了从较低的电压供应中获得所需的刺激幅度，需要在脑起搏器电路中加入电容堆或电荷泵等结构，即DC-DC转换电路。虽然这些结构可以帮助增大电压，但是所增加的电压的传送却是以产生额外功耗为代价，并且是非线性的。这使得对DC-DC转换电路所消耗能量的预测变得很困难。因DC-DC电路最本质的作用是将较低的电池供电电压抬升到所需的刺激电压幅度，所以刺激幅度对DC-DC转换电路的效率的影响一定是最大的，这就为在除幅度外的其他刺激参数、阻抗和供电电压固定的基础上建立效率与刺激幅度的关系式提供了充分的合理性。然而，这种做法仍然存在误差。首先，这种非线性拟合本身就存在误差。既然刺激幅度是对效率产生影响的最大的因素，那么这个误差也将会是最大的误差来源。其次，从刺激参数的角度，除幅度外的刺激参数对效率的影响都没有考虑，即便通过上面的分析，幅度对效率的影响是最大的，而其他刺激参数对效率的影响可能很小，但是当这些刺激参数组合在一起的时候情况又会有所不同。举个例子，可以在频率为150Hz，脉宽90μs，阻抗1000Ω，供电电压2.8V下绘制一条效率和幅度之间的关系曲线，而当频率、脉宽、阻抗和供电电压改变后又可以再拟合一条效率和幅度关系曲线，这两条曲线之间必然会存在偏差。表面上看，只是改变了对效率影响因素很小的刺激参数(频率和脉宽)，但实际上改变的是刺激参数的组合，而这个组合当中包含了幅度对效率的影响。因此这也是误差来源之一。再次，DBS电池供电电压的变化也将会对效率造成影响。

如上所述，引入DC-DC转换电路的效率是为了解决这部分电路对能量的消耗，相当于效率是用来表征其对能量消耗的一种方式。而DC-DC电路之所以会有能量损耗源于其工作原理或者说它存在的原因，就是要将较低的供电电压抬升到所需的较高的刺激电压。由于随着DBS电池的使用，电池供电电压会逐渐下降。为了达到相同的刺激幅度，供电电压越小就需要DC-DC转换电路将其乘以更大的倍数，也就是将其电压值抬高得更多，那么势必会使DC-DC转换电路产生更大的能量损耗，效率也会随之变化。所以，电池供电电压的变化会对效率产生影响。举个例子，在电池寿命初期，供电电压较大，为达到需要的刺激幅度，经过DC-DC转换电路电压被抬高了2倍；而随着电池的使用，在其寿命末期，供电电压会变小，为了达到同样的刺激幅度，电压则需要被抬高3倍，从而在DC-DC转换电路上产生更大的能量消耗。也就是说，在DBS电池寿命的初期和末期，由于电池供电电压有较大的变化，DC-DC转换电路的效率也会发生改变。简言之，DC-DC转换电路的效率会随电池供电电压的变化而变化。虽然如此，但由电池放电曲线，在电池寿命的很长一段时间内供电电压都是几乎不变的，也就是斜率的绝对值很小几乎为0，只有当供电电压小于2.5V以后，斜率的绝对值会较大，供电电压衰减的会比较快，但这时DBS电池电量也所剩无几，就是说供电电压衰减较快的这一阶段只占了DBS电池寿命中很少的一段时间。综上所述，DBS电池的放电特点决定了供电电压对DC-DC转换电路效率的影响会很小，几乎可以忽略。

空载功耗是指脑起搏器电路上除了负载实际功耗(包含负载理论功耗和DC-DC转换电路的效率)之外的电路部分所消耗的能量，包括时钟电路和内存管理模块的能量损耗。

利用最小二乘法建立的回归模型的功耗计算原理公式(1)如下：

以下将上述公式(1)分解，对每一部分进行分析如下。

1，负载理论功耗公式(2)以及由公式(2)可以得到公式(3)如下：

其中，U为电池电压(供电电压)，单位为V；I为刺激电流(负载消耗/功耗)，单位为μA；amp为刺激幅度，单位为V；pw为刺激脉宽，单位为μs；cf为刺激频率，单位为Hz；im为阻抗(负载)，单位为Ω；T为刺激周期，T＝1/cf，单位为s。公式(2)中乘以2是因为本发明实施例采用双通道脑起搏器为例。如果采用单通道脑起搏器则不需要乘以2。

2，效率公式(4)如下：

其中效率用的是频率150Hz，脉宽90μs，负载1000Ω，2.8V供电电压下采集数据拟合的效率-幅度关系曲线和公式，其中负载理论功耗代入式(3)获得。实测负载功耗和空载功耗则是通过编写计算机能够识别的脚本利用自动测试系统采集的实际数据。最终确定的是效率和幅度之间的定量关系式，即η＝f(amp)。

参见图3，为效率-幅度关系曲线。利用最小二乘法拟合的效率公式(5)为：

η(amp)＝-0.0007×amp⁴+0.0178×amp³-0.1590×amp²+0.6172×amp-0.0473 (5)

(R²＝0.9539)，

其中，amp为刺激幅度，单位为V。

3，空载功耗

频率为30Hz的空载功耗公式(6)，频率为150Hz的空载功耗公式(7)以及任意频率和幅度下的空载功耗公式(8)分别为：

Ncd(30)＝-0.0091×amp³+0.6963×amp²+1.0049×amp+9.7751 (6)

(R²＝0.9776)

Ncd(150)＝0.0903×amp³-0.6880×amp²+15.059×amp+14.909 (7)

(R²＝0.9800)

Ncd(amp,cf)＝Ncd(30)+(Ncd(150)-Ncd(30))/(150-30)·(cf-30) (8)

其中，Ncd(30)和Ncd(150)分别表示频率为30Hz和150Hz的空载功耗；Ncd(amp,cf)表示任意幅度和频率下的空载功耗。

由于空载功耗与频率之间呈线性关系，分别用30Hz和150Hz下的空载功耗-幅度曲线的最小二乘法拟合公式来计算其他任意频率的空载功耗。因此最终确定的空载功耗是刺激幅度和刺激频率的函数，即Ncd＝f(amp,cf)。

根据上面分别对负载理论功耗、DC-DC转换电路的效率和空载功耗的理论分析和拟合公式，可以获得在任意刺激参数组合(刺激幅度、频率、脉宽)、阻抗和供电电压下的功耗计算公式。基于电路原理的最小二乘回归模型最终确定的功耗计算公式(9)如下：

其中，amp：幅度；pw：脉宽；im：阻抗；cf：频率；Ncd：空载功耗；η：效率。

以下对基于电路原理的最小二乘回归模型的功耗预测误差来源和分布特点进行分析。

根据公式(1)，预测任意刺激参数组合下的功耗需要计算三个量：负载理论功耗、效率和空载功耗。由于负载理论功耗是根据脑起搏器电路原理和能量守恒定律计算出来的，不太可能会产生误差。因此对最小二乘预测模型的误差来源和影响范围进行分析主要从两个方面：效率和空载功耗。

1，方法固有误差来源——效率

回顾式公式(4)和公式(5)以及图3，由于效率曲线仅在频率150Hz，脉宽90μs，阻抗1000Ω和电池供电电压2.8V这样一种参数组合下利用最小二乘法绘制效率-幅度变化曲线拟合的效率计算公式，下面分析当刺激参数改变时对效率计算产生的误差。在这里没有考虑到的可能会影响效率计算的可变参数有：频率、脉宽、阻抗和供电电压。另外，虽然拟合曲线已经考虑了幅度变化对效率的贡献，但是由于最小二乘法的拟合方式本身就存在一定的误差，也就是说在上述的特定参数组合下的效率-幅度变化曲线是否能用四阶多项式去拟合，或者说由于这种拟合所产生的效率误差有多大，呈现什么样的分布特点，这一部分会作为幅度对效率误差的影响考虑进来。

(1)分析不同幅度和电池供电电压下的效率误差的分布特点。效率误差指的是在不同幅度和供电电压下利用效率计算公式(5)预测的效率与利用自动测试系统在对应的参数组合下采集实测负载功耗和空载功耗并根据式(4)计算的实际效率之差的绝对值。简言之，这个效率误差指的是在某一参数组合下的模型预测效率和实际效率之间的绝对误差。

参见图4，不同幅度和供电电压下的效率误差是在频率150Hz、脉宽90us，阻抗1K下采集的数据，自变量是幅度和供电电压，因变量是效率误差。由此可以得出两个主要结论：①这6条曲线的整体趋势极为相近，都呈现“U”型。这体现了幅度对效率误差的贡献，不同幅度对效率误差的影响非常明显，幅度在2-7V时效率计算误差很小，小于0.05，而在幅度值的两侧特别是在1V和9.5V附近效率误差较大，最高达到0.2，差不多是幅度在2-7V时效率误差的4倍。②比较不同曲线效率误差趋势之间的差异，即为供电电压对效率误差的影响。可见这个影响很小，可以忽略。只有当电池供电电压<2.2V时，效率误差才会变得越来越大，不过在脑起搏器实际使用时当电池电压<2.2V时基本就认为电量已消耗殆尽不再使用。

(2)分析不同频率、脉宽和阻抗下的效率计算误差的分布特点。既然不同幅度对效率误差影响差异很大，并呈现中间低两端高的“U”型趋势，在这一部分分别取幅度为4.5V和0.9V作为代表分别来对不同频率、脉宽和阻抗下的效率误差进行分析。忽略电池供电电压的影响，均取为2.8V。

①幅度amp＝4.5V时，由式(5)，得到η(4.5)＝0.84533125，这是由最小二乘法拟合公式得到的刺激幅度为4.5V时的效率值，显然这个值是不随频率、脉宽和阻抗的变化而变化的，这与实际情况并不相符，会存在一定误差，下面就来计算这个误差，并对其分布特点和影响范围进行分析。在不同频率、脉宽和阻抗下的实际效率则通过在对应参数组合下利用自动测试系统采集实测负载功耗和空载功耗并根据式(4)进行计算得到。再将计算效率(0.84533125)和实际效率作差取绝对值，从而得到对应参数组合下的效率误差。

参见图5和图6。图5为不同频率、脉宽和阻抗下的效率误差。其中频率取2,30,70,110,130,150,170,200,250(单位Hz)；脉宽取30,70,110,150,200,350,450(单位μs)；阻抗取300,500,1000,2000,2500(单位Ω)。共315个参数组合(包括阻抗)，并且图5和图6是按照频率、脉宽、阻抗从低到高进行绘制的。以阻抗为300Ω为例，如下表1所示：

表1效率误差采集数据分布

从图7-图9的对应不同参数组合下的效率误差分布图中，可以看出在不同的频率、脉宽和阻抗的组合下效率误差的分布情况。从中可以很明显的看到，高频、高脉宽、低阻抗效率误差较大。这意味着，在高频高脉宽低阻抗时利用基于电路原理的最小二乘回归模型预测的功耗值较为不准确，并且误差会很大。这在之后会做详细分析。

②幅度amp＝0.9V时，由式(2-5)，得到η(0.9)＝0.39190693。利用同样的方法分析在幅度为0.9V时不同频率、脉宽和阻抗下的效率误差的大小和分布特点。

参见图10-图11，可以得出以下结论：(1)效率计算误差随刺激参数和阻抗的变化趋势跟4.5V时相近，但误差值整体要更大，差不多是4.5V时的4倍，这与图4是相吻合的；(2)低频(<30Hz)低阻抗(<500Ω)、频率脉宽都较高时效率误差较大；(3)阻抗小于1000Ω时效率误差大部分在0.2-0.4(4.5V时效率误差范围是0.05-0.1，0.9V时是4.5V时的4倍)；(4)阻抗大于1000Ω时效率误差大部分在0.1-0.3(4.5V时效率误差范围是小于0.05)。总结而言，0.9V时效率误差随频率、脉宽和阻抗的变化趋势与4.5V时相近，只是相应误差值抬高了4倍。

2.方法固有误差来源——空载功耗

对于空载功耗，我们认为空载功耗和频率呈线性变化，脉宽对空载功耗没有影响。以下分析不同频率、脉宽和幅度下空载功耗计算误差的分布情况。由图12-13可以得出结论：(1)空载功耗计算误差随刺激频率的增大而增大，在频率为110-150Hz之间较为平稳(在刺激幅度为4.5V，脉宽90μs，供电电压2.8V下绘制)；(2)刺激幅度对空载功耗计算误差的影响较为明显，整体呈现“U”型趋势，这种趋势与刺激幅度对效率预测误差的影响是类似的。当幅度值在2-5V时空载功耗计算误差较小，而在幅度值的两侧误差较大。

当效率和空载功耗分别有Δη和Δb的误差时，会对功耗Cd产生ΔCd的功耗预测误差。以下是效率和空载功耗计算误差对最终预测功耗值误差产生影响的理论推导。

$\left\{\begin{array}{c}\frac{a}{\mathrm{\&eta;}\left(amp\right)}+b=Cd---\left(10\right)\\ \frac{a}{\mathrm{\&eta;}+\mathrm{\&Delta;}\mathrm{\&eta;}}+(b+\mathrm{\&Delta;}b)=Cd+\mathrm{\&Delta;}Cd---\left(11\right)\end{array}\right.$

其中，a为负载理论功耗，b为空载功耗，η为效率，Cd为负载功耗。

负载理论功耗的计算公式(13)为：

$a=\frac{2\&times;1000\&times;{\mathrm{amp}}^2}{im}\&CenterDot;\frac{pw\&CenterDot;cf}{U}---\left(13\right)$

根据对效率误差和空载功耗计算误差在不同参数组合下的分布特点的分析，Δη在低频、低阻抗和高频、高脉宽时较大并且随幅度变化呈现“U”型趋势，即在刺激幅度为2-7V时Δη较小，而在幅度值的两侧Δη较大；Δb随频率的增大而增大，与幅度呈“U”型趋势，即在刺激幅度为2-5V时Δb较小，在幅度值的两侧Δb较大。结合式(12)，Δη项后面乘以的是负载理论功耗a，由式(13)，a与幅度的平方、频率和脉宽成正比，与阻抗成反比。这意味着通过因子a的作用将会放大Δη的影响因素(不同参数组合)对ΔCd的作用。例如：Δη在低阻抗时较大，而阻抗越小a越大，那么ΔCd在低阻抗时将会很大；同样，Δη在高频、高脉宽时较大，而频率和脉宽越大a值也会越大，那么ΔCd在高频、高脉宽时将会很大。而对于空载功耗误差Δb直接会对功耗预测误差ΔCd产生线性作用，也就意味着，Δb的影响因素(不同参数组合)对ΔCd的作用是一样的，并没有被放大或缩小，因此Δb和ΔCd随不同参数组合的分布规律一致。

因此，效率误差和空载功耗计算误差在不同参数组合下的分布特点可以代表最小二乘回归模型最终的功耗预测误差的分布情况。或者说，基于电路原理的最小二乘回归模型功耗预测误差的分布符合这样的特点：只有在典型刺激参数和阻抗范围内，幅度2-5V，频率30-150Hz，脉宽小于300μs，阻抗大于500Ω，可以保持相对较好的功耗预测精度。而在高频高脉宽、低阻抗、幅度值不在2-5V内时功耗预测精度很低。

以下介绍采用最小二乘法构建的最小二乘回归模型的电量状态的长期监测方法。

在某种确定的刺激参数组合下的某段时间内的脑起搏器电池电量消耗就是在这种参数组合下的功耗对时间的积分，即式(14)：

$D=\begin{array}{cc}{\&Integral;}_0^{t}Cd& dt\end{array}---\left(14\right)$

其中，D为脑起搏器在时间t内的电池电量消耗；Cd为功耗，通过功耗预测模型获得。

对于DBS植入者而言，术后刺激参数的调节是改善症状非常重要的一个环节。可以说几乎每个DBS植入者都需要进行刺激参数调节以找到更适合他症状的刺激参数，这个过程很可能要一直伴随着他。而每一次刺激参数的调整，功耗值都将会发生改变，并且随着时间的延长，电量消耗将会越来越大。电量状态的长期监测指的是记录患者每一次调节的刺激参数，利用功耗预测模型估计该刺激参数组合下的功耗，再将这个功耗值对两次调参的时间间隔进行积分，就得到了在一段时间内的DBS电量消耗值，每一次调参都给出预测的DBS剩余电量和该预测方法的可能误差。

在脑起搏器电量的实际监测过程中需要考虑电池从购买到生产加工结束这段期间的自放电现象以及从脑起搏器出厂到开机前数月的电池静态功耗而产生的电量损耗。将这两部分电量损耗考虑进去，在设定刺激参数前，先计算出脑起搏器开机前的总的电量损耗。

根据脑起搏器电池放电特点，自放电每年消耗电池电量的1％，静态功耗为9-10μA，为了保险起见，这里取静态功耗为10μA。脑起搏器开机前剩余电量计算式(15)为：

F＝BC-BC×1％×t₁/12-10×24×30t₂/1000 (15)

其中，BC为电池初始容量，单位：mAh；t1为电池从购买到生产加工结束的时间，单位：月；t2为从脑起搏器出厂到开机前时间，单位：月。

模拟一次调参后，DBS电池剩余电量(单位：mAh)计算式(16)为：

$Q=F-D/1000=BC-BC\&times;1\%\&times;t_1/12-10\&times;24\&times;30\&times;t_2/1000-\begin{array}{cc}{\&Integral;}_0^{t}Cd& dt\end{array}/1000---\left(16\right).$

模拟5次调参后，DBS电池剩余电量(单位：mAh)计算式(17)为：

$Q=F-\underset{i=1}{\overset{5}{\&Sigma;}}{D}_i/1000=BC-BC\&times;1\%\&times;t_1/12-10\&times;24\&times;30\&times;t_2/1000-\underset{i=1}{\overset{5}{\&Sigma;}}\left(\begin{array}{cc}{\&Integral;}_0^{t_i}{\mathrm{Cd}}_i& dt\end{array}\right)/1000---\left(17\right).$

模拟n次调参后，DBS电池剩余电量(单位：mAh)计算式(18)为：

$Q=F-\underset{i=1}{\overset{n}{\&Sigma;}}{D}_i/1000=BC-BC\&times;1\%\&times;t_1/12-10\&times;24\&times;30\&times;t_2/1000-\underset{i=1}{\overset{n}{\&Sigma;}}\left(\begin{array}{cc}{\&Integral;}_0^{t_i}{\mathrm{Cd}}_i& dt\end{array}\right)/1000---\left(18\right).$

参见图14-图15是模拟五次调参过程的电量状态监测情况，两幅图的不同在于纵坐标，也就是电池剩余电量的评估依据不同，图14是电池剩余容量，图15是以百分比的形式来体现电池剩余容量。五次脑起搏器刺激参数调整情况已在图中标出。由于取的刺激参数都是临床常用的典型数值组合，并且功耗值都较小，因此电量剩余状态的预测值和实际值并没有偏差很大。

本发明进一步通过公式(19)计算电池寿命:

$P=\frac{BL}{\begin{array}{cc}{\&Integral;}_0^{t}Cd& dt\end{array}}=\frac{7200\&times;1000}{Cd\&times;24\&times;365}---\left(19\right)$

其中，P为使用寿命(年)；BC为电池总容量(mAh)；Cd为功耗(μA)，通过功耗预测模型获得。为了保险起见，本发明实施例，脑起搏器电池初始容量取为7200mAh(电池生产厂家给出的数据是7500mAh)。

参见如下表2，本发明实施例随机抽取4组刺激参数组合，进行对应的功耗预测和寿命评估。

表2功耗预测和寿命评估情况

接下来介绍本发明的支持向量机(SVM)预测模型。

统计学习理论(Statistic Learning Theory,SLT)是一种专门研究小样本情况下机器学习规律的基本理论和数学架构。由Vapnik提出的基于结构风险最小化的学习机器支持向量机(Support Vector Machine,SVM)作为一种非常有潜力的回归分类技术，是一种基于统计学习理论的模式识别方法。函数回归估计是一种常见的机器学习问题，SVM在这个问题上通过同时控制函数的复杂度和逼近精度，获得了很好的推广能力，用SVM逼近函数的方法即为支持向量回归(Support Vector Regression,SVR)。

SVM方法是建立在统计学习理论的VC维理论和结构风险最小化原理基础上的，根据有限的样本信息在模型的复杂性和学习能力之间寻求最佳折中，是一种专门研究小样本情况下机器学习规律的理论。在这种理论下的统计推理规则不仅考虑了对渐进性能的要求，而且追求在现有有限信息的条件下得到最优结果。SVM理论与神经网络等其他学习方法相比，具有核的参数能够自动地通过优化的方法计算出来，并且避免了局部最小点、过学习等缺陷。总结而言，SVM方法用于回归估计问题不但能够很好地解决以往困扰很多学习方法的小样本、过学习、高维数、局部最优等实际难题，而且有很强的泛化能力。

一方面，基于电路原理的最小二乘回归模型功耗预测的精度有限，并且只有在典型刺激参数组合和阻抗附近有相对于它自己的较好的预测精度。而在非典型参数下功耗预测精度则很低。另一方面，由于在不考虑不同阻抗和电池供电电压影响的前提下，光是在电压模式下刺激参数组合就有五十余万种，如果再将阻抗考虑进来，这个数值将会以倍数的形式增长。而自动测试系统每采集一个参数组合的功耗值需要至少半分钟，将所有刺激参数组合下的功耗都采集出来很难实现。鉴于该问题有刺激幅度、频率、脉宽和阻抗四个变化的刺激参数，属于高维、小样本(最多只能采集2万余种参数组合)、非线性问题，很难直接进行数据拟合。于是尝试采用支持向量机进行预测，并且获得了很好的预测效果，可使功耗预测误差远小于最小二乘回归模型，同时其误差分布更均匀、更平坦，预测功耗的绝对误差和相对误差趋势互补，更有利于进行后期的电量监测和寿命评估，不至因误差的长期积累给寿命预测带来太大的影响。

下面介绍SVM用于回归预测的基本思想和理论基础。

神经网络是基于经验风险最小化原则，但是真正能够直接影响模型实际预测效果的是期望风险或者说实际风险，而结构风险直接决定了期望风险。经验风险最小难以保证实际风险最小，这是神经网络易出现过拟合的原因。在解释过拟合之前，首先要理解训练数据集精度和测试数据集精度是两个不同的概念。所谓过拟合，简单的理解就是，如果一个分类器，其在训练数据集精度越高，在测试数据集精度越低，就说明出现了过拟合。

不同于神经网络，支持向量机SVM是基于结构风险最小化原则，而结构风险是由经验风险和置信范围共同决定的。其中经验风险表征了模型的拟合精度，置信范围则代表模型的推广能力。参见图16，若想期望风险最小，则需要经验风险和置信范围达到一个平衡，也就是在点h*处，可以使结构风险最小，既不会出现欠学习也不会出现过学习现象，使得模型的拟合精度和预测精度都较高。简言之，SVM的基于结构风险最小化原则中包含两个优化指标：经验风险和置信区间。两者共同决定了支持向量的实际风险。对于线性可分数据，SVM的原理是固定经验风险，最优化VC置信范围。这就需要最小化VC维，最终倒向了最基本的问题，就是最大化Margin。从此，该问题变成了纯数学问题。而对于线性不可分数据，则需要加入惩罚因子或者提高维度，将其转变为线性可分问题。相比之下，神经网络等机器学习方法都是固定置信范围，最小化经验风险，因此易出现拟合精度很高而预测精度很低，也就是过拟合的情况。SVM基于结构风险最小化，能够在保证拟合精度的同时最大化置信范围以保证其推广能力。所以，SVM的泛化能力要比神经网络好。

参见图17，支持向量机的基本思想是引入特征空间的概念，通过一个非线性映射把原数域上的非线性问题转化为特征空间上的线性问题，同时还引入核函数的概念，把特征空间上的线性问题又返还到原数域上实施，而不涉及特征空间和非线性映射的具体形式，这样就能在不提高计算复杂度的同时获得最好的推广能力。

为了在高维特征空间中构造最优分类超平面，并不需要以显式形式来考虑特征空间，而只需要能够计算支持向量与特征空间中向量的内积，也就是以核函数的方式来计算。使用核函数带来的好处是，可以将输入空间中的非线性分类面转化为高维的特征空间F中的线性超平面来处理。参见图18，为核函数原理图，体现了将输入空间的非线性分类映射到特征空间的线性分类的过程。

总结而言，利用SVM可以把样本x通过非线性映射ψ(x)映射到高维特征空间F，并在F中求解最优回归函数。简单地说，支持向量机可以实现将低维空间的非线性回归转化为高维空间中的线性回归。而核函数则通过计算支持向量与特征空间中向量内积的方式将高维空间中的线性问题又返还到原输入空间实施。因此，在最优回归函数中采用适当的核函数K(x,xi)代替高维空间中的向量内积ψ(x_i)·ψ(x)，就可以实现某一非线性变换后的线性拟合，而计算复杂度却没有增加。可以说对于支持向量机算法而言，核函数扮演着重要的角色。因此核函数的选择非常重要，将在后面具体分析。

支持向量回归分线性回归和非线性回归。对于线性回归，考虑用线性回归函数(20)

y＝f(x)＝ω·x+b (20)

对于非线性回归，则把样本x通过非线性映射ψ(x)映射到高维特征空间F，并在F中求解最优回归函数，并且采用适当的核函数K(x,xi)代替高维空间中的向量内积ψ(x_i)·ψ(x)，从而实现某一非线性变换后的线性拟合。

f(x,ω)＝(ω·ψ(x))+b (21)

其中，ψ(x)是将样本点映射到高维空间的非线性变换。

SVM回归机可表示为

$min\frac{1}{2}\left|\right|\mathrm{\&omega;}|{|}^2+C\underset{i=1}{\overset{l}{\&Sigma;}}({\mathrm{\&xi;}}_i+{\mathrm{\&xi;}}_i^{*})---\left(22\right)$

其中，ξ_i和为松弛变量，在不能完全满足式(22)的前两个约束条件(不包括ξ_i和)下引入。

满足约束：

$\left\{\begin{array}{c}f(x_i,\mathrm{\&omega;})-y_i\&le;\mathrm{\&epsiv;}+{\mathrm{\&xi;}}_i,\\ y_i-f(x_i,\mathrm{\&omega;})\&le;\mathrm{\&epsiv;}+{\mathrm{\&xi;}}_i^{*},\end{array}\right.---\left(23\right)$

${\mathrm{\&xi;}}_i\&GreaterEqual;0,{\mathrm{\&xi;}}_i^{*}\&GreaterEqual;0,i=1,2,...,l.$

优化函数(22)为二次型，约束条件是线性的，因此是个典型的二次规划问题，可用拉格朗日乘子法求解。引入拉格朗日乘子α_i,η_i,可得

$\begin{array}{c}L(\mathrm{\&omega;},b,{\mathrm{\&xi;}}_i,{\mathrm{\&xi;}}_i^{*},{\mathrm{\&alpha;}}_i,{\mathrm{\&alpha;}}_i^{*},{\mathrm{\&eta;}}_i,{\mathrm{\&eta;}}_i^{*})=\frac{1}{2}\left|\right|\mathrm{\&omega;}|{|}^2+C\underset{i=1}{\overset{l}{\&Sigma;}}({\mathrm{\&xi;}}_i+{\mathrm{\&xi;}}_i^{*})-\\ \underset{i=1}{\overset{l}{\&Sigma;}}{\mathrm{\&alpha;}}_i(\mathrm{\&epsiv;}+{\mathrm{\&xi;}}_i-y_i+\mathrm{\&omega;}\&CenterDot;\mathrm{\&psi;}(x_i)+b)-\underset{i=1}{\overset{l}{\&Sigma;}}{\mathrm{\&alpha;}}_i^{*}(\mathrm{\&epsiv;}+{\mathrm{\&xi;}}_i^{*}+y_i-\mathrm{\&omega;}\&CenterDot;\mathrm{\&psi;}(x_i)-b)-\\ \underset{i=1}{\overset{l}{\&Sigma;}}({\mathrm{\&eta;}}_i{\mathrm{\&xi;}}_i+{\mathrm{\&eta;}}_i^{*}{\mathrm{\&xi;}}_i^{*}).\end{array}---\left(24\right)$

在最优解处有

$\left\{\begin{array}{c}\frac{\&part;L}{\&part;\mathrm{\&omega;}}=\mathrm{\&omega;}-\underset{i=1}{\overset{l}{\&Sigma;}}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})\&CenterDot;\mathrm{\&psi;}\left(x_i\right)=0,\\ \frac{\&part;L}{\&part;b}=\underset{i=1}{\overset{l}{\&Sigma;}}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})=0,\\ \frac{\&part;L}{\&part;{\mathrm{\&xi;}}_i}=C-{\mathrm{\&alpha;}}_i-{\mathrm{\&eta;}}_i=0,\\ \frac{\&part;L}{\&part;{\mathrm{\&xi;}}_i^{*}}=C-{\mathrm{\&alpha;}}_i^{*}-{\mathrm{\&eta;}}_i^{*}=0.\end{array}\right.---\left(25\right)$

将式(25)代入式(24)，将线性可分条件下的原问题变换为其对偶问题，求解如下的极大化问题：

$\begin{array}{c}\max W({\mathrm{\&alpha;}}_i,{\mathrm{\&alpha;}}_i^{*})=\underset{i=1}{\overset{l}{\&Sigma;}}{y}_i({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})-\mathrm{\&epsiv;}\underset{i=1}{\overset{l}{\&Sigma;}}({\mathrm{\&alpha;}}_i+{\mathrm{\&alpha;}}_i^{*})-\\ \frac{1}{2}\underset{i,j=1}{\overset{l}{\&Sigma;}}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})({\mathrm{\&alpha;}}_j-{\mathrm{\&alpha;}}_j^{*})K(x,x_i).\end{array}---\left(26\right)$

求解上述模型的Lagrange对偶问题，得到回归函数：

$f\left(x\right)=\underset{i=1}{\overset{l}{\&Sigma;}}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})K(x,x_i)+b.---\left(27\right)$

其中，K(x_i,x)＝ψ(x_i)·ψ(x)为核函数。

在这个模型中有三个非常重要的参数，分别是惩罚系数C，允许的误差范围ε(在模型训练时ε选为0.01)、核参数γ。惩罚系数C相当于是对误差的宽容度，C越大说明你越不能容忍出现误差；核参数γ隐含地决定了数据映射到新的特征空间的分布。一般地，C越大，拟合误差越小，预测误差越大；ε越大，预测误差越大，拟合误差与ε之间没有明显的单调增大或者减小的关系。对于这三个参数的选择没有明确的规则，一般通过试验给出。由于将允许误差范围ε选定为0.01，还需进行选择的参数有惩罚因子C和核参数γ，关于C和γ的最佳组合的选取会在后面详细介绍。

以下介绍SVM预测模型建模方法流程。SVM建模方法通过Matlab编程实现。参见图19，SVM预测模型建模方法包括以下步骤：

步骤S10，确定所选的核函数，进入步骤S11；

步骤S11，初始化核函数的参数，进入步骤S12；

步骤S12，获取实验样本数据，进入步骤S13；

步骤S13，将样本数据用于学习或预测对比，如果是学习，进入步骤S14，

如果是预测对比，则进入步骤S18；

步骤S14，对样本数据进行预处理，进入步骤S15；

步骤S15，判断所选参数特征是否具有代表性，如果是，进入步骤S16，

如果否，则返回步骤S14；

步骤S16，算法训练建立模型，进入步骤S17；

步骤S17，采用网络搜索方法优化参数，进入步骤S18；

步骤S18，将预测数据与实验数据进行对比，进入步骤S19；

步骤S19，判断误差是否在允许范围内，如果是，进入步骤S20，如果

否，则返回步骤S16；以及

步骤S20，作为预测模型。

上述方法中，主要的三步是步骤S10，核函数的选择；步骤S14，数据预处理；以及步骤S17，参数优化。

所述步骤S10中，在支持向量机中使用的核函数主要有四类：

(1)线性核函数 K(x_i,x_j)＝x_i·x_j

(2)多项式核函数 K(x_i,x_j)＝[γ(x_i·x_j)+coef0]^d,γ>0

(3)径向基核函数(高斯、RBF) K(x_i,x_j)＝exp(-γ||x_i-x_j||²),γ>0

(4)Sigmoid函数 K(x_i,x_j)＝tanh(v(x_i·x_j)+coef0)

对于支持向量机算法而言，核函数的选择非常关键。通常而言，径向基核函数(高斯核函数、RBF核函数)是合理的首选。这个核函数将样本非线性地映射到一个更高维的空间，与线性核不同，它能够处理分类标注(输入特征)和属性(输出功耗值)的非线性关系。并且，线性核是RBF的一个特例，因此，使用一个惩罚因子C的线性核与某些参数(C，γ)的RBF核具有相同的性能。同时，Sigmoid核的表现很像一定参数的RBF核。第二个原因，超参数(hyperparameter)的数量会影响到模型选择的复杂度(因为参数的选择只能靠试验)，而多项式核比RBF核有更多的超参数。最后，RBF核有更少的数值复杂度(numericaldifficulties)，而多项式核涉及高阶运算。此外，Sigmoid核在某些参数下不是合法的(例如：不是两个向量的内积)。当然，也存在一些情形RBF核是不适用的。特别地，当特征维数非常大的时候，很可能只能适用线性核。对于线性核相当于没有将样本映射至高维空间，线性区分(或回归)在原始特征空间中被完成，这是最快的选择。

综合上述分析，由于无论是小样本还是大样本，高维还是低维等情况，RBF核函数均适用，具有较宽的收敛域。并且对于非线性模型来说，RBF核函数是较为理想的回归依据函数，因此选择径向基RBF核函数。

值得注意的是，上面这些核函数中都包含可调节的参数，比如多项式核函数中的参数γ，coef0,d；RBF核函数中的参数γ以及Sigmoid核函数中的参数v，coef0等。根据核函数参数以及式(22)中惩罚因子C和式(23)中允许的误差范围ε选取的不同，模型的预测精度可能会有很大差异。这就是SVM方法的参数选择问题，又称模型选择问题。将在下面进行具体分析。

所述步骤S12中，将需要采集的数据分为三组，共14561个样本：

第一组：共630个。幅度(V)：0-0.9(step：0.1)；频率(Hz)：2,30,70,110,130,150,170,200,250；脉宽(μs):30,70,110,150,200,350,450

第二组：共13175个。幅度(V)：1-4(step:0.1)；频率(Hz):2,10,30,50,70,90,110,130-185(step:5),200-250(step:10)；脉宽(μs)：30-120(step:10),150-450(step:50)。

第三组：共756个。幅度(V)：4.5-10(step:0.5)；频率(Hz)：2,30,70,110,130,150,170,200,250；脉宽(μs)：30,70,110,150,200,350,450。

以上样本的阻抗均为1000Ω。用这14561个样本在不考虑阻抗和电池供电电压的前提下去预测53万余个刺激参数(幅度、频率、脉宽)组合。(仅考虑电压模式)

利用自动测试系统采集三个维度的数据，即幅度、频率和脉宽。由于增加一个维度会给数据采集和预测精度带来影响，因此并没有加入阻抗信息，会在后面作为修正因数考虑。

下面详细解释一下为什么增加一个维度会给数据采集和预测精度带来影响。从数据采集的角度，现在需要采集的样本数量为14561个，如果增加阻抗信息，例如加入300,500,1000,2000,3000(单位Ω)这五个阻抗值，那么意味着将需要采集14561×5＝72805个样本，自动测试系统每采集一个样本需要至少30s的时间，原本只需5天的数据采集时间，加入阻抗信息后至少需要25天。可见增加一个维度给数据采集工作带来了一定的困难。从模型预测精度的角度，对于支持向量机本身而言，输入样本的特征越多，也就是维度越高，SVM预测精度会越低。

所述步骤S14中，样本数据预处理是指对样本数据的输入特征幅度、频率和脉宽的数值进行处理。适用于本问题的样本预处理方式主要有两种：归一化和取序列号。归一化是较为常用的数据预处理方式，而取序列号是笔者根据所研究问题的实际情况提出的。参见图20-图25，通过仿真图形和对比表可以看出，取序列号的处理方案会使模型预测精度大大改善。归一化方案分为两种，一种是对样本的所有输入特征和输出功耗都进行归一化处理；另外一种是只对样本的所有输入特征进行归一化处理，而输出功耗则不进行处理，使用其原始数据。图像仿真说明这两种不同的归一化方式，在模型预测精度上有较大的差异，前者预测精度非常低，后者预测精度比前者高很多，且仅略低于取序列号的方案。下面分别对这三种样本数据预处理方案进行分析，并绘制功耗预测的仿真图形，进行预测误差对比。

所述归一化处理为按照如下归一化方式对数据进行处理：

X₁＝2((X-X_min)/(X_max-X_min))-1 (28)

其中，X为原数据；X₁为X经归一化处理后的值；X_max和X_min分别为原数据X所在数据组中的最大值和最小值。

所述取序列号的方法为：幅度序列号＝幅度×20；脉宽序列号＝(脉宽-30)/10；频率序列号按照从自然数0-61以1为增量进行编号。如下表3所示。

表3对输入样本特征值取序列号预处理

参见表4，为分别采用不对样本预处理、归一化和取序列号的方式模型预测精度对比：

表4

标准偏差

其中，x_i为预测值，x为实际值，n为测试集样本数。

通过对以上3组6幅Matlab仿真图和采用不同样本预处理方案模型预测精度的对比表2，可以得出如下结论：

第一，在支持向量机算法中，如果不对训练和预测样本进行预处理，模型预测精度会非常差，甚至远不如最小二乘回归模型。如图20和图21，以及模型预测精度对比表2中的第二列。

第二，由图22-图25，对样本输入特征(幅度、频率、脉宽)取序列号的方式是最适用于本文所研究问题的样本预处理方案，通过取序列号可使模型预测精度最高。对其原因进行简要分析：①幅度是等比例放大，变为原来的10倍；②频率进行编号处理，降低较大数据值的同时，较小数据并没有变得更小，在保证数据之间关系的同时缩小了数据之间的差距，并且增强了数据之间的内部特征规律和联系，同时使数据范围变小且更有效，数据分布得更合理；③对脉宽做了一个线性变换pw＝30+10×n(pw为实际脉宽值，n为对应的序列号)，由于实际可调脉宽值在30-450us之间，可使数据值有效减小，在减小时也并未打破数据之间的内部特征关系。并且不是直接缩小10倍，而是有一个截距，从而保证了数据的完整性，使得变换后的数据从0开始，并且数值也更小。综上，对样本的输入特征通过取序列号这样的映射关系可以完整的保存不同特征值内部之间的关系，而不会打破特征值内在的关联。另外，根据仿真图形，取序列号和只对样本的所有输入特征进行归一化处理，而输出功耗不进行处理使用其原始数据这两种样本预处理方案的模型预测精度相当，由于篇幅等原因这里没有给出仿真图形对比。

第三，图24和图25中的归一化方案采用的是对样本的所有输入特征和输出功耗都进行归一化处理的方式。可以看出这种样本预处理方案相比于取序列号预测精度略差。对其原因进行简要分析。这种处理方式和只对样本的所有输入特征进行归一化处理，而输出功耗则不进行处理使用其原始数据的方案有一个很大的差别，就是这种方式由于对训练样本的输出功耗进行了归一化，那么SVM经训练学习后得到的功耗预测值相当于也是归一化后，这时就需要对其进行反归一化。即如果先将整个数据集(包括训练数据和预测数据)进行归一化，然后训练SVM模型，再用模型去做预测，再将得到的预测结果进行反归一化，得到预测的功耗数据。在这个过程本身存在一个问题，就是反归一化时由于并不知道预测值对应的功耗的最大和最小值，而是仍然用的实测训练数据集中功耗的最大最小值，除非这两者的最大最小值相同，否则必然会存在偏差。

因此，对于SVM算法的训练，要想得到好的预测效果，最重要的就是数据，即对于原始数据的特征提取和处理。采用不同的样本预处理方案可能会对模型预测精度产生意想不到的较大影响。仿真结果说明，对样本输入特征(幅度、频率、脉宽)取序列号的数据预处理方式可使模型预测精度最高，其预测样本的标准差是最小的。并且尽量不要对样本输出特征即功耗值做任何处理(例如：归一化和反归一化)，这会使预测精度大打折扣。

所述步骤S17中，在所选取的核函数的基础上，对其参数进行优化选择。核函数的选择和参数优选是决定SVM回归预测精度的关键，这往往需要一定的先验知识，目前还没有一般性的结论，这里采用网格搜索的方式对惩罚系数C和核参数γ进行优选。

优化的目标函数是均方误差mse，参加优选的参数是惩罚因子C和核参数γ，优选范围为：log₂C＝-15:15，log₂γ＝-15:15，C和γ的组合共有961种。参数优化方案是利用网格搜索的方式选择使均方误差mse最小的一组参数作为最优的C与γ。在最优参数组合处，SVM的经验风险和置信范围接近最佳的组合(结构风险最小化示意图16中h*处)，既不会出现过学习现象，也不会出现欠学习现象，此时SVM具有很好的推广能力。同时，训练时随机生成训练样本和测试样本，并将3折交叉验证和参数优选相结合，以防止过拟合，提高模型的泛化能力。

图26为某一次参数优选的示意图(刺激幅度0-0.9V，阻抗1000Ω，供电电压2.8V；样本数据预处理方式采用取序列号的方案)。图中横坐标为log₂C，纵坐标为log₂γ，不同颜色的线表示均方误差mse等值线，越接近红色(暖色)表示mse越大，越接近蓝色(冷色)表示mse越小(图中有标出mse的具体数值)。图26中圈出的位置为mse最小处，此时log₂C＝14,log₂γ＝-10，即best C＝16384，bestγ＝9.76563×10^-4。

参见图27，为本发明提供的脑起搏器电量监测系统10，其包括：一控制模块100以及与该控制模块100连接的计算模块101、数据获取模块102、通讯模块103，显示模块104、数据输入模块105以及存储模块106。

所述控制模块100控制控制整个脑起搏器电量监测系统10的工作。所述数据获取模块102用于获取所述脑起搏器的刺激参数，例如：幅度amp，脉宽pw，阻抗im，以及频率cf。所述数据获取模块102可以通过通讯模块103与所述脑起搏器的体外程控仪进行通讯，从而获取脑起搏器的刺激参数。所述数据输入模块105用于输入电池总容量、电池从购买到加工结束的时间以及脑起搏器从出厂到开机的时间等信息。所述计算模块101用于通过上述公式根据刺激参数计算所述脑起搏器的功耗、剩余电量以及使用寿命。所述显示模块104将所述计算模块101计算得到的结果显示给用户。参见图28，所述显示模块104通过一图形用户界面将监控数据显示给用户，其中，显示数据包括：电池总容量、电池从购买到加工结束的时间、脑起搏器从出厂到开机的时间、功耗、剩余电量、误差以及刺激参数幅度amp，脉宽pw，阻抗im，以及频率cf。所述存储模块106用于存储数据信息。

进一步，所述脑起搏器电量监测系统10也可以没有数据获取模块102和通讯模块103，完全通过数据输入模块105来手动输入数据信息。

进一步，所述脑起搏器电量监测系统10也可以通过通讯模块103将所述计算模块101计算得到的结果发送至手机等客户端。

进一步，所述脑起搏器电量监测系统10还可以包括一比较模块107和一提示模块108。所述比较模块107用于将所述计算模块101计算得到的结果与一安全阈值进行比较，从而在危险时通过该提示模块108进行提示。

以上已经给出了本发明的多个实施方式，可以理解的是，在不偏离本公开内容精神以及范围的情况下，可以做出各种变化、替换、以及改变，这些实施方式也在本发明的保护范围内。

Claims (12)

1.一种脑起搏器的电量监测方法，该方法包括：

获取脑起搏器的刺激参数；

通过最小二乘算法功耗预测模型计算该脑起搏器的电池相关参数获得第一组电池功耗参数；

通过支持向量机功耗预测模型计算该脑起搏器的电池相关参数获得第二组电池功耗参数；以及

将该第一组电池功耗参数和第二组电池功耗参数中的至少一组显示给用户。

2.根据权利要求1的脑起搏器的电量监测方法，其特征在于，所述刺激参数包括：幅度amp，脉宽pw，阻抗im，以及频率cf；所述通过最小二乘算法功耗预测模型计算该脑起搏器的电池相关参数的方法为：通过公式(1)-(2)计算脑起搏器的功耗Cd：

η(amp)＝-0.0007×amp⁴+0.0178×amp³-0.1590×amp²+0.6172×amp-0.0473 (2)

其中，η(amp)表示任意幅度下的效率。

3.根据权利要求2的脑起搏器的电量监测方法，其特征在于，所述公式(1)中的空载功耗通过公式(3)-(5)计算，所述公式(1)中的负载理论功耗通过公式(6)计算：

Ncd(30)＝-0.0091×amp³+0.6963×amp²+1.0049×amp+9.7751 (3)

Ncd(150)＝0.0903×amp³-0.6880×amp²+15.059×amp+14.909 (4)

Ncd(amp,cf)＝Ncd(30)+(Ncd(150)-Ncd(30))/(150-30)·(cf-30) (5)

负载理论功耗＝1000×n×amp²·pw/(2.8×im/cf) (6)

其中，Ncd(30)和Ncd(150)分别表示频率为30Hz和150Hz的空载功耗，Ncd(amp,cf)表示任意幅度和频率下的空载功耗，公式(6)中当脑起搏器为双通道时n＝2，当脑起搏器为单通道时n＝1。

4.根据权利要求2的脑起搏器的电量监测方法，其特征在于，所述通过最小二乘算法功耗预测模型计算该脑起搏器的电池相关参数的方法进一步包括：根据公式(7)计算脑起搏器电量消耗D；根据公式(8)计算脑起搏器开机前剩余电量F；根据公式(9)计算脑起搏器n次调节参数后的剩余电量Q；以及根据公式(10)计算脑起搏器电池寿命P：

$\begin{array}{cc}D={\&Integral;}_0^{t}Cd& dt\end{array}---\left(7\right)$

F＝BC-BC×1％×t₁/12-10×24×30t₂/1000 (8)

$Q=F-\underset{i=1}{\overset{n}{\&Sigma;}}{D}_i/1000=BC-BC\&times;1\%\&times;t_1/12-10\&times;24\&times;30\&times;t_2/1000-\underset{i=1}{\overset{n}{\&Sigma;}}\left(\begin{array}{cc}{\&Integral;}_0^{t_i}{\mathrm{Cd}}_i& dt\end{array}\right)/1000---\left(9\right)$

$P=\frac{BC}{\begin{array}{cc}{\&Integral;}_0^{t}Cd& dt\end{array}}=\frac{7200\&times;1000}{Cd\&times;24\&times;365}---\left(10\right),$

其中，公式(8)中，BC为电池初始容量，单位：mAh；t1为电池从购买到生产加工结束的时间，单位：月；t2为从脑起搏器出厂到开机前时间，单位：月。

5.根据权利要求1的脑起搏器的电量监测方法，其特征在于，所述支持向量机功耗预测模型通过以下方法建立：

步骤S10，确定所选的核函数，进入步骤S11；

步骤S11，初始化核函数的参数，进入步骤S12；

步骤S12，获取实验样本数据，进入步骤S13；

步骤S13，将样本数据用于学习或预测对比，如果是学习，进入步骤S14，如果是预测对比，则进入步骤S18；

步骤S14，对样本数据进行预处理，进入步骤S15；

步骤S15，判断所选参数特征是否具有代表性，如果是，进入步骤S16，如果否，则返回步骤S14；

步骤S16，算法训练建立模型，进入步骤S17；

步骤S17，采用网络搜索方法优化参数，进入步骤S18；

步骤S18，将预测数据与实验数据进行对比，进入步骤S19；

步骤S19，判断误差是否在允许范围内，如果是，进入步骤S20，如果否，则返回步骤S16；以及

步骤S20，作为预测模型。

6.根据权利要求5的脑起搏器的电量监测方法，其特征在于，所述步骤S10中，所选的核函数为：径向基核函数K(x_i,x_j)＝exp(-γ||x_i-x_j||²),γ>0；所述步骤S12中，获取实验样本数据的方法为：利用自动测试系统采集三个维度的数据，幅度、频率和脉宽；所述步骤S14中，对样本数据进行预处理的方法为归一化或取序列号。

7.根据权利要求6的脑起搏器的电量监测方法，其特征在于，所述归一化公式(1)为：X₁＝2((X-X_min)/(X_max-X_min))-1(1)，其中，X为原数据；X₁为X经归一化处理后的值；X_max和X_min分别为原数据X所在数据组中的最大值和最小值；所述取序列号的方法为：幅度序列号＝幅度×20；脉宽序列号＝(脉宽-30)/10；频率序列号按照自然数以1为增量进行编号。

8.根据权利要求1的脑起搏器的电量监测方法，其特征在于，所述将该第一组电池功耗参数和第二组电池功耗参数中的至少一组显示给用户的方法包括：比较该第一组电池功耗参数和第二组电池功耗参数的误差；以及将误差较小的一组参数显示给用户。

9.根据权利要求1的脑起搏器的电量监测方法，其特征在于，所述将该第一组电池功耗参数和第二组电池功耗参数中的至少一组显示给用户的方法为将该第一组电池功耗参数和第二组电池功耗参数一起显示给用户。

10.一种脑起搏器电量监测系统，其包括：

一控制模块，所述控制模块控制控制整个脑起搏器电量监测系统的工作；

一计算模块，所述计算模块用于通过上述权利要求1至9所述的任意一种方法进行计算；

一显示模块，所述显示模块用于将所述计算模块计算得到的结果显示给用户；

一数据输入模块，所述数据输入模块用于输入脑起搏器的数据信息；以及

一存储模块，所述存储模块用于存储数据信息。

11.根据权利要求10的脑起搏器电量监测系统，其特征在于，进一步包括一数据获取模块和一通讯模块；所述数据获取模块用于通过该通讯模块与所述脑起搏器的体外程控仪进行通讯，从而获取脑起搏器的刺激参数。

12.根据权利要求10的脑起搏器电量监测系统，其特征在于，进一步包括一比较模块和一提示模块；所述比较模块用于将所述计算模块计算得到的结果与一安全阈值进行比较，从而在危险时通过该提示模块进行提示；且所述示模块通过一图形用户界面将所述计算模块计算得到的结果显示给用户。