1.   引　言

时间可逆指系统特征不随时间逆序而发生变化的性质, 反之则称时间不可逆(time irreversibility, TIR)^[1]. 时间不可逆是衡量系统非平衡性特征的重要方法, 其量化需要计算正反序列或是序列对称向量的联合概率差异. 当前TIR分析主要采用符号化方式^[2], 由于排列类型^[3]有效表征了序列的结构特征, 且无需模型设定, 其在TIR分析中受到广泛的关注^[4–9].

Martinez等^[4]利用排列概率的分布差异量化时间不可逆, 有效地揭示了经济、气象、生理等复杂系统的非平衡性特征. Zanin等^[5,6]通过经济信号验证了正反序列的排列概率差异在复杂系统分析中的优势. Yao等^[1,7–9]系统研究了排列时间不可逆(permutation TIR, pTIR)在复杂系统分析中的影响因素, 其中特别关注了空排列的影响以及等值元素的物理意义. 这些研究为排列时间不可逆理论的理解和现实应用提供了重要的基础. 排列类型作为序列结构的符号化表达方式, 在简化空间特征的同时, 也影响时间不可逆的分析结果. 首先, 原始排列和幅度排列两种基本排列类型表征的空间结构存在本质差异^[9], 其中幅度排列是序列空间结构的直接表征形式. 在时间不可逆的量化中, 特别是基于对称向量概率差异性的简化计算, 原始排列可能会带来偏差^[8]. 其次, 空排列(forbidden permutation)的出现对于pTIR排列概率差异的计算是一个不利的因素^[1]. Amigo等^[10–12]广泛研究了空排列的特征, 并结合现实应用提出了多种非线性分析方法. 后续学者从多角度研究了空排列相关的系统信息, 如混沌动力学、相关性等^[13,14]. 然而在pTIR中, 空排列会导致部分基于除法的参数(如Kullback-Leibler距离, 即相对熵)不适用于排列概率差异的计算^[1,7–9,15]. 除此之外, 等值元素在排列分析中通常会被忽略, 其在pTIR分析中起到了关键的作用. 等值是构建全面排列结构的关键因素^[9], 对排列的概率分布和排列分析都会产生重要的影响, 甚至造成完全相反的结果^[1,16,17]. 等值会生成自对称向量, 在时间不可逆的分析中包含明确的物理意义, 即时间可逆和时间对称^[1,7–9]. 为优化等值元素的排列类型, Bian等^[18]提出了等值坐标最小化的排列构建方法, Zunino^[19]和Cuesta等^[20]通过排列熵论证了等值元素对复杂度量化的影响. 遗憾的是, 排列类型是向量空间结构的粗粒化表示方式, 无法体现元素之间的绝对位置关系, 因此当前对排列分析的优化并不能解决其在精确度方面不足的问题.

为此, 本文提出一种基于模糊排列的时间不可逆(fuzzy permutation time irreversibility, fpTIR)分析方法, 将元素之间的差值转化为排列类型的隶属度, 实现对向量结构特征的精准表征, 提高时间不可逆性分析的精确度. 首先介绍时间不可逆的概念和幅度排列的构建方式, 然后详述向量排列隶属度和fpTIR的分析方法. 接下来, 构建模型序列及其代替数据, 验证fpTIR的有效性, 然后和模糊排列熵(fuzzy permutation entropy, fPEn)作对比, 研究PhysioNet公开数据库中的心衰(congestive heart failure, CHF)患者、健康年轻人和健康老年人心率的复杂特征. 在此基础上对fpTIR的分析结果进行讨论和总结.

2.   理　论

2.1   时间不可逆概念

在统计概念上, 如果系统在时间逆尺度上的统计特征不变, 该系统是时间可逆的, 反之则是时间不可逆. Weiss^[21]指出, 给定平稳信号$X\left(t\right)$, 如果$\left\{X\right({t}_{1}), X({t}_{2}), \cdots , X({t}_{m}\left)\right\}$与$\{X({-t}_{1}), X(-{t}_{2}), \cdots , X(-{t}_{m})\}$对于所有的${t}_{1}, {t}_{2}, \cdots , {t}_{m}$和$m$都有相同的联合概率分布, 则$X\left(t\right)$是时间可逆的, 否则$X\left(t\right)$时间不可逆. Kelly^[22]的研究表明, 对于所有的$n$和$m$, 如果$\left\{{X(t}_{1}\right), {X(t}_{2}), \cdots {, X(t}_{m})\}$和$\{{X(-t}_{1}+ n), {X(-t}_{2}+n), \cdots {, X(-t}_{m}+n)\}$有着相同的联合概率分布, 则序列是时间可逆的. 特别是当$n={t}_{1}+ {t}_{m}$时, $\left\{{X(t}_{1}\right), {X(t}_{2}), \cdots {, X(t}_{m})\}$和其对称形式$\left\{{X(t}_{m}\right), \cdots , {X(t}_{2}\left){, X(t}_{1}\right)\}$有相同的概率分布.

2.2   模糊幅度排列

原始排列(original permutation)和幅度排列(amplitude permutation)是两种基本的排列类型^[9], 其中原始排列为排序序列元素在原始序列中的位置, 幅度排列为原始序列元素在排序序列中的位置. 原始排列可以定位序列排序后的元素坐标, 广泛应用于科研和数值分析软件, 如MATLAB的“sort”、PYTHON的“argsort”和R语言的“order”函数等. 幅度排列是序列空间结构的直接表征形式, 能有效避免排列时间不可逆在信号分析中可能存在的偏差^[1,9], 是本文采用的排列类型.

设长度为L的向量$X\left(i\right)=\left\{x\left({i}_{1}\right), x\left({i}_{2}\right), \cdots , x\left({i}_{n}\right), \cdots , x\left({i}_{L}\right)\right\}$, 将其元素按照从小到大的顺序排序得$X\left(j\right)=\{x\left({j}_{1}\right), x\left({j}_{2}\right), \cdots , x\left({j}_{n}\right), \cdots , x({j}_{L})\}$, 其中$x({j}_{1}) < x ({j}_{2}) \cdots < x ({j}_{n}) < \cdots < x({j}_{L})$. 幅度排列由原始向量$X\left(i\right)$的元素$x\left({i}_{n}\right)$在排列向量$X\left(j\right)$中的位置${j}_{in}$构成, 即${\pi }_{j}= ({j}_{1 n}, {j}_{2 n}, \cdots , {j}_{in}, \cdots , {j}_{Ln})$. 等值元素是构建排列类型的关键, 其存在可能产生时间可逆或时间对称的自对称向量, 对排列分析有重要的影响^[1,7–9]. 为有效表征向量的等值元素, 可采用等值坐标统一化的方式对幅度排列进行修正^[9]. 首先将等值元素先按照出现次序排列, 如$x\left({j}_{1}\right)= x\left({j}_{2}\right) < \cdots < x\left({j}_{3}\right)=x\left({j}_{4}\right)=x\left({j}_{5}\right)$. 然后, 对等值元素在排序中的坐标统一化, 即将等值元素的坐标修改为同组的相同值, 如最大值$x\left({j}_{2}\right)= x\left({j}_{2}\right) < \cdots < x\left({j}_{5}\right)=x\left({j}_{5}\right)=x\left({j}_{5}\right)$或Bian等^[18]的最小值$x\left({j}_{1}\right)=x\left({j}_{1}\right) < \cdots < x\left({j}_{3}\right)=x\left({j}_{3}\right)=x\left({j}_{3}\right)$, 并修改等值元素在幅度排列中的坐标.

排列类型是向量结构的粗粒化表征方法, 无法体现元素之间的绝对差异. 为解决该问题, 首先定义标准向量: 计算排序向量的差分向量, 如果其标准差为零, 则为标准向量, 否则为特异向量. 随着差分向量标准差的增大, 向量的特异性增强, 也就是与标准向量之间的差异增大. 为了量化向量的特异性, 在幅度排列中引入隶属度, 构建模糊排列. 对于升序排列后的向量$X\left(j\right)$, 计算相邻元素的差值${x}' ({j}_{i}) = x({j}_{i+1}) - x({j}_{i})$, 构建升序差值向量${X}' (j) = \left\{{x}' ({j}_{1}),~ {x}' ({j}_{2}),~ \cdots , {x}' ({j}_{i}),~ \cdots ,~ {x}' ({j}_{m-1})\right\}$, 计算${X}'(j)$的标准差$\sigma \left(j\right)$. 如果$X\left(j\right)$包含等值元素则将其认定为一个元素, 从而避免对排列隶属度计 算的影响. 根据模糊理论^[23], 本文利用负指数函 数^[24,25]计算排列向量$X\left(j\right)$隶属排列类型的概率${p}_{j}$, 如(1)式所示, 其中$\alpha$为可调的控制参数:

如果向量$X\left(j\right)$包含两个元素, 可将$\sigma \left(j\right)$修正为向量标准差和其所在序列标准差的比值; 如果$X\left(j\right)$全为等值元素, 等价为一个元素, 隶属度${p}_{i}$设定为1. 模糊排列为排列类型及其隶属度的组合形式${\pi }_{f}= ({\pi }_{j};{p}_{j})$. 由此可知, 标准向量升序序列$\sigma \left(j\right)= 0$, 即${p}_{j}= 1$; 特异向量$\sigma \left(j\right) > 0$, 则${p}_{j} < 1$, 并且${p}_{j}$随着特异性的增强而降低. 以5元素的向量$X\left(1\right)=\left\{{\mathrm{9, 3}}, {\mathrm{7, 1}}, 5\right\}$和$X\left(2\right)=\{8.5, 3.1, 6.5, 1.3, 5.2\}$为例, 其幅度排列都是${\pi }_{j}= (5, 2, 4, 1, 3)$, 升序差分向量分别为${X}'\left(1\right)=\left\{2, 2, 2, 2\right\}$和${X}'\left(2\right)=\{1.8, 2.1, 1.3, 2\}$, 标准差为$\sigma \left(1\right)= 0$和$\sigma \left(2\right)= 0.35$. 将(1)式的控制参数设为$\alpha = 1$, 则$X\left(1\right)$和$X\left(2\right)$的隶属度为${p}_{1}= 1$和${p}_{2}= 0.7$, 两者的模糊排列分别为${\pi }_{1}= ({\mathrm{5, 2}}, {\mathrm{4, 1}}, 3;1)$和${\pi }_{2}= ({\mathrm{5, 2}}, {\mathrm{4, 1}}, 3;0.7)$.

2.3   基于模糊排列的时间不可逆

对于序列$X\left(t\right) = \{x\left(1\right), x\left(2\right), \cdots , x\left(t\right), \cdots , x(L)\}$, 首先进行维度$m$和延迟$\tau$的向量重构, 得向量序列${X}_{\tau }^{m}\left(t\right) = \{{x}_{\tau }^{m}\left(1\right), {x}_{\tau }^{m}\left(2\right), \cdots , {x}_{\tau }^{m}\left(t\right), \cdots , {x}_{\tau }^{m}(L- (m-1)\tau )\}$, 其中向量为${x}_{\tau }^{m}\left(t\right) = \{x\left(t\right), x\left(t+\tau \right), \cdots , x\left(t+(m-1)\tau \right)\}$. 然后, 通过向量${x}_{\tau }^{m}\left(t\right)$的幅度排列和隶属度生成模糊排列类型, 构建模糊排列序列${\pi }_{\tau }^{m}\left(t\right) = \{{\pi }_{{\mathrm{f}}}\left(1\right), {\pi }_{{\mathrm{f}}}\left(2\right), \cdots , {\pi }_{{\mathrm{f}}}\left(t\right), \cdots , {\pi }_{{\mathrm{f}}}(L- (m-1)\tau \left)\right\}$. 时间不可逆为序列$X\left(t\right)$和其时间逆序$X(-t)$的模糊排列概率差异, 其中模糊排列的概率$p\left(\pi \right)$为幅度排列隶属度${p}_{j}$的累加与排列数量的比值, 如(2)式所示. 考虑到空排列的存在, 本文利用基于减法的Y_s计算模糊排列概率差异^[1,7–9,15], 如(3)式所示, 其中${p}_{{\mathrm{f}}} (\pi)$和${p}_{{\mathrm{b}}} (\pi)$分别为正反序列模糊排列的概率分布. 模糊排列时间不可逆的分析流程如图1所示.

图 1  模糊排列时间不可逆的流程图

Fig. 1.  Flowchart of fuzzy permutation time irreversibility.

下载: 全尺寸图片 幻灯片

作为对比, 利用香农熵计算模糊排列概率分布的平均信息量^[26], 即模糊排列熵(fuzzy permutation entropy, fPEn), 如(4)式所示. 通过${\mathrm{l}}{\mathrm{n}}\left[k\right(m\left)\right]$对fPEn进行归一化, 其中$k\left(m\right)$为维度$m$时最大的排列类型数量^[18], 如当$m$为2, 3, 4和5时, $k\left(m\right)$分别为3, 13, 73和501.

3.   结　果

本节首先构建模型序列, 根据代替数据理论对fpTIR的有效性进行验证, 然后对不同生理状态的心率信号进行分析.

3.1   模型序列验证

通过logistic (${x}_{t+1}=r\cdot {x}_{t}(1-{x}_{t})$, $r$ = 4)和Henon(${x}_{t+1}={y}_{t}+1-\alpha {x}_{t}^{2}$, ${y}_{t+1}=\beta {x}_{t}$, $\alpha$ = 1.4和$\beta$ = 0.3)模型构建混沌序列, 利用一阶自回归模型(AR1, ${x}_{t+1}=\delta {x}_{t}+{\xi }_{t}$, 其中${\xi }_{t}$为高斯信号, $\delta$ = 0.3)构建线性序列, 然后采用Schreiber的改进幅频扰乱(improved amplitude adjusted Fourier transform, iAAFT)方法^[27,28]为每组模型序列构建500组代替数据. 根据代替数据理论, 如果模型序列的fpTIR在代替数据的2.5%和97.5%分位数之间, 则模型序列和代替数据没有显著差异; 反之, 若模型序列的fpTIR大于代替数据的97.5%分位数或是小于2.5%分位数, 则模型序列和代替数据有显著差异. Logistic, Henon和AR1模型序列及其代替数据的fpTIR和fPEn如图2所示.

图 2  Logistic, Henon和AR1模型序列及其代替数据的fpTIR和fPEn. 隶属度控制参数$\alpha$= 1, 向量重构参数为$m$ = 2—5, $\tau$ = 1. 图中fpTIR-97.5, fpTIR-2.5, fPEn-97.5和fPEn-2.5分别表示代替数据fpTIR和fPEn的97.5%和2.5%分位数

Fig. 2.  The fpTIR and fPEn of the logistic, Henon, and AR1 series and their surrogates. Control parameter of membership degree is $\alpha$ = 1, dimension and delay are $m$ = 2–5, $\tau$ = 1. The fpTIR-97.5, fpTIR-2.5, fPEn-97.5 and fPEn-2.5 denote the 97.5% and 2.5% of fpTIR and fPEn of surrogate data.

下载: 全尺寸图片 幻灯片

由图2可知, fpTIR和fPEn都可有效地表征模型序列的复杂特征, 然而两者在混沌序列中存在差异. 如图2(a), (b)所示, logistic和Henon混沌序列的fpTIR高于其代替数据的97.5%分位数, 而其fPEn却低于代替数据的2.5%分位数. fpTIR和fPEn对模糊排列概率分布的分析方式不同, 这是两者在混沌序列代替数据分析中存在差异的根本原因. 通过iAAFT生成的代替数据具有较高的随机性、不可预测性, 其模糊排列的概率分布更加均匀、差异性更小, 即更低的非平衡性和更高的信息量. 因此, 混沌序列的fpTIR显著高于代替数据, 而其fPEn却显著低于代替数据. 如图2(c), (d)所示, 线性AR1模型的fpTIR和fPEn都在其代替数据的2.5%和97.5%分位数之间. 根据代替数据理论, logistic和Henon序列是非线性的, 而AR1模型序列是线性的, 这和三组模型序列的特征一致.

3.2   心率信号fpTIR分析

人体心率受到多种内部(神经、体液、激素等)和外部(声音、光照、温湿度等)因素的影响, 表现出高度的动态变化和复杂特征. 接下来, 利用PhysioNet数据库^[29]中的心衰患者、健康老年人和健康年轻人心率对fpTIR进行分析. 44组心衰心率((55.5±11.4)岁, 年龄22—79岁)信号采集自chfdb和chfdb2两个数据集^[30], 20名健康老年人((74.5±4.4)岁, 年龄68—85岁)和20名健康年轻((25.8±4.3)岁, 年龄21—34岁)心率采集源自Fantasia数据库^[31], 其中健康老年和健康年轻心率包含相同数量的男性和女性受试者. 三组心率信号的复杂动力学特征符合复杂度丢失理论^[32]: 健康年轻机体的心率调节系统有高度的复杂特征, 随着老年人心率调节功能的降低, 其复杂度也会随之降低, 而心衰会导致心率呈现更低的复杂特征.

首先研究fpTIR在心衰、健康老年和健康年轻心率信号分析中的特点. 幅度排列构建的维度$m$设为2和3, 延迟$\tau$为1和2, 模糊排列隶属度$\alpha$ = 0.1—1 (步长0.1), 3种心率的fpTIR和pTIR如图3所示. 从图3可知, 心率的fpTIR和pTIR符合复杂度丢失理论, 其中fpTIR能够更加准确地表征心率信号的非平衡性特征. 当$m = 2$, $\tau= 1$ 时, fpTIR对3种心率的识别效果优于pTIR; 当$m = 3$, $\tau = 1$, 心衰心率的pTIR异常地高于健康老年人, 而fpTIR正确反映了心率复杂特征的变化, 即健康年轻心率>健康老年心率>心衰心率. 3种心率的fpTIR并没有受到隶属度参数$\alpha$的显著影响, 特别是当$m= 3$时, 心率的fpTIR在[0.1, 1.0]控制参数区间内并没有发生明显变化.

图 3  心衰、健康老年和健康年轻心率的fpTIR和pTIR (均值±标准误差), 幅度排列构建维度$m$ = 2和3、延迟$\tau$ = 1和2, 模糊排列隶属度$\alpha$ = 0.1—1 (步长0.1)

Fig. 3.  The fpTIR and pTIR (mean ± standard error) of the CHF, healthy elderly, and healthy young heartbeats, dimension and delay of amplitude permutations are $m$ = 2 and 3, $\tau$ = 1 and 2, $\alpha$= 0.1–1 with step size of 0.1.

下载: 全尺寸图片 幻灯片

根据心率的fpTIR分析结果, 选择排列隶属度控制参数$\alpha$ = 0.1, 并结合模糊排列的概率分布特征, 进一步对比性分析fpTIR和fPEn在3种心率信号的表现, 如图4所示. 由图4(a), (b)可知, 心衰、健康老年和健康年轻心率的fpTIR和fPEn呈现相反的结果: fpTIR和复杂度丢失理论一致, 而fPEn和复杂度丢失理论相反(心衰患者>健康老年人>健康年轻人). 和代替数据相同, fpTIR和fPEn在心率信号中的矛盾性结果是因为两者对排列概率的分析方式不同. 当$m$ = 2时, 有上升(12)、下降(21)和等值(11) 3种幅度排列, 其中等值为时间可逆. fpTIR为上升和下降的概率差异, fPEn则为3种排列包含的平均信息量. 如图4(c)所示, 心衰患者心率的升降概率差异最小, 即最小的fpTIR, 而3种排列的概率差异也是最小, 即最大的平均信息量和fPEn. 与之相对应, 健康年轻心率升降概率差异和3种排列的概率差异都最大, 因此具有最大的fpTIR和最小的fPEn. 当$m =3$时, 共有13种幅度排列, 其中自对称的131, 212和111意味着时间可逆. 如图4(d)所示, 3种心率的模糊排列概率分布相对复杂, 但和$m= 2$ 的情况一致: 心衰心率对称排列的概率差异和所有排列类型的概率差异最小, 因此具有最小的fpTIR和最大的fPEn, 而健康年轻心率有最大的fpTIR和最小的fPEn. fpTIR正确地表征了与心率相关的复杂度丢失理论, 符合心率调节的生理、病理特征, 不过fPEn对3种心率有更好的区分效果, 并且对维度和延迟的选择也不敏感, 因此fpTIR和fPEn以及它们之间的关系都值得进一步研究.

图 4  心衰、健康老年和健康年轻心率的fpTIR、fPEn和排列概率分布(均值±标准误差), 模糊排列隶属度$\alpha$ = 0.1　(a), (b) $m2\tau 1$表示幅度排列维度2、延迟1, 并以此类推; (c), (d) 幅度排列的维度和延迟为$m$ = 2, $\tau$ = 2和$m$ = 3, $\tau$ = 2

Fig. 4.  The fpTIR, pTIR and probabilities of permutations (mean ± standard error) of the CHF, healthy elderly, and healthy young heartbeats, in membership degree of amplitude permutations, $\alpha$ = 0.1: (a), (b) $m2\tau1$ represents dimension of 2 and delay of 1, and so on; (c), (d) dimension and delay of amplitude permutations are $m$ = 2, $\tau$ = 2 and $m$ = 3, $\tau$ = 2.

下载: 全尺寸图片 幻灯片

通过模型数据验证和心率信号分析, fpTIR能够有效地表征复杂系统的特征, 其利用向量元素之间的绝对差异对排列类型进行优化, 提高了复杂系统非平衡性分析的精确度.

4.   讨　论

首先对fpTIR和fPEn在复杂系统分析中的差异进行讨论. fpTIR是正反序列排列或对称排列的概率差异^[1,4–8], fPEn是所有排列概率分布包含的平均信息量^[3,26,33], 这是导致两者在复杂系统分析中存在差异性结果的根本原因. 因此, 排列概率分布差异越小, 信号的熵值就越高, 而TIR则越低. PhysioNet心率信号^[30,31]的模糊排列概率分布特征明显, 更加突显了fpTIR和fPEn的差异性、甚至相反的结果. 从数值分析角度看, fpTIR和fPEn可能存在更加复杂的关系. 在极端情况下, 如果所有对称排列的概率分布相同, fpTIR为0, fPEn则随着排列概率分布的变化而变化; 在另一个极端情况下, 如果所有排列均为单排列(单排列的对称排列或逆序中的对应排列即空排列)^[1,7–9,15], 则fpTIR为1, fPEn仍随着排列概率分布的差异而变化. 本文时间不可逆和熵值的差异性和前期的研究结果一致, 在其他复杂生理信号(如睡眠脑电^[1]和心率^[16])分析中也有体现. 在复杂系统分析中, 时间不可逆和熵值的联合分析能够更加全面地探索系统特征, 有助于从不同的角度理解系统特征.

模糊排列分析具有更高的精确度, 能够更加准确地提取复杂系统的特征, 然而代价则是较高的实现复杂度和较低的抗干扰能力. 模糊排列的构建需要计算向量的绝对差值, 并在排列类型中引入隶属度, 这一过程增加了排列类型构建的复杂度. 模糊排列分析需要对排列的隶属度进行针对性的考虑, 这在现实信号分析中对算法编程和硬件设计提出了更高的要求. 此外, 排列隶属度是向量元素绝对差值的转化值, 信号微弱的变化会改变向量排列的隶属度, 因此模糊排列方法对干扰和噪声有更高的敏感度. 相比之下, 传统的排列分析是典型符号时间序列分析方法, 易于实现和抗噪声性能强是主要的优点. 因此, 在现实系统的分析中, 若信号采集干扰较小, 模糊排列方法能够提升复杂系统分析的精确度; 反之, 若信号噪声较大, 则传统的排列分析方法或许更为合适.

5.   结　论

本文通过模糊排列提取序列精确的符号动力学特征, 提高了复杂系统时间不可逆的分析效果, 总结如下.

1) 模糊排列有效提取了向量精确化的空间结构特征, 利用排列隶属度实现了对特异向量的有效表征. 根据代替数据理论, fpTIR可以有效地分析系统的复杂特征. 在复杂生理信号的分析中, fpTIR更加准确地表征了心率的复杂度丢失理论, 并且具有较高的稳定性. 需要注意的是, 模糊排列分析对算法编程和硬件设计实现都有较高的要求, 且对噪声更加敏感.

2) 在模型序列和心率的分析中, fpTIR和fPEn存在一定的差异性. 熵值和时间不可逆都是统计物理概念, 可通过信息论方式实现复杂系统的量化分析, 但是由于两者对概率分布的处理方式存在区别, 在信号分析中可能产生不同乃至相反的结果. fpTIR和fPEn的关联性分析连接了非平衡性和熵值复杂度, 为从更广泛的角度研究复杂系统特征提供了有价值参考.