机械通气模型参数辨识与仿真

张玉欣

(北华大学电气与信息工程学院,吉林 吉林 132021)

呼吸系统力学参数是评估人体呼吸功能的重要指标,也是呼吸机参数个性化、精准设置的重要依据,可为医生研究病理、判断病人病情、分析治疗效果、临床诊断等提供一定的理论支持.呼吸系统参数值因人而异,并且随着患者病程发展而变化.因此,研究呼吸方程,实现对病人进行机械通气过程呼吸系统参数实时监测与辨识具有重要意义[1-3].人类呼吸系统的力学结构通常被等效为电学形式.比如,PRC模型中,呼吸肌形成的压力等效为电压源,气道阻力等效为电阻,顺应性等效为电容[4];RIC模型中,气道阻力等效为电阻,顺应性等效为电容,惯性等效为电感[5].机械通气系统也被等效为气动系统,建立的呼吸方程包括气流方程、压力方程和体积方程[6-7].在以上定义的呼吸方程中,肺弹性系数和气道阻力等呼吸力学参数均为常数,无法描述其非线性特性.KANAE等[8]、NAKAMICHI等[9]提出了二阶非线性常微分方程的呼吸模型,肺弹性系数和气道阻力用高阶多项式、模糊逻辑等非线性函数系统描述,给出了描述非线性动态呼吸过程的基本框架.

医用呼吸机是一种流体机械,属于管道流体力学的研究范畴,吉志丽等[10]、温志芳[11]在基于流体力学理论建立的时间切换机械通气过程数学模型中,设定呼吸机参数和患者呼吸力学参数进行仿真试验,流速、流率、肺内压、气道压、潮气量的试验结果与临床测量数据基本一致,可用于呼吸系统模型的参数辨识.本文利用流率、气道压、潮气量的流体力学模型采样数据作为数据样本集,利用数值积分算法和最小二乘法对二阶呼吸系统方程进行参数辨识;比较复合梯形递推法、辛普森递推法、科茨递推法、龙贝格递推法等数值积分法的辨识结果,确定最优数值积分方法.

1 呼吸方程

KANAE等[8]、NAKAMICHI等[9]提出的无自主呼吸病人机械通气时二阶非线性常微分呼吸方程为

(1)

式中:Paw(t)、V(t)、F(t)分别为呼吸机管道末端测得的气道压力、潮气量、气体流量;fE(V)为肺弹性系数;fE(V)V(t)为肺内压(图中用PA表示);Peea为呼吸末端肺泡压力;e(t)为误差值;其余参数为常系数.

式(1)的矩阵形式:

Paw(t)=φT(t)θ+e(t),

(2)

2 呼吸方程的参数辨识

数值积分运算递推公式如下:

(3)

式中:n=0,1,…,N-17;l=1,2,3.当m=0,k=0,1,2,3,4时为梯形序列;当m=1,k=0,1,2,3时为辛普森序列;当m=2,k=0,1,2时为科茨序列;当m=3,k=0,1时为龙贝格序列.

呼吸方程是连续时间函数模型,首先对每个呼吸周期以等间隔T进行采样,得到离散数据的呼吸样本集F(n)、V(n)、Paw(n),其中,n=0,1,2,…,N-1.

呼吸方程的数值积分方程形式为

Paw(n)=φT(n)θ+e(n),

(4)

式中:φ(n)为φ(t)的数值积分结果;e(n)为e(t)的数值积分结果.

记p=(Paw(0),Paw(1),…,Paw(N-1))T,Φ=(φ(0),φ(1),…,φ(N-1))T,e=(e(0),e(1),…,e(N-1))T,式(4)可统一写为

p=Φθ+e.

(5)

(6)

本文比较梯形、辛普森、科茨和龙贝格4种递推数值积分算法对呼吸方程的参数辨识结果.定义积分区间二分次数为4,即积分长度为16.

3 试验与分析

3.1 仿真试验

为了避免直接采集肺内静态压力给病人造成不适,本文肺内压、流率、潮气量、气道压等试验数据均通过软件仿真的方法获得[12-13].具体做法:首先利用MATLAB软件得出时间切换机械通气模型流率、潮气量、肺内压、气道压函数曲线,均为周期函数,图1仅为4个周期的流率、潮气量、肺内压、气道压;其次对流率、潮气量、气道压函数曲线前3个周期以0.05 s进行等间隔采样,共计采样得到60组试验数据.

图1 基于流体力学理论的机械通气模型仿真

试验步骤:

第1步 将60组采样数据作为式(3)4种数值积分递推公式的输入数据.

第3步 根据4组参数辨识结果拟合出4条肺内压曲线.

第4步 分别计算4条拟合肺内压曲线与图1中肺内压曲线的平均相对误差.

第5步 比较4个平均相对误差大小,最小值对应的数值积分算法为最佳算法.

表1 fE(V)=k1+k2V(t)时4种数值积分算法参数辨识结果

将表1中4种数值积分算法辨识出的参数分别代入fE(V)V(t),拟合出1个周期的肺内压曲线,见图2.

图2 肺内压力拟合曲线(fE(V)=k1+k2V(t))

图2中4条肺内压拟合曲线与图1基于流体力学模型肺内压曲线的平均相对误差依次为1.887 5、0.813 3、0.611 3、0.587 7.

图3 肺弹性系数拟合曲线(fE(V)=k1+k2V(t))

表2 fE(V)=k1+k2V(t)+k3V2(t)时4种数值积分算法参数辨识结果

将表2中4种数值积分算法辨识出的参数分别代入fE(V)V(t),拟合出1个周期肺内压的曲线,见图4.

图4 肺内压拟合曲线(fE(V)=k1+k2V(t)+k3V2(t))

图5 肺弹性系数为二阶时肺弹性系数拟合曲线 (fE(V)=k1+k2V(t)+k3V2(t))

3.2 结果分析

本文定义的肺弹性系数为多项式,由图3、图5肺弹性系数拟合曲线可以看出肺弹性系数随潮气量的变化情况,而基于流体力学理论的机械通气模型中肺弹性系数为常数,本文方法对肺弹性系数以及肺内压的辨识更精准,可以为医生的治疗、诊断提供参考.

4 小 结

本文定义的呼吸方程中的肺弹性系数为多项式,通过数值积分法和最小二乘法实现了呼吸方程的参数辨识.仿真试验结果表明,龙贝格递推数值积分法的辨识效果最好.本文方法避免了测量肺内静态压力给病人造成的不适感,可在无实测呼吸数据的情况下,可以将基于流体力学理论的机械通气模型肺内压仿真采样数据作为输入数据,辨识呼吸方程各呼吸特性参数值,以及病人的肺内压力值.本文方法可为医生设置呼吸机参数提供参考依据,同时也可为医生进行病理、生理研究以及临床诊断等提供理论支持.但本文方法并没有考虑同一个病人不同周期呼吸数据的差异,因此,在今后的研究中,将重点解决这个问题,实现个性化、差异化、高准确率的呼吸特性参数辨识.