非线性科学及其在医学中的应用

摘要

简要介绍了非线性科学中的分形维方法和信息熵方法及其用于房颤信息的处理和分析，计算了未房颤时P波的及房颤时f波的格子维和信息熵，特别是，比较了未房颤时和发生房颤时的相关计算结果的显著差别，从而很容易判断是否发生了房颤。

正文

《1 引言》

1 引言

非线性科学是一门新兴的先进的热门学科, 是诸多非线性理论的总称, 其研究对象涉及复杂大系统。在这类系统中, 整体不能简单地归结为部分之和, 其各个部分之间的相互作用常常是非线性的。

传统的科学方法是通过将复杂的研究对象分解为若干个细小的组成部分进行研究。这种传统方法已广泛用于各个领域 (自然科学, 社会科学, 工程技术, 等等) , 取得了丰硕成果。遗憾的是, 传统方法的固有弱点在于“只见树木, 不见森林”, 缺乏对所研究对象的全面了解和综合分析, 因而用于研究复杂大系统是无法胜任的。

人体正是一个复杂大系统, 其各部分之间的强耦合相互作用是非线性的, 因而非线性科学在医学领域大有用武之地。

笔者介绍分形和分形维、信息熵等非线性理论及其在医学中的应用和笔者的一些有关研究成果。

《2 分形与分形维》

2 分形与分形维

《2.1分形概述》

2.1分形概述

分形 (fractal) 是数学家B.B.Mandelbrot于1975年最先创用的。迄今为止, 分形尚无严格的数学定义, B. B. Mandelbrot曾建议把分形定义为:局部以某种方式与整体相似的、即具备自相似性的几何形体。分形几何理论首先出现在地貌、地表形态、地质结构、材料结构等研究领域, 随后在宇宙星系的研究等方面也得到了广泛应用。

这里所说的自相似性既包括严格自相似性, 也包括统计自相似性。属于前一类的有Koch曲线、Peano曲线、Cantor栅栏 (Cantor集) 等。现实世界的几何体大多不具备严格的自相似性, 却可以有统计意义上的自相似性, 例如, 曲曲弯弯的海岸线便具有以Koch曲线为基本单元的统计自相似性。人和动物体内的肺、血管、胆管等的复杂分支结构也具有这种统计自相似性。

人体内的分形结构产生于胚胎发育过程, 与生命同在。例如, 健康人的心率即便在静止状态下也有很大波动, 心率的时序曲线呈现出忽高忽低的波浪状, 看起来似乎是完全随机的, 但如果观察较长时段的这种变化图像就会发现, 从一批心率数据中还是能找出其规律性。例如, 观察一段几个小时的心率时序图即可见到, 各个较短时段的心率时序图像都同整个时段的总图像具有统计自相似性。

维数是几何对象的重要基本特征量, 分形理论中的这一特征量称为分形维 (分形维数, 分数维) 。

《2.2分形维》

2.2分形维

在传统的Euclid几何学中, 维数是一个很重要的概念。其涉及的Euclid维都是整数:点为零维, 线为1维, 面为2维, 体 (空间) 为3维;在Einstein相对论中, 则讨论了时-空4维世界;热力学所研究的n维相空间, 其维数n可以是大于4的整数;至于科学技术研究中涉及的n维因素空间, 影响事件进行过程的因素每增加一个, 即相当于维数增1, 这儿的n也是整数。

以上所述, 涉及的都是规则的几何图形。然而对于广泛的复杂几何图形, 例如, 对Koch曲线、Sierpinski地毯等, 整数维的概念已不能适合所有场合, 必须将维数的概念加以拓展, 使之可以具有分数值, 即引入分形维 (或称分维数) , 例如, Koch曲线的分形维是1.261 8, Sierpinski 地毯的分形维是1.892 8, Vicsek图形的分形维是1.465 0, 等等^[1]。

分形维的种类很多, 例如, 相似维, 容量维 (Kolmogorov维) , 关联维等等^[2,3,4], 常用的前二者的定义如下:

1) 相似维

一般地, 当把某个图形的长度 (或标度) 放大1/r倍时 (r>1) , 得到N (r) 个和原图形相似的小图形。如果 N=r^-D, 则指数D就具有维数的意义, 称为相似维^[5], 记作

$D_{s} = \frac{\ln Ν (r)}{\ln (1 / r)} 。 (1)$ $D_{s} = \frac{\ln Ν (r)}{\ln (1 / r)} 。 (1)$

例如, 将直线段、正方形、立方体的边二等分 (即r=1/2) , 如图1所示, 则这时该直线段、正方形、立方体可看成为分别由2¹, 2², 2³个严格相似的小图形组成。这里的相似维1, 2, 3分别与有关图形的Euclid维 (即经验维) 是一致的。

《图1》

图1 线、面、体的划分与其维数

Fig.1 Partition of a line, a square, a cubic and their dimension values

2) 容量维

设所要考虑的图形是d维Euclid空间R^d中的有界子集, 并用半径为ε的d维小球来覆盖该图形, 以N (ε) 表示这类小球的个数的最小值, 则容量维D_c为^[6]

$D_{c} = \lim_{ε \to 0} \frac{\ln Ν (ε)}{\ln (1 / ε)} 。 (2)$ $D_{c} = \lim_{ε \to 0} \frac{\ln Ν (ε)}{\ln (1 / ε)} 。 (2)$

实际操作中, 常将上述覆盖用的d 维小球改为d维小盒, 由式 (2) 得到的便是盒子维 (也称盒维数) 。

当研究2维平面上的图形时, 用于覆盖的可以是2维小格子, 这时由式 (2) 得到的是格子维D_g。

笔者采用简单有效的格子维对未房颤时和发生房颤时的ECG曲线进行划分, 随之求得分形维。

《3 房颤的分形维分析》

3 房颤的分形维分析

ECG的时域波形一般由P波、QRS-T复合波等组成。房颤发生时, 其中的P波演变为 f 波 (房颤波) 。利用频数直方图技术将ECG图像中的QRS-T波剔除掉, 留下的便是代表心房活动的P波 (未房颤时) 或f波 (房颤时) ^[7], 如图2所示。

《图2》

图2 P波和f波的图像

Fig.2 The curves of P wave and f wave

《3.1数据的归一化处理》

3.1数据的归一化处理

图2所示的P波曲线、f波曲线的自变量单位是秒 (s) , 因变量单位是毫伏 (mV) , 两者的量纲不同, 于是, 分别提取每一心动周期的P波数据、f波数据, 并对之进行归一化处理, 使其转化为无量纲的物理量, 而后求出分形维。

采用的归一化方法如下:

1) 时间值的归一化以T表示一个心动周期的时间值 (s) , 令均匀采样点的序号数i 依次取值为0, 1, 2, …, N -1, 共N个采样点, (N -1) 个采样间隔, 各采样点的时间值t (i) 的计算式为

$t (i) = \frac{i Τ}{Ν - 1} ‚ i = 0, 1, 2, \dots ‚ Ν - 1 。 (3)$ $t (i) = \frac{i Τ}{Ν - 1} ‚ i = 0, 1, 2, \dots ‚ Ν - 1 。 (3)$

而t (i) 的归一化取值n (i) 为

$n (i) = \frac{t (i)}{Τ} = \frac{i}{Ν - 1}, i = 0, 1, 2, \dots ‚ Ν - 1 。 (4)$ $n (i) = \frac{t (i)}{Τ} = \frac{i}{Ν - 1}, i = 0, 1, 2, \dots ‚ Ν - 1 。 (4)$

2) 幅值的归一化令该心动周期内的幅值y (i) 的最大值和最小值分别为 y_max和 y_min, 则

$\begin{array}{l} y_{\max} = \max {y (i)}, i = 0, 1, 2, \dots, Ν - 1, \\ y_{\min} = \min {y (i)}, i = 0, 1, 2, \dots, Ν - 1 。 (5) \end{array}$ $\begin{array}{l} y_{\max} = \max {y (i)}, i = 0, 1, 2, \dots, Ν - 1, \\ y_{\min} = \min {y (i)}, i = 0, 1, 2, \dots, Ν - 1 。 (5) \end{array}$

于是, 归一化幅值 Y (i) 为

$Y (i) = (y (i) - y_{\min}) / (y_{\max} - y_{\min}) 。 (6)$ $Y (i) = (y (i) - y_{\min}) / (y_{\max} - y_{\min}) 。 (6)$

经此处理后, Y (i) , n (i) 都化为[0, 1] 区间内的无量纲值。

《3.2求分形维的方法》

3.2求分形维的方法

利用3.1节的方法, 分别对每一心动周期的P波、f波的相关数据进行归一化处理, 得到归一化后的P波数据或f波数据, 据此绘出归一化P波曲线、归一化f波曲线。

先讨论如何计算P波曲线的分形维。

P波曲线是处于2维Euclid平面上, 于是, 用边长为a的小方格将归一化P波曲线所在的平面进行划分, 以N (a) 表示包含该曲线的小方格数目。通过改变方格子的边长a的值, 可得到一系列的数据对{a, N (a) }, 再据此算出相应的一系列的数据对{ln a, ln N (a) }, 并对之进行最小二乘法拟合, 便得出

$\ln Ν (a) = l n k - D \ln a 。 (7)$ $\ln Ν (a) = l n k - D \ln a 。 (7)$

式 (7) 的指数形式为

$Ν (a) = k a^{- D} 。 (8)$ $Ν (a) = k a^{- D} 。 (8)$

显然, D相当于维数, 即所要求的分形维 (确切地说, 是格子维) 。

f波的分形维的求法也是如此。

《3.3计算结果》

3.3计算结果

房颤的分形维分析中所使用的ECG数据均来自MIT/BIH心电数据库, 未发生房颤时的和房颤时的数据各半, 共32例ECG曲线, 每例各取其50个心动周期进行分析, 即利用上述方法分别求出每例ECG曲线的50个D值, 再算出这50个D的均值D_mean 。

表1所列的16例是未房颤时的和16例发生房颤时的分维数均值D_mean。由表1可见, 房颤时f波的D_mean值小于未房颤时P波的D_mean值。

表1P波、f波的分形维

Table 1 The fractal-dimension values of P waves and f waves

《表1》

P波				f波
编号	D_mean	编号	D_mean	编号	D_mean	编号	D_mean
1	1.840 2	9	1.843 6	1	1.673 4	9	1.653 4
2	1.852 4	10	1.872 4	2	1.668 9	10	1.672 6
3	1.852 8	11	1.833 9	3	1.687 5	11	1.661 5
4	1.848 9	12	1.881 7	4	1.691 8	12	1.690 3
5	1.868 4	13	1.846 2	5	1.631 3	13	1.651 7
6	1.885 8	14	1.864 3	6	1.690 9	14	1.673 1
7	1.869 1	15	1.855 0	7	1.676 6	15	1.632 9
8	1.839 4	16	1.843 1	8	1.696 6	16	1.671 1

在利用此16 (例) ×50 (周期) P波数据和16 (例) ×50 (周期) f波数据、即总共 (16+16) ×50 = 1 600个数据进行房颤诊断时, 若根据P波某个周期的数据算出的分形维D的值显示为房颤属性, 则这样的结果为假阳性;而若根据f波某个周期的数据算出的分形维D的值显示为非房颤属性, 则这样的结果为假阴性。相应于这1 600个数据, 计算结果的分布情况见表2。

表2算出的D值的分布

Table 2 The distribution of D values

《表2》

数据总数	阳性数		阴性数		正确率^* /%
数据总数	真阳性	假阳性	真阴性	假阴性	正确率^* /%
1 600	646	121	679	154	82.8

* 正确率 = (646 + 679) /1 600

《4 信息熵[8,9,10,11,12]》

4 信息熵[8,9,10,11,12]

《4.1信息熵的概念》

4.1信息熵的概念

熵是描述复杂大系统的一个很有用的物理量, 熵的大小是无序度的一种度量, 适于用来研究随机过程。从Clausius的热力学熵, 到Boltzmann的统计熵, 再到Shannom的信息熵, 熵分析法已广泛应用于物理、化学、生物学、医学、经济学等众多自然科学和社会科学领域。

信息熵 (又名Shannom熵) 是C.E. Shannom于1948年在其论文“通信的数学理论 (A mathematical theory of communication) ”中提出来的, 其定义为

$Η (p_{j}) = - \sum_{j = 1}^{Ν} p_{j} \lg p_{j} ‚ j = 1, 2, \dots ‚ Ν ‚ (9)$ $Η (p_{j}) = - \sum_{j = 1}^{Ν} p_{j} \lg p_{j} ‚ j = 1, 2, \dots ‚ Ν ‚ (9)$

式中, p_j 是第j个消息的出现概率, H (p_j) 是信息熵 (即信源发出的平均信息量) 。

信息熵是分析随机过程的通信系统所含的信息量、衡量其不确定性 (即随机性) 的重要方法。

《4.2信息通信系统与医学》

4.2信息通信系统与医学

医学诊断治疗过程和护理过程实质上可以看成是对病人所患疾病进行相关信息的采集、传输、存储、处理的过程。

在临床医疗实践中, 当医务人员首次接触病人时, 即开始对病人的信息进行采集和处理。医生通过各种诊断手段, 运用自身的专业知识和技能, 对病人所患疾病进行分析综合, 进行诊断, 拟定治疗方案并实行之。护士通过各种途径对病人身体状况进行判断, 做出护理诊断, 依据医嘱进行各项护理活动, 执行护理程序。整个医疗活动过程既有病人 (信源) 的生理和病理信息的发送过程, 又有医生护理人员 (信宿) 对病人信息的接收过程, 而具体的医疗活动行为可视为广义上的信息流通的信道。在抽象层次上, 可以把临床诊断过程中的信息流视为一种通信系统, 其框图如图3所示^[11,12]。

在这个系统中, 信源信息X的发生并不是一次性单向传递过程, 其总量是在医务人员主动行为的影响下形成的, 实际是一个不断增加的值。或者说, 它同时也是时间的函数, 现实中表现为医生病案记录中所包含的病人全部信息, 其中, 既包括阳性结果, 也包括阴性结果;既包括有效信息, 也包括冗余信息。信宿收到的信息Y中既有信源信息X (这是主体) , 也有噪声Z, 必须从Y中将Z剔除, 从而得到与X相应的有效信息, 据以计算信息熵H。

《4.3房颤的信息熵分析》

4.3房颤的信息熵分析

将信息熵理论用于房颤信息的处理和分析。

房颤的ECG检测过程是一种医学信息的流通传递系统, 如图3所示。图3中, 信源发出的信息是生理、病理信号, 而信宿收到的是从ECG信号获得的并剔除了传输噪声的P波、f波信息。

《图3》

图3 医学诊断过程的通信系统框图

Fig.3 The medical diagnosis as a communication system and its block diagram

在根据式 (9) 计算信息熵时, 必须知道各个 p_j 的值 (j=1, 2, … , N) 。求 p_j 的方法是:由很窄的狭长条块将每一心动周期的、经过归一化处理的P波曲线和f波曲线予以分割, 再根据各条块中含有的图形像素数d (j) _{j = 1, 2, …, N} 来计算概率

$p_{j} = \frac{d (j)}{\sum_{j} d (j)}, j = 1, 2, \dots ‚ Ν 。 (10)$ $p_{j} = \frac{d (j)}{\sum_{j} d (j)}, j = 1, 2, \dots ‚ Ν 。 (10)$

对于P波、f波, 算得的信息熵如表3所示。

房颤的信息熵分析中所使用的ECG数据均来自MIT/BIH心电数据库, 未发生房颤的和房颤的各16例。对每例ECG曲线, 均取其50个心动周期, 分别计算出每个周期的信息熵, 得到50个H值, 而后求其平均值H_mean。表3所列的是16例未房颤时的和16例发生房颤时的H_mean。由表3可见, 房颤时f波的H_mean值显著大于未房颤时P波的H_mean。

相应于1 600个数据的计算结果见表4。

表3P波、f波的信息熵

Table 3 The Shannom entropy of P waves and f waves

《表3》

P波				f波
编号	H_mean	编号	H_mean	编号	H_mean	编号	H_mean
1	5.736 2	9	5.695 5	1	8.355 7	9	8.056 9
2	6.177 6	10	5.901 1	2	8.284 3	10	8.212 6
3	5.702 0	11	5.952 2	3	8.183 3	11	8.278 0
4	6.156 4	12	6.063 9	4	8.197 8	12	8.097 6
5	5.969 7	13	5.966 4	5	8.370 7	13	8.338 1
6	6.077 9	14	6.059 7	6	8.409 5	14	8.457 3
7	5.966 4	15	5.829 0	7	8.432 5	15	8.376 6
8	6.059 7	16	5.898 9	8	8.208 5	16	8.224 4

表4算出的H值的分布

Table 4 The distribution of H values

《表4》

数据总数	阳性数		阴性数		正确率^* /%
数据总数	真阳性	假阳性	真阴性	假阴性	正确率^* /%
1 600	685	103	697	115	86.4

* 正确率 = (685 + 697) /1 600

《6 结语》

6 结语

利用分形维分析法对P波、f波的信息进行了处理和分析, 利用2维盒子维 (即格子维) 计算了16例未房颤时每一心动周期的 (每例取50个周期) 、16例房颤时每一心动周期 (每例也取50个周期) 的分形维D, 据以计算相应的50个D的均值D_mean。

未房颤时与房颤时的值之间的明显差异表明, 用分形维分析很容易判断是否发生了房颤。

利用信息熵 (shannom熵) 对P波、f波的信号进行处理和分析, 用上述计算D_mean的同样方法, 可以得到相应的32个均值H_mean。

未房颤时与房颤时的H_mean值间的显著差异表明, 用信息熵分析法很容易判断是否发生了房颤。

展示更多