小波分解结合自适应神经模糊推理系统的呼吸预测研究

朱 丹，王 伟，付东山

（天津医科大学肿瘤医院放疗科，国家肿瘤临床医学研究中心，天津市“肿瘤防治”重点实验室，天津市恶性肿瘤临床医学研究中心，天津 300060）

高精度、高剂量、高疗效、低损失始终是当前肿瘤精确放射治疗所追求的目标。在实时图像引导呼吸运动跟踪或者门控放射治疗中，从图像采集、信号处理到射束调整需要一定的时间，使治疗射束照射滞后于肿瘤运动，产生系统时间延迟，影响胸腹部肿瘤放疗效果，增加正常组织的损伤。呼吸运动预测是实时跟踪放疗中的一项关键技术。

国内外学者进行了大量研究和探索，已经提出了多种呼吸预测算法。Murphy等[1]探索了神经网络(Neural Network,NN)预测呼吸运动的能力，研究发现神经网络算法的精确性优于线性算法。但神经网络算法存在一些不足[2-4]，其容易发生局部最小化问题，影响预测准确度；也可能因为输入层和隐藏层的节点过多产生过拟合问题。Sheng等[5]研究出了一种混合预测算法，应用于CyberKnife®放射外科系统 (Accuray Incorporated,Sunnyvale,CA,USA)的Synchrony®呼吸同步追踪系统中，该混合算法结合了线性算法、模式匹配和模糊逻辑3种方法。线性算法采用最小二乘实现线性预测；模式匹配是在历史呼吸数据中寻找和当前呼吸周期最相似的一段呼吸信号作为参考，以此来预测呼吸运动；模糊逻辑是将历史呼吸数据按位置信息分成大、中、小3部分，每一部分再按速度信息分为正速度、零速度和负速度3部分，构成3×3模糊控制规则表进行呼吸预测。Synchrony®呼吸同步追踪系统已广泛应用于临床实践，具有较好的实时跟踪能力。Kakar等[6]采用自适应神经模糊推理系统（Adaptive Neuro Fuzzy Inference System,ANFIS）预测呼吸信号，该方法结合了神经网络的学习能力和模糊逻辑的推理能力，其基于临床病例数据的实验结果表明，ANFIS的均方根误差在亚毫米级别。史少华等[7]提出了一种基于分离有限状态模型的呼吸预测算法，该算法将呼吸基线与呼吸起伏分离，并分别进行预测，预测误差不会因为呼吸运动基线的明显漂移而受到影响。

在放疗过程中，患者呼吸运动的基线、幅度、频率经常发生变化，增加了呼吸运动的复杂度及其预测的难度。对非规则呼吸的预测仍是跟踪或门控放疗研究的重点和难点。呼吸运动的基线变化会导致患者解剖结构发生显著的时空变化，引起胸腹部肿瘤放疗中几何和剂量学的误差[8]。本文提出一种基于小波分解和自适应神经模糊推理系统的呼吸运动预测算法(WANFIS)。该WANFIS算法利用小波分解将呼吸信号分成基线、低频、高频3部分，并分别采用3种不同方法对3部分信号进行预测。基线部分较平滑，采用线性拟合方法进行预测；低频部分是呼吸信号去除基线后的呼吸幅度部分，采用ANFIS进行预测，ANFIS模型采用位置信息作为输入参数，通过位置的模糊划分构造模糊集合，建立模糊规则控制表；高频部分主要是噪声，采用简单移动平均法进行预测。最后综合3部分预测值作为呼吸运动的预测结果。通过30例非规则呼吸运动临床数据的回顾性分析，比较 NN、Synchrony、ANFIS、WANFIS 4种算法的归一化均方根误差（Normalized Root Mean Square Error，nRMSE）、平均值（Mean）和标准偏差（STDEV）。

1 资料与方法

1.1 呼吸样本数据来源及预测过程 本文所采用的病例数据来源于天津医科大学肿瘤医院射波刀中心CyberKnife治疗胸腹部肿瘤患者光学传感器记录的体表呼吸运动数据，采样频率为26 Hz。选择30例有基线漂移或者振幅突变的非规则呼吸运动数据对照比较NN、Synchrony、ANFIS以及WANFIS 4种呼吸预测算法的精确性。

根据已报道的Synchrony呼吸同步追踪系统的延迟时间[9]，将预测步长τs定为3个采样点。实时预测时，先对N个历史呼吸数据{St-N+1,…,St-1,St}进行归一化预处理，再应用预测算法进行预测，最后对预测值进行反归一化得到τs个采样点后的呼吸运动预测位置St+τs。每新采集一个呼吸数据，按照先进先出原则更新历史呼吸数据。

1.2 小波分解结合自适应神经模糊推理系统的预测方法（WANFIS） WANFIS算法流程图如图1所示。该算法首先通过小波分解将历史呼吸信号分解成基线低频高频部分，再采用不同方法分别进行预测。基线部分反映呼吸信号的基线漂移情况，比较平滑，采用线性拟合算法进行预测；低频部分是去除基线后的呼吸幅度部分，采用ANFIS算法进行预测；高频部分主要是噪声干扰，采用简单移动平均法进行预测。最后综合3部分预测值作为呼吸运动的预测结果。

图1 WANFIS算法流程图Fig 1 Flow chart of WANFIS algorithm

1.2.1 小波分解 小波变换是一种时间频率局部化分析方法，通过伸缩平移运算对信号逐步进行多尺度细化。目前利用小波函数进行分解和重构的主流方法是Mallat算法[10]，分解示意图如图2所示：

图2 Mallat小波分解示意图Fig 2 Schematic diagram of Mallat wavelet decomposition

当j=0 时，cA0即为原信号 s；当j≥1 时，cAj为第j层分解得到的低频系数。cAj与低通滤波器卷积并下采样可得到第j+1层分解后的近似系数cAj+1：

cAj与高通滤波器卷积并下采样可得到第j+1层分解后的细节系数cDj+1：

式(1)和式(2)中：h（m）和g（m）分别是低通滤波器L_D和高通滤波器H_D的冲击响应序列。每分解一次，k的取值范围减小一半。

下采样过程的作用是使信号长度减半，减少小波变换的复杂度。原信号s的长度选择2的整数次幂，本文取N=1 024。对小波分解得到的高低频系数进行单支重构可以得到原信号在该系数对应尺度下的高低频分量，其长度与原信号一致。原信号经过一次分解得到第1层的低频分量和高频分量，分解出的高频分量不再分解，对低频分量继续分解得到第2层的低频分量和高频分量，依次类推，从而得到信号多个尺度层次上的分量。

频带为0～ω的信号经过一次小波分解成高低频两部分，高频部分频带为ω/2～ω，低频部分频带为0～ω/2[11-12]。根据Nyquist采样定理，本文原信号s的频带为0～13 Hz，经过一次小波分解后低频分量对应的频带为0～7.5 Hz，高频分量对应的频带为7.5～13 Hz。多尺度分解依次类推，Mallat塔式分解示意图如图3所示。

图3 Mallat塔式分解示意图Fig 3 Schematic diagram of Mallat pyramidal decomposition

在小波分解中，存在两个影响因素：一是小波基的选取，本文选取5阶Daubechies小波；二是分解层数的确定，本文首先将原信号分解到第8层，将这层的低频分量作为基线部分，然后将原信号去除基线的其余部分进行一次小波分解，分解得到的低频分量作为低频部分，高频分量作为高频部分。图4和图5是两个呼吸运动信号经过小波分解成了基线、低频和高频3部分的示例。

图4 小波分解呼吸信号示例1Fig 4 First example of wavelet decomposition of respiratory signal

图5 小波分解呼吸信号示例2Fig 5 Second example of wavelet decomposition of respiratory signal

在小波分解过程中，为了消除边界效应，需要选择合适的边界延拓方式[13]。在利用小波分解得到基线部分时，采用0阶平滑延拓；在利用小波分解得到低频部分和高频部分时，采用1阶平滑延拓。

1.2.2 线性拟合和简单移动平均 基线部分采用线性拟合方法进行预测。线性拟合预测是根据信号过去n个已知抽样值序列，通过线性方程，预测τs个采样点以后的信号位置的方法：

式（3）中：系数α通过最小二乘法优化，使拟合误差平方和最小。

简单移动平均法是将信号过去n个已知抽样值序列的平均值作为τs个采样点以后的信号位置预测值的方法。利用该方法预测信号分解得到的高频部分，以减小噪声对呼吸预测的干扰。本文取n=6。

1.2.3 自适应神经模糊推理系统（ANFIS） ANFIS是一种基于Takagi-Sugeno模型的模糊推理系统，它将模糊控制的模糊化、模糊推理和反模糊化三个基本过程全部用神经网络结构实现，利用神经网络的学习机制自动地从输入输出样本数据中生成ifthen模糊规则，构成自适应神经模糊控制器，通过在线学习不断调整模糊推理控制规则，实现系统的自适应能力。

图6是ANFIS的一种模型结构[14]，有3个输入参数，一个输出参数f。输入参数来自呼吸信号低频部分，输出f是对应时间τs个采样点以后的位置预测值ANFIS模型结构中每一层的结点具有相同的函数。

图6 ANFIS模型结构图Fig 6 Model structural diagram of ANFIS

式（4）中：{a，b，c}是隶属函数参数，该层每个结点输出相应模糊集的隶属度。

第2层：这一层每个结点的输出是所有输入参数对模糊集隶属度的代数积，代表一条规则的激励强度wr。模型构成2×2×2模糊集合，共生成8条模糊规则（r=1，2，...，8）。

第3层：将各条规则的激励强度归一化：

第4层：这一层每个结点的输出为每条规则的输出，是每条规则的归一化激励强度与对应函数的乘积（r=1，2，...，8），其中：

……

fr是规则r对应的函数，{br，1，br，2，br，3，br，4} 是规则函数参数（r=1，2，...，8）。

第5层：计算所有规则的总输出，得到模型输出：

依据模型输入输出参数构造训练集：

矩阵共k+1行（k=N-3τs-1），前三列为模型的输入数据，第4列为模型的输出数据。应用训练集训练预测模型，采用梯度下降法调节隶属函数参数，采用最小二乘法调节规则函数参数。进行预测时，将最新数据｝输入已训练好的模型，得到τs个采样点以后的预测值：

2 结果

本文基于MATLAB平台对NN、Synchrony、ANFIS、WANFIS 4种算法进行对照比较。图7表示了呼吸信号的基线出现明显漂移状态下四种算法的预测结果，WANFIS 与 NN、ANFIS、Synchrony相比，预测结果更加平滑，更接近实际呼吸运动。WANFIS减小了基线漂移的影响。图8表示了呼吸信号发生振幅突变情况下4种算法的预测结果。WANFIS的预测结果优于其他3种算法。下面采用nRMSE、Mean、STDEV 3种不同的指标定量比较分析四种算法预测呼吸运动的精确性和鲁棒性。每例呼吸信号的预测误差统计长度都为2 000个采样点。

nRMSE代表了呼吸运动预测曲线与实际呼吸信号间的距离，并且排除了呼吸信号本身幅度大小对结果的影响：

式（9）中：xi是信号x第i个实际值是第i个预测值是x的平均值。

表1列出了30例呼吸信号4种算法的nRMSE。其中，28个呼吸病例的WANFIS的nRMSE小于Synchrony，24个小于ANFIS。

图7 基线漂移呼吸信号4种算法的预测结果Fig 7 Prediction results of four algorithms for baseline drift respiratory signal

图8 振幅突变呼吸信号4种算法的预测结果Fig 8 Prediction results of four algorithms for amplitude jump respiratory signal

表1 4种算法的nRMSE比较Tab 1 nRMSE comparison of four algorithms

表2统计了4种算法预测误差的平均值和标准偏差。WANFIS预测误差的平均值和标准偏差明显小于NN和Synchrony；对于30个病例中的25个，WANFIS小于ANFIS。

表2 四种算法预测误差的平均值和标准偏差比较Tab 2 Comparison of Mean and STDEV of prediction error among four algorithms

3 讨论

本研究提出的WANFIS算法不仅具有神经网络的学习能力和模糊逻辑的推理能力，还结合了小波变换的分解能力，将呼吸信号分解成基线、低频、高频3部分，根据每部分的特征分别进行预测。基于30例临床数据回顾性分析，将WANFIS算法与NN、Synchrony、ANFIS这3种典型预测算法进行对照比较。30例临床数据是有基线漂移或者振幅突变的非规则呼吸运动。本文提出的WANFIS算法的nRMSE平均值为0.09，小于NN的0.17、Synchrony的0.11以及ANFIS的0.11。

Kakar等[6]比较了ANFIS对经过呼吸训练和不经过呼吸训练的呼吸信号的预测结果，发现ANFIS更适合预测规则稳定的呼吸信号。对于不规则呼吸信号，ANFIS预测误差过大。因此本文将小波分解和ANFIS结合起来提出WANFIS算法，将呼吸信号的基线分离出来，从不规则呼吸中提取出相对规则的信号成分采用ANFIS进行预测，减小了基线漂移对预测精度的影响，使预测结果更加准确。史少华等[7]利用求若干个周期呼吸数据平均值的方法提取基线，来解决基线漂移问题。但该方法不能立即检测到基线的变化，无法应用于临床。WANFIS通过小波分解可以即时探测到基线的变化，能够应用于临床实时跟踪中。另外，如果预测结果含有太多高频成分，可能会在放疗呼吸运动补偿中导致机器颤动[5]。WANFIS将呼吸信号高频部分分离出来通过简单移动平均进行预测，相当于对呼吸信号滤波去噪，减小了噪声对呼吸预测的影响，使预测的呼吸曲线更加平滑。这是小波分解的另一个作用。

WANFIS算法的精确性和鲁棒性都优于NN、Synchrony以及ANFIS 3种典型算法，其预测结果更接近实际呼吸运动，能更有效地补偿系统时间延迟。本研究的下一步工作是将WANFIS算法结合临床实践，应用于呼吸实时跟踪放疗和呼吸门控放疗。