胎儿心电（fetal electrocardiogram，FECG）是研究胎儿心脏电生理活动的一项客观指标，反映了胎儿在孕期中的成长和健康状况。关于FECG的提取方法主要有单通道提取和多通道提取两类。由于腹壁混合心电信号中除了各种噪声以外，还含有母体心电（maternal electrocardiogram，MECG）和胎儿心电两个心电信号。单通道提取属于欠定混合信号的分离问题，必须寻求辅助的先验信息支持。目前，单通道提取方法还存在许多关键性问题需要解决，研究成果也比较少。多通道提取方法可以通过导联数的增加，实现正定或超定混合信号的盲源分离。因此，盲源分离算法应用于腹壁混合心电信号的分离已成为研究热点。如奇异值分解（singular value decomposition，SVD）^[1]、主成分分析（principal component analysis，PCA）^[2-3]、独立成分分析（independent component analysis，ICA）^[4-6]和周期成分提取（periodic component extraction，PCE）^[7-8]等方法，均被用于FECG的提取。这些方法在瞬时线性混合模型下，基于源信号相互独立的假设，对混合心电信号进行分离。实际上，一些研究成果表明，腹壁混合心电信号中MECG和FECG之间并不总是相互独立的，而是具有弱相关性^[9]。这使得传统的盲源分离算法在胎儿心电提取上往往达不到预期的分离效果。此外，支持向量机（support vector machine，SVM）^[10]、人工神经网络（artificial neural network，ANN）^[11-14]等学习算法也被用于多通道FECG提取。但是，这些算法需要大量的训练过程才能获得较好的学习结果，心电信号的非线性和非平稳性，使得学习成果的普适性受到限制。

针对MECG与FECG的弱相关性影响盲源分离算法提取FECG效果的问题，本文基于经验模态分解（empirical mode decomposition，EMD）方法，消除MECG干扰，减弱相关性。采用准周期分量提取方法（quasi-periodic component extraction，QPCE）提取固有模态函数（intrinsic mode function，IMF）中的FECG信息，然后重构FECG，从而建立一个新的多通道FECG提取方法，命名为EMD–QPCE方法。通过EMD–QPCE方法对DaISy数据库、ADFECGDB数据库和Challenge 2013数据库中真实的腹壁混合心电信号进行FECG提取实验，验证方法的有效性。并与PCA、ICA和QPCE等传统盲源分离算法进行比较。结果表明，EMD–QPCE方法具有更好的成功提取率和提取质量。

1 EMD–QPCE胎儿心电提取方法

本文方法包含两部分主要内容：一是，基于EMD的去噪和MECG信息消除；二是，基于QPCE方法的IMF信息提取和FECG重构。EMD–QPCE方法提取FECG的基本过程是：将多通道腹壁混合心电信号分别进行EMD分解，得到各通道各自的一系列IMF，将各通道信号中相同阶数的IMF进行QPCE提取，得到FECG的各阶IMF，最后相加得到提取的最终FECG。

1.1 腹壁信号的EMD处理

EMD是1998年Huang等^[15]提出的一种基于信号自身数据驱动的自适应信号分解方法。EMD通过迭代筛选过程，将信号分解为由高频到低频的一系列固有模态函数IMF和最后余项，原始信号则由IMF和最后余项叠加重构。固有模态函数IMF必须满足两个条件^[15]：一是，IMF的极值点个数必须与过零点个数相等或最多相差一个；二是，IMF的极大值包络和极小值包络的均值等于零。

本文应用EMD进行腹壁混合心电信号处理的内容包括4个方面：

1）应用EMD阈值去噪方法，去除高频随机噪声；

2）基于FECG心率特征，剔除基线漂移和运动伪迹；

3）依据FECG和MECG的基频和幅值差异，消除IMF中MECG的QRS波群（MQRS）干扰；

4）利用FECG的QRS波群（FQRS）中显著的IMF，检测FECG的特征周期（准周期）。

1.1.1 高频噪声去除和运动伪迹基线漂移剔除

1）阈值法去高频随机噪声

基于EMD分解的高频噪声主要集中在低阶IMF中的特性，EMD阈值去噪方法首先要确定含高频噪声的IMF分界，对低于分界阶数的IMF进行阈值去噪^[16-18]。Donoho等^[19]提出了一个小波分解系数中随机噪声的估计方法；之后，Boudraa等^[17]将其应用于EMD的IMF中，式（1）为IMF中随机噪声标准差估计公式：

$\sigma _{{\rm{noise}}}^i = \frac{{{\rm{median}} \left( {\left| {s_{{\rm{imf}}}^i(n) - \bar s_{{\rm{imf}}}^i(n)} \right|} \right)}}{{0.674\;5}}$

(1)

式中， ${\rm{median}}( \cdot )$ 代表取中值， $s_{{\rm{imf}}}^i$ 为IMF的符号表示，i为IMF的阶数， $\sigma _{{\rm{noise}}}^i$ 代表第i阶 $s_{{\rm{imf}}}^i$ 中的噪声标准差估计， $\bar s_{{\rm{imf}}}^i$ 为第i阶 $s_{{\rm{imf}}}^i$ 的平均值。

另外， $s_{{\rm{imf}}}^i$ 的总体标准差可由式（2）计算：

$ {\;\;\;\;\;\;\;\;\sigma _{{\rm{imf}}}^i} = \sqrt {\frac{1}{{N - 1}}\sum\limits_{n = 1}^N {{{\left( {s_{{\rm{imf}}}^i(n) - \bar s_{{\rm{imf}}}^i(n)} \right)}^2}} } $

(2)

式中，N为信号的长度。

由于高频噪声存在于低阶IMF中，故使用式（3）确定噪声能量占比最小的IMF阶数。

${\;\;\;\;\;\;\;\;\;\;\;\;\;I_{{\rm{Htd}}}} = \min \left\{ {\mathop {\arg \min }\limits_{1 < i < M} \left\{ {\sigma _{{\rm{noise}}}^i/\sigma _{{\rm{imf}}}^i} \right\}} \right\}$

(3)

式中， ${I_{{\rm{Htd}}}}$ 为噪声标准差与总体标准差之比的第1个极小值对应的IMF阶数，M为IMF的最大阶数。则阶数小于等于 ${I_{{\rm{Htd}}}}$ 的IMF为含有高频随机噪声的IMF，本文将采用通用阈值中的硬阈值法对其进行去噪^[17,19]。

通用阈值计算：

$\varLambda _{{\rm{Td}}}^i = \sigma _{{\rm{noise}}}^i\sqrt {2\ln (N)} $

(4)

硬阈值去噪：

$ {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\tilde{s}}_{\rm{imf}}^{i}(n)=\left\{\!\begin{array}{l}{s}_{\rm{imf}}^{i}(n),\;\left|{s}_{\rm{imf}}^{i}(n)\right|\ge {\varLambda }_{\rm{Td}}^{i};\\ 0,\;\left|{s}_{\rm{imf}}^{i}(n)\right|<{\varLambda }_{\rm{Td}}^{i}\end{array}\right.$

(5)

式中， $i \le {I_{{\rm{Htd}}}}$ ， $\tilde s_{{\rm{imf}}}^i$ 为第i阶去噪后的IMF。

2）基线漂移和运动伪迹的去除

根据统计结果，FECG的心率通常在120～180 bpm范围内，相应的心频为2～3 Hz^[20]。频率低于FECG最小心频的IMF，应不含有FECG信息或者FECG信息可以忽略不计，则重构时可以舍弃这些IMF。

假设腹壁信号的采样频率为 ${f_{{\rm{sample}}}}$ ，FECG的心频为 ${f_{{\rm{fecg}}}}$ ，对应的周期长度为 ${T_{{\rm{fecg}}}}$ （采样点数），则：

${T_{{\rm{fecg}}}} = {f_{{\rm{sample}}}}/{f_{{\rm{fecg}}}}$

(6)

式中，选取 ${f_{{\rm{fecg}}}} = 120$ /60 Hz，为FECG心频统计范围的最小值，以保证不丢失FECG的有用信息。

记 $T_{{\rm{imf}}}^i$ 为第i阶IMF的平均周期，则对应于平均周期小于等于 ${T_{{\rm{fecg}}}}$ 的最大IMF阶数为：

${I_{{\rm{Ltd}}}} = \mathop {\arg \max }\limits_i \left\{ {T_{{\rm{imf}}}^i \le {T_{{\rm{fecg}}}}} \right\}$

(7)

阶数大于 ${I_{{\rm{Ltd}}}}$ 的IMF将被舍弃。

3）信号重构

高频噪声和基线漂移滤除后的信号重构如下：

$ {s_{{\rm{rec}}}}(n) = \sum\limits_{i = 1}^{{I_{{\rm{Htd}}}}} {\tilde s_{{\rm{imf}}}^i} + \sum\limits_{i = {I_{{\rm{Htd}}}} + 1}^{{I_{{\rm{Ltd}}}}} {s_{{\rm{imf}}}^i} $

(8)

1.1.2 MQRS信息消除和FECG特征周期提取

尽管MECG和FECG的频谱有很大程度的重叠（MQRS的主要频谱在10～30 Hz之间，而FQRS的主要频谱在15～40 Hz之间），但是MECG与FECG的心率不同，两者的QRS波群的基频也不同，在频域上MQRS基频为1.4 Hz左右，而FQRS基频为2.2 Hz左右^[21-22]。这意味着，经EMD分解后，MECG和FECG的主要能量分布在不同的IMF中。事实上，由于MQRS的幅值远大于FQRS，可达5～20倍^[20]，能量主要集中在其方差最大的IMF上；而FQRS的能量通常集中在低1阶的IMF中。图1为DaISy数据库第1导联（Ch1）的EMD分解结果。由图1可以看出：标准差 ${\sigma _{{\rm{imf}}}}$ 的最大值出现在imf4上，其中，FQRS波很小；imf5已不包含FQRS波；相比之下，imf3中的FQRS明显。

1）MQRS信息消除

按第1.1.1节中式（2）计算每个IMF的标准差 $\sigma _{{\rm{imf}}}^i$ ，当 $\sigma _{{\rm{imf}}}^i$ 达到第1个极大值时，将其对应的IMF确定为MQRS为主的IMF，对应阶数记为I。给定一个阈值，通过阈值和峰值剔除方法，将 $s_{{\rm{imf}}}^I$ 中MQRS范围外的区域置1，MQRS内区域置0，从而建立一个剔除函数 $g(n)$ ，然后用 $g(n)$ 与 $s_{{\rm{imf}}}^I$ 及各低阶IMF相乘，消除MQRS信息。具体方法如下：

计算I：

$I = \min \left\{ {\mathop {\arg \max }\limits_i \left\{ {\sigma _{{\rm{imf}}}^i} \right\}} \right\}$

(9)

针对 $s_{{\rm{imf}}}^I$ 定义一个阈值 ${{\varLambda }_{{\rm{T}}{{\rm{d}}_I}}}$ ：

${\varLambda _{{\rm{T}}{{\rm{d}}_I}}} = C\sigma _{{\rm{imf}}}^I$

(10)

式中： $\sigma _{{\rm{imf}}}^I$ 为 $s_{{\rm{imf}}}^I$ 的标准差；C为系数，本文取为2。

剔除函数 $g(n)$ 的建立：

首先，令 $g(n) = s_{{\rm{imf}}}^I(n)$ ，按式（11）循环迭代，直到没有可剔除的极值为止。

$ g(n)=\left\{\! {\begin{array}{l}0,\;\left|g(n)\right|\!<\!{\varLambda }_{{\rm{Td}}_{I}}\cup \left|g(n-1)\right|\!<\!\left|g(n)\right|\!>\!\left|g(n+1)\right|;\\ g(n),\;{\text{其他}}\end{array}} \right.$

(11)

然后，令

$g(n) = \left\{\! \begin{array}{*{20}{c}} 0,\;{g(n) \ne 0;} \\ \!\!1,\;{g(n) = 0} \\ \end{array} \right.$

(12)

将 $g(n)$ 乘以阶数小于等于I的IMF，消除MECG的影响。即：

${\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;s_{{\rm{imf}}}^i(n)} = s_{{\rm{imf}}}^i(n) \cdot g(n)\begin{array}{*{20}{c}} , \end{array}i \le I$

(13)

鉴于阶数大于I的IMF中FECG成分可以忽略，FECG提取和信号重构时，只考虑I及以下阶数的IMF，称I为截止阶数。

2）FECG特征周期提取

由于FECG的基频高于MECG，消除MECG信息后的IMF中， $s_{{\rm{imf}}}^{I - 1}$ 的FQRS最强。对 $s_{{\rm{imf}}}^{I - 1}$ 进行自相关计算，并在1～ ${T_{{\rm{fecg}}}}$ 范围内检测自相关系数的最大峰值点，对应的时延被确定为FECG的特征周期 ${T_\tau }$ 。

设 $R(\tau )$ 为 $s_{{\rm{imf}}}^{I - 1}$ 相对于时延 $\tau $ 的自相关系数，定义为：

$ R(\tau ) = \frac{{\displaystyle\sum\limits_{n = 1}^N {\left( {s_{{\rm{imf}}}^{I - 1}(n) \cdot s_{{\rm{imf}}}^{I - 1}(n + \tau )} \right)} }}{{\sqrt {\displaystyle\sum\limits_{n = 1}^N {{{\left( {s_{{\rm{imf}}}^{I - 1}(n)} \right)}^2}} } \sqrt {\displaystyle\sum\limits_{n = 1}^N {{{\left( {s_{{\rm{imf}}}^{I - 1}(n + \tau )} \right)}^2}} } }}, $

则

${T_\tau } = \mathop {\arg \max }\limits_{1 < \tau < {T_{{\rm{fecg}}}}} \left\{ {R(\tau )} \right\}$

(14)

1.2 QPCE方法

FECG具有准周期特性，在已知FECG的特征周期 ${T_\tau }$ 的情况下，将周期成分提取方法用于提取FECG的IMF是一种比较好的选择。为了区别于标准周期成分提取方法，本文将准周期成分提取方法命名为QPCE（quasi-periodic component extraction）。

假设，由m个混合信号构成的矩阵为 ${{X}} = ({{{x}}_1}, {{{x}}_2}, \cdots ,{{{x}}_m})^{\rm{T}}$ ，其中，x为列向量表示的信号，是由n个相互独立的源信号 ${{{s}}_1},{{{s}}_2}, \cdots ,{{{s}}_n}$ 线性叠加而成的，且s为信号值构成的列向量。设 ${{S}} = {({{{s}}_1},{{{s}}_2}, \cdots ,{{{s}}_n})^{\rm{T}}}$ ，则有：

${{X}} = {{AS}}$

(15)

式中， ${{A}}$ 为 $m \times n$ 阶系数矩阵，称为混合矩阵。需要找到一个正交矩阵 ${{W}} = ({{{w}}_1},{{{w}}_2}, \cdots ,{{{w}}_n})$ （其中w为列向量），使得

${{Y}} = {{{W}}^{\rm{T}}}{{X}}$

(16)

为源信号 ${{S}}$ 的估计。考虑 ${{Y}} = {({{{y}}_1},{{{y}}_2}, \cdots ,{{{y}}_n})^{\rm{T}}}$ 中的第i个信号为目标信号FECG（其中y为列向量），具有准周期特性，其特征周期为 ${T_\tau }$ 。那么有：

${{{y}}_i} = {{w}}_i^{\rm{T}}{{X}}$

(17)

并且

${\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\varepsilon ({{{w}}_i})} = \frac{{\displaystyle\sum\limits_{n = 1}^N {{{({{{y}}_i}(n) - {{{y}}_i}(n + {T_\tau }))}^2}} }}{{\displaystyle\sum\limits_{n = 1}^N {{{({{{y}}_i}(n))}^2}} }}$

(18)

应为最小。

将式（17）代入式（18）可得：

${\;\;\;\;\;\;\;\;\;\;\;\;\;\varepsilon ({{{w}}_i})} = 2 - \frac{{{{w}}_i^{\rm{T}}\left( {{{{R}}_x}({T_\tau }) + {{{R}}^{\rm{T}}_x}{{({T_\tau })}}} \right){{{w}}_i}}}{{{{w}}_i^{\rm{T}}{{{R}}_x}(0){{{w}}_i}}}$

(19)

式中， ${{{R}}_x}({T_\tau }) = {{X}}(n){{X}}^{\rm{T}}{(n + {T_\tau })}$ ， ${{{R}}_x}(0) = {{X}}^{\rm{T}}(n){{X}}{(n)}$ 。

由式（18）可知 $\varepsilon ({{{w}}_i}) > 0$ ，所以 $\varepsilon ({{{w}}_i})$ 最小可以等价于式（19）中的第2项最大。设其最大值为 ${\lambda _{{\rm{max}}}}$ ，并记 ${{A}} = {{{R}}_x}({T_\tau }) + {{{R}}^{\rm{T}}_x}{({T_\tau })}$ 和 ${{B}} = {{{R}}_x}(0)$ ，则有：

$\frac{{{{w}}_i^{\rm{T}}{{A}}{{{w}}_i}}}{{{{w}}_i^{\rm{T}}{{B}}{{{w}}_i}}} = {\lambda _{{\rm{max}}}}$

(20)

或

${{w}}_i^{\rm{T}}{{A}}{{{w}}_i} = {\lambda _{{\rm{max}}}}{{w}}_i^{\rm{T}}{{B}}{{{w}}_i}$

(21)

易知 ${{A}}$ 为对称矩阵， ${{B}}$ 为正定矩阵。用 ${({{w}}_i^{\rm{T}})^{ - 1}}$ 左乘式（21）得：

${{A}}{{{w}}_i} = {\lambda _{{\rm{max}}}}{{B}}{{{w}}_i}$

(22)

那么，式（22）为标准的广义特征值最大值问题。

将矩阵 ${{B}}$ 进行分解，其特征值和特征向量分别记为 ${{\varSigma }} = {\rm{diag(}}{\sigma _1},{\sigma _2}, \cdots ,{\sigma _m})$ （ ${\sigma _i} > 0$ ）和 ${{Q}} = ({{{q}}_1},{{{q}}_2}, \cdots ,{{{q}}_m})$ （其中q为列向量），则有：

${{B}} = {{Q\varSigma }}{{{Q}}^{\rm{T}}} = {{{P}}^{\rm{T}}}{{P}}$

(23)

式中， ${{P}} = {{{\varSigma }}^{1/2}}{{{Q}}^{\rm{T}}}$ 和 ${{{P}}^{ - 1}} = {{Q}}{{{\varSigma }}^{ - 1/2}}$ 为非奇异变换。

将式（23）代入式（22），并令 ${{{u}}_i} = {{P}}{{{w}}_i}$ ， ${{{w}}_i} = {{{P}}^{ - 1}}{{{u}}_i}$ 和对称矩阵 ${{S}} = {({{{P}}^{ - 1}})^{\rm{T}}}{{A}}{{{P}}^{ - 1}}$ 。则式（22）又可化为标准的特征值最大值问题。

${{S}}{{{u}}_i} = {\lambda _{{\rm{max}}}}{{{u}}_i}$

(24)

求得 ${{S}}$ 的最大特征值 ${\lambda _{{\rm{max}}}}$ 和对应的特征向量 ${{{u}}_B}$ ，则最佳提取向量 ${{{w}}_B}$ 为：

${{{w}}_B} = {{{P}}^{ - 1}}{{{u}}_B}$

(25)

1.3 IMF提取与FECG重构

设腹壁导联信号组成的矩阵 ${{X}} = {({{{x}}_1},{{{x}}_2}, \cdots ,{{{x}}_m})^{\rm{T}}}$ ，其中，每一列代表一个腹壁导联信号，m为导联个数。IMF提取与FECG重构步骤如下：

1）用第1.1节中的方法对各个导联信号进行处理，得出消除MQRS影响的IMF、IMF截止阶数I和FECG特征周期 ${T_\tau }$ 。若各导联的I不同，则取最大值。 ${T_\tau }$ 则由各导联提取的特征周期综合分析确定（或可通过其他方法获得）。

2）首先，将各导联同阶的IMF组成向量组如下：

${{s}}_{{\rm{IMF}}}^i(n) = {[\!\!\begin{array}{*{20}{c}} {s_{{\rm{im}}{{\rm{f}}^1}}^i(n)}&{s_{{\rm{im}}{{\rm{f}}^2}}^i(n)}& \cdots &{s_{{\rm{im}}{{\rm{f}}^m}}^i(n)} \end{array}\!\!]^{\rm{T}}}$

(26)

式中：imf的上角标m为导联数；s的上角标i为第i阶， $i = 1,2, \cdots ,I$ 。

然后，对式（26）中的每组信号，用第1.2节中的QPCE方法，按准周期 ${T_\tau }$ 提取FECG的IMF，记为 $s_{{\rm{Fimf}}}^i$ ，并重构FECG信号：

${s_{{\rm{fecg}}}}(n) = \sum\limits_{i = 1}^I {s_{{\rm{Fimf}}}^i(n)} $

(27)

计算 ${s_{{\rm{fecg}}}}(n)$ 的标准差为 ${\sigma _{{\rm{fecg}}}}$ ：

${\;\;\;\;\;\;\;\;\;\;\;\;\;\;\sigma _{{\rm{fecg}}}} = \sqrt {\frac{1}{{N - 1}}\sum\limits_{n = 1}^N {s_{{\rm{fecg}}}^2(n)} } $

(28)

标准化式（27）得：

${s_{{\rm{fecg}}}}(n) = {s_{{\rm{fecg}}}}(n)/{\sigma _{{\rm{fecg}}}}$

(29)

2 性能评价方法 2.1 仿真信号的性能评价方法

关于FECG提取效果的评价，到目前为止尚没有统一的方法。对于仿真信号，由于已有干净的FECG信号，通常采用均方根误差（root mean square error，RMSE）和相关系数（correlation coefficient，CC）作为评价指标。用字母H统一表示评价指标相关变量。

${\;\;\;\;\;\;\;\;\;\;\;H_{{\rm{RMSE}}}} = \sqrt {\frac{1}{N}\sum\limits_{n = 1}^N {{{(f(n) - F(n))}^2}} } $

(30)

${\;\;\;\;\;\;H_{{\rm{CC}}}} \!=\! \frac{{\displaystyle\sum\limits_{n = 1}^N {\left[ {(f(n) \!-\! \overline {f(n)} )(F(n) \!-\! \overline {F(n)} )} \right]} }}{{\sqrt {\displaystyle\sum\limits_{n = 1}^N {{{(f(n) \!-\! \overline {f(n)} )}^2}} \displaystyle\sum\limits_{n = 1}^N {{{(F(n) \!-\! \overline {F(n)} )}^2}} } }}$

(31)

式中， $f(n)$ 为仿真FECG信号， $F(n)$ 为从混合信号中提取的FECG信号， $\overline {f(n)} $ 和 $\overline {F(n)} $ 为相应信号的均值。

2.2 真实信号的性能评价方法

对于临床获得的真实腹壁混合心电信号，不可能有干净的FECG信号作为参考，尽管ADFECGDB数据库中已有直接测量头皮的FECG信号，但是仍然含有噪声，不能作为干净的FECG信号使用，基于头皮FECG计算的RMSE和CC不能真实反映实际情况。为此，本文参考相关文献[23-24]采取统计指标和可视化方法进行评价。

首先，确定成功提取的定性标准，成功提取的信号需满足以下3个条件：

1）MECG能够有效被分离；

2）FECG的R波能够在波形中清晰显现；

3）从波形上能够得到准确的胎儿心率。

对于满足上述3个条件的成功提取信号，借鉴文献[23]的R峰探测指标，定义R峰提取灵敏度Se、阳性预测值PPV和测量准确度F₁指标表达式如下：

${H_{{\rm{Se}}}} = \frac{{{N_{{\rm{TP}}}}}}{{{N_{{\rm{TP}}}} + {N_{{\rm{FN}}}}}}$

(32)

${H_{{\rm{PPV}}}} = \frac{{{N_{{\rm{TP}}}}}}{{{N_{{\rm{TP}}}}{\rm{ + }}{N_{{\rm{FP}}}}}}$

(33)

${H_{{{\rm{F}}_{\rm{1}}}}} = 2\frac{{{H_{{\rm{PPV}}}} \cdot {H_{{\rm{Se}}}}}}{{{H_{{\rm{PPV}}}} + {H_{{\rm{Se}}}}}} = \frac{{2{N_{{\rm{TP}}}}}}{{2{N_{{\rm{TP}}}} + {N_{{\rm{FN}}}} + {N_{{\rm{FP}}}}}}$

(34)

式中， ${N_{{\rm{TP}}}}$ 、 ${N_{{\rm{FP}}}}$ 、 ${N_{{\rm{FN}}}}$ 分别为真阳性即正确的FQRS的数目TP、假阳性即错误的FQRS的数目FP、假阴性即丢失的FQRS数目FN。

对应的心拍胎心率FHR（fetal heart rate）由式（35）计算：

${H_{{\rm{FHR}}}} = \frac{{60{f_{{\rm{sample}}}}}}{{{N_{{\rm{RR}}}}}}$

(35)

式中， ${N_{{\rm{RR}}}}$ 为相邻R峰（包括真阳性和假阳性）之间的采样点数， ${f_{{\rm{sample}}}}$ 为采样频率。

平均胎心率由式（35）获得的每个心拍胎心率平均计算：

${\overline H_{{\rm{FHR}}}} = \frac{1}{K}\sum\limits_{i = 1}^K {H_{_{{\rm{FHR}}}}^i} $

(36)

平均心率误差MFHRE定义如下：

${\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;H_{{\rm{MFHRE}}}} = \frac{{\left| {{{\overline H}_{{\rm{FH}}{{\rm{R}}_{\rm{D}}}}} - {{\overline H}_{{\rm{FH}}{{\rm{R}}_{\rm{E}}}}}} \right|}}{{{{\overline H}_{{\rm{FH}}{{\rm{R}}_{\rm{E}}}}}}}$

(37)

式中， ${\overline H_{{\rm{FH}}{{\rm{R}}_{\rm{D}}}}}$ 为真实心电数据库中专家标注的胎心率平均值， ${\overline H_{{\rm{FH}}{{\rm{R}}_{\rm{E}}}}}$ 为提取的FECG的胎心率平均值（消除丢失R峰的影响）。

3 实验结果及讨论 3.1 DaISy数据库实验结果

DaISy数据库^[25]由比利时学者Lathauwer提供，包含一组8导联数据。其中，前5导联信号来自于孕妇腹部电极，后3导联信号来自于孕妇胸部电极。采样频率为250 Hz，采样时间为10 s。图2为1～5导联腹壁信号波形。

将主成分分析（PCA）、准周期成分提取（QPCE）、独立成分分析（ICA）3种方法与本文提出的EMD–QPCE方法进行对比，以评价本文方法的好坏。其中，由于ICA方法提取的FECG在大小和方向上的不确定性，在进行比较之前，需要将各种方法的提取结果处理为单位方差的标准化形式。使用上述方法对DaISy数据的前5导联进行多通道FECG提取，结果绘于图3。由图3可知，DaISy数据库的信噪比较高，各种方法提取的效果都很好，而本文的EMD–QPCE方法提取结果（图3（d））的噪声最小。

3.2 ADFECGDB数据库实验结果

ADFECGDB数据库^[26]共有r01、r04、r07、r08和r10共5个数据集，每个数据集包含1个直接头皮测量信号和4个腹壁测量信号。该数据库是波兰西里西亚医科大学妇产科系2012年发布的。采样频率1 000 Hz，采样时长5 min，本文选择长度10 s的信号进行FECG提取实验。

图4～8分别给出了本文提出的EMD–QPCE方法与现有的主成分分析（PCA）、准周期成分提取（QPCE）、独立成分分析（ICA）方法对ADFECGDB数据库中每个数据集的FENG提取结果。图4～8中：为了更好地反映本文方法的提取准确性，在（a）和（b）中标记出了FECG的R峰位置，（a）中的标记来源于原始数据的专家标记，（b）中标记来源于EMD–QPCE提取结果。由图4～8可知：对于r01和r08，上述4种方法都能提取出FECG信号，但是，PCA、QPCE、ICA方法都有很大程度的MECG残留，本文EMD–QPCE方法的MECG残留较少。对于r04、r07和r10，只有EMD–QPCE方法成功提取出了FECG信号。

由此表明：本文EMD–QPCE方法能够消除MECG的相关性影响，明显提高FECG的成功提取率。然而，EMD–QPCE方法在FQRS与MQRS重叠之处还是会产生FQRS的丢失。EMD–QPCE方法提取的R峰更清晰，具有较高的R峰位置准确性，因此不影响胎心率的计算。

图4～8中，对EMD–QPCE方法提取的FECG的R峰位置进行了标记。通过与相应的专家标记结果对比，R峰统计和平均心率计算结果分别列于表1和2中。表1中，前5行是对每个数据集分别统计TP、FP、FN、Se、PPV、F₁的结果，最后一行是由5个数据集进行总体统计的结果。表2中，最后一行为根据5个数据集分别提取的平均胎心率再计算得到的平均值。由图4～8和表1、2的结果可知：对ADFECGDB数据库全部数据集，EMD–QPCE方法的FECG成功提取率达到100%，而其他3种方法的成功提取率仅为40%。EMD–QPCE方法的总体统计结果中，灵敏度Se为96.1%，阳性预测值PPV为99%，准确度F₁为97.5%；平均胎心率误差MFHRE的5组平均计算结果为0.06%，最大不超过0.23%。充分证明了本文方法的有效性和优越性。

3.3 Challenge 2013数据库实验结果

Challenge 2013数据库^[27]分为两个部分，即训练数据库Set–A和测试数据库Set–B。数据库中的每组数据都只有4个腹壁测量信号，采样频率1 000 Hz，采样时长5 min。Set–A包含有专家对FECG的R峰标注，以便对提取结果进行检验；Set–B中没有标注。本文分别在这两个数据库中选取了前20组数据，用PCA、QPCE、ICA和本文EMD–QPCE方法进行FECG提取实验。

表3为Set–A中20组数据的PCA、QPCE、ICA和本文EMD–QPCE方法成功提取结果的统计。

由表3可知，EMD–QPCE方法的成功提取率为60%，远高于PCA方法的25%和QPCE、ICA方法的30%。

表4为Set–B中20组数据的PCA、QPCE、ICA和本文EMD–QPCE方法成功提取结果统计。

由表4可知，EMD–QPCE方法的成功提取率为50%，ICA方法为20%，PCA、QPCE方法的成功提取率仅为15%。事实上，以往的提取方法对Challenge 2013数据库的成功提取率都不高。本文选择Challenge 2013数据库，主要源于Set–A具有专家标注结果，便于对实验结果进行检验。表3和4结果表明，EMD–QPCE方法的成功提取率明显高于传统的PCA、QPCE和ICA方法。

综上，值得指出的是：在选择的20组Set–A数据和20组Set–B数据中，只要PCA、QPCE和ICA 3种方法中有一个方法能够成功提取FECG时，EMD–QPCE方法一定能成功提取出FECG，如Set–A04（图9）和Set–B08。对于3种方法均不能提取FECG的情况，只要信号中的FQRS可分辨，EMD–QPCE方法仍能提取FECG，如Set–A19（图10）。EMD–QPCE方法未能提取FECG是因为数据集的各个通道均无可分辨的FECG信号，如Set–B04（图11）。

图9为Set–A04中EMD–QPCE、PCA、ICA、QPCE提取的FECG波形。其中，图9（a）为Set–A04第4通道的原始数据波形和专家标注结果，图9（b）标出了EMD–QPCE提取FECG的R峰位置。

由图9可以看出：EMD–QPCE提取的FECG更清晰，R峰位置与专家标注结果一致。ICA和QPCE方法虽然也能提取出FECG，但是，噪声比较大，甚至会影响FECG的R峰识别。

图10为Set–A19中EMD–QPCE、PCA、ICA、QPCE方法的提取结果。其中，图10（a）给出的是第2通道的原始数据波形和专家标注结果。

由图10可以看出：EMD–QPCE提取的FECG波形中R峰准确清晰。而其他方法都未能提取出FECG，然而，其波形中可以看到存在明显的FECG，这说明在Set–A19中，FECG和MECG可能具有较强相关性，导致PCA、ICA、QPCE等盲源分离算法不能有效地将其分离。

图11为Set–B04的原始信号波形和去除基线漂移后的波形。

由图11可知，Set–B04的各个通道均观察不到有效的FECG信息。对于这类信号，EMD–QPCE方法也不能提取出FECG。

由于Set–B的数据中没有FQRS的标注，因此仅对Set–A中成功提取的FECG进行总体统计。在本文选取的20组Set–A数据记录中，由EMD–QPCE方法能成功提取FECG的数据记录有：Set–A01、Set–A03、Set–A04、Set–A05、Set–A08、Set–A10、Set–A12、Set–A13、Set–A15、Set–A17、Set–A19和Set–A20共计12个数据记录。12个数据记录中，专家标记的FECG的R峰总数为263个；EMD–QPCE提取的FECG的R峰总数为243个，其中，准确提取数TP为240，错误提取数FP为3，与标记数对比丢失的R峰数FN为20。计算结果表明，EMD–QPCE的总体灵敏度Se为92.31%，阳性预测值PPV为98.77%，准确度F₁为95.43%，平均胎儿心率误差MFHRE仅为0.595%。由此充分反映出EMD–QPCE方法具有较高的提取精度。

4 结　论

针对多通道腹壁胎儿心电信号的提取，本文提出了EMD–QPCE方法。该方法在去噪和基线漂移滤除后，进行EMD分解，基于FQRS和MQRS的基频差异，在IMF层面上去除MQRS的干扰和相关性，采取QPCE方法对IMF组合进行FECG信息提取，然后依据提取的各阶IMF信息重构FECG。通过对DaISy数据库、ADFECGDB数据库和Challenge 2013数据库的真实胎儿心电信号进行不同算法的提取实验，证明EMD–QPCE方法的成功提取率显著高于传统的PCA、QPCE和ICA方法。对于PCA、QPCE和ICA方法能够成功提取出的FECG真实信号，EMD–QPCE也能够成功提取并具有更好的质量。对于PCA、QPCE和ICA方法不能成功提取，但FQRS可分辨的信号，EMD–QPCE也能够成功提取。这证明了EMD消除MECG干扰和相关性的方法是有效的。统计结果表明：EMD–QPCE方法的FECG提取灵敏度Se为92.31%，阳性预测值PPV为98.77%，准确度F₁为95.43%，平均胎儿心率误差MFHRE为0.595%，具有非常好的准确度和精度，以及应用的潜力和价值。

值得指出的是，EMD–QPCE方法在FQRS与MQRS有严重重叠之处，会产生FQRS的丢失，而提取出的其他R峰都很准确和清晰，能够正确计算胎儿心率。关于R峰丢失问题，是在消除MQRS影响时引起的，有待进一步研究和解决。