在产科的临床医疗领域，对孕妇进行孕期监测被认为是检测胎儿是否健康发育的最为有效手段之一，因此，产生了多种孕期监测的方法.目前应用最为广泛的是基于超声多普勒原理的胎心率(Fetal Heart Rate, FHR)监护仪^[1]，该方法通过探头主动发射超声信号作用于胎儿，再根据胎儿心脏的跳动对超声信号的不同反射可以获得胎儿的瞬时心率，从而确定胎儿的健康状况.但也因此被认为是一种“入侵式”的监护方法，基于安全的理由不能长时间使用.而基于胎儿心电(Fetal Electrocardiogram, FECG)监护的方法^[2-3]通过母亲腹部的表面电极被动采集胎儿的心电信号^[4-5]，与传统方法相比，它更加灵敏以及具有更好的安全性，能进行长时间监护，同时采集到的胎儿心电信号能更全面反映胎儿的健康状况，所以可以有效地发现多种胎儿疾病，包括胎儿缺氧、脐带缠绕等.

但是在实际的应用中，由于存在母体心电(Maternal Electrocardiogram, MECG)、肌电信号、体表子宫电信号以及各种噪声的干扰^[6]，并且胎儿心电信号本身又极其微弱，幅值仅为母体心电信号的10%^[7]，通过孕妇腹部采集到的胎儿心电信号存在明显的观测困难.因此，如何从这样的混叠信号中提取出干净的胎儿心电信号是产科临床医学中的一大难题.

自20世纪90年代以来，信号处理中的盲分离(Blind Source Separation, BSS)技术^[8]得到迅速发展，该技术能在对源信号缺乏先验知识且混叠未知的情况下，根据信号的独立性、稀疏性等统计特性，通过运算处理，最终分离出各路源信号.基于胎心监护中混叠的母体心电与胎儿心电具有独立性，且在大部分时段满足相对稀疏的特征，所以，盲分离是提取胎儿心电的一种有效方法^[9].

盲分离的信号模型主要分为线性混叠模型和卷积混叠模型.其中，线性混叠模型较为简单，是指多路源信号经过线性变化组合而成.目前，典型的盲分离及独立分量分析主要就是针对线性混叠模型开展的，一般假设源信号满足线性混叠模型，这类问题也称为常规的盲分离问题.然而，在实际临床产科应用中，由于人体复杂的信号采集环境且孕妇胎位因人而异，通过电极获得的观测信号往往会出现不同程度的畸变或者非线性混叠, 这直接导致信号模型不满足线性混叠.简单的处理是调整电极的位置，使得观测信号近似满足线性混叠模型，或者重新选取更加满足线性混叠模型的观测信号.为此，本文提出了一种高效、准确的基于时域稀疏信号相关性计算，以判断采集到的孕妇腹部心电是否满足线性混叠模型的盲分离分析方法.同时，针对心电信号的采集过程中可能遇到各种干扰信号而导致判别方法的结果出现误差甚至与真实情况相反的情况，本文认为在进行相关分析前需要对采集到的原始信号进行初步的滤波处理，从而保证本文提出的判别方法更加准确可靠.

1 信号相关性的判别

主要针对满足时域稀疏性的源信号(在本文中指的就是胎儿心电与母体心电)进行线性混叠度分析，通过计算源信号两两间的相关性，从而得到一组更加符合线性混叠盲分离模型的观测信号，为线性混叠的盲源分离问题提供一种简单有效的验证及判据, 主要步骤如下.

1.1 获取稀疏区域

通过现有的母体心电以及胎儿心电的监测及定位方法^[10-11]，分别得到观测信号中母体心电与胎儿心电R波位置，如图 1所示.其中，R_m与R_f分别是母体心电与胎儿心电的R波位置.再根据成人与胎儿心电信号的长度信息估计出母体心电与胎儿心电的PQRS波段所占的时域范围，从而找到没有发生混叠的时段.分别截取出各路观测信号中没有出现混叠的母体心电与胎儿心电相对稀疏的区域，组成新的信号k_i(n), n=1, 2, …, M; i=1, 2, …, N.其中M为由稀疏区域组成的新信号的长度，N表示一共有N路观察信号.

1.2 预处理

在度量信号中心电成分的相关性前，需要对新截取的稀疏性信号k_i(n)做预处理，其中包括了去均值以及归一化处理，处理的公式分别为

$ {{\bar k}_i}\left( n \right) = {k_i}\left( n \right) - \left( {1/M} \right)\sum\nolimits_{j = 1}^M {{k_i}\left( j \right)} ,\left( {n = 1,2, \cdots ,M} \right), $

(1)

$ {K_i}\left( n \right) = {{\bar k}_i}\left( n \right)/\sqrt {\sum\nolimits_{j = 1}^M {{{\left[ {{{\bar k}_i}\left( j \right)} \right]}^2}} } ,\left( {n = 1,2, \cdots ,M} \right), $

(2)

其中, ${\bar k_i}\left( n \right) $表示去均值后的信号，K_i(n)表示去均值且归一化后的信号.

1.3 相关性度量

得到预处理后的信号K_i(n)后，可以通过公式(3)计算出观测信号中心电信号的线性相关度为

$ C\left( {i,j} \right) = \sum\nolimits_{n = 1}^M {{K_i}\left( n \right){K_j}\left( n \right)} ,\left( {i < j;n = 1,2, \cdots ,M} \right). $

(3)

利用穷举法以及求均值的方法，通过式(4)计算所有的N路观测到的孕妇腹部混叠心电信号中任意组合的J路信号间的平均相关性.

$ {C_{{\rm{ave}}}} = \left( {1/C_N^J} \right)\sum\nolimits_{i = 1,j = 1}^J {C\left( {i,j} \right)} ,\left( {i < j} \right), $

(4)

其中，C_N^J表示从N路信号中任意取J路不同组合信号的所有可能情况的数量.

最后，利用阈值判断的方法，就可以实现对组合所得的J路观测信号进行是否符合线性混叠模型的判断.判断的规则为

C_ave>T：选取的多路观测信号满足线性混叠模型;

C_ave≤T：选取的多路观测信号不满足线性混叠模型.

其中，判断阈值T可以根据信号的实际情况进行选取.当相关性测量值C_ave大于阈值T时，表明该J路观测信号的组合满足线性混叠模型，可以利用常规盲分离算法进行胎心电信号的提取; 否则，需要调整孕妇腹部的采集电极或者选择另外一种组合的观测信号，直到相关性测量值满足判断的规则.

2 观测信号去噪

通过上述的方法以及文献[12]中的基于相关性的判断方法，可以快速有效地判断观测信号中混叠的胎儿心电与母亲心电是否满足常规的线性叠加盲分离模型.然而，原始的心电采集信号往往淹没在各种噪声中.而这些噪声的存在往往破坏了信号的稀疏性以及线性叠加的特性.因此，在不考虑噪声干扰的情况下，对信号进行相关性的判别所得到的相关性测量值将与真实的心电观测信号间的相关度存在较大的偏差，甚至可能出现判断结果与真实的信号相关性相反的情况.此时，仅仅利用相关性判别无法准确知道观测信号中的混叠心电信号是否满足线性叠加盲分离模型.

为了去除噪声信号的影响，需要对采集到的原始信号作去噪处理.在心电采集过程中常见的非心电噪声有肌电信号、工频信号、基线漂移等等，因此针对这些不同的噪声，研究人员提出了各种不同的去噪方法^[13-15].在此，本文针对心电信号采集中最常见的噪声——工频干扰进行探讨并提出相关的去噪方法.

工频干扰主要指的是50 Hz及其谐波(100 Hz、150 Hz, …)的干扰.基于这些干扰信号的频域特征，本文将利用梳状滤波器^[16]完成去除工频干扰的目的.令f₁=50 Hz，工频谐波可以通过式(5)和(6)表示：

$ {\omega _1} = \left( {2{\rm{ \mathsf{ π} }}{f_1}} \right)/{f_s}\left[ {{\rm{radians}}/{\rm{sample}}} \right], $

(5)

$ {\omega _i} = i{\omega _1},i = 1,2,3, \cdots ,M. $

(6)

其中，f_s为信号的采样率; M为根据奈奎斯特采样定理采样信号可能存在的不同频率谐波的个数.而梳状滤波器的传递函数如式(7)所示：

$ \begin{array}{l} \;\;\;\;H\left( z \right) = \frac{{N\left( z \right)}}{{N\left( {{\rho ^{ - 1}}z} \right)}} = \\ \frac{{1 + {b_1}{z^{ - 1}} + {b_2}{z^{ - 2}} + \cdots + {b_M}{z^{ - M}}}}{{1 + \rho {b_1}{z^{ - 1}} + {\rho ^2}{b_2}{z^{ - 2}} + \cdots + {\rho ^M}{b_M}{z^{ - M}}}}, \end{array} $

(7)

其中，z= e^±jω₀，N(z)为多项式表达式，其中0 < ρ < 1.因此，通过梳状滤波器可以有效地去除心电采集中的工频干扰，使得采集到的信号更加接近真实的心电信号.

3 基于FastICA的胎心电提取技术

独立分量分析方法(Independent Component Analysis，ICA)^[17]是一种基于源信号统计独立性的分析方法，目前该理论发展较为成熟，在语音信号分离、特征提取、生物电信号处理、CDMA通信系统以及金融数据分析等信号处理领域得到广泛的应用^[18-19].其数学模型表示为

$ {\mathit{\boldsymbol{x}}_i}\left( t \right) = \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{s}}_i}\left( t \right), $

(8)

其中，$\mathit{\boldsymbol{s}}{_i}\left( t \right) = {[{\rm{ }}s{_1}\left( t \right), {\rm{ }}s{_2}\left( t \right), \ldots , {\rm{ }}s{_N}\left( t \right)]^{\rm{T}}} $为由N路源信号组成的矩阵，A为N×N维未知的混合矩阵，$ \mathit{\boldsymbol{x}}{_i}\left( t \right) = {\left[ {{x_1}\left( t \right), {x_2}\left( t \right), \ldots , {x_N}\left( t \right)} \right]^{\rm{T}}}$为由N路观测信号组成的矩阵.在本文中，观测信号指的是从孕妇腹部采集到的心电信号，其中包含了胎儿心电、母体心电以及各种采集中遇到的干扰信号.在此，本文所感兴趣的源信号就是胎儿的心电信号.对源信号的估计为

$ {\mathit{\boldsymbol{z}}_i}\left( t \right) = \mathit{\boldsymbol{W}}{\mathit{\boldsymbol{x}}_i}\left( t \right) = \mathit{\boldsymbol{WA}}{\mathit{\boldsymbol{s}}_i}\left( t \right) $

(9)

其中，W为N×N维的解混矩阵，$ \mathit{\boldsymbol{z}}{_i}\left( t \right) = {[{\rm{ }}z{_1}\left( t \right), {\rm{ }}z{_2}\left( t \right), \ldots , {\rm{ }}z{_N}\left( t \right)]^{\rm{T}}}$表示对源信号的估计.因此ICA问题就在于如何获得矩阵W，使之满足W = A^-1.

ICA的优化算法大致可分为两类，即批处理算法(Batch processing)与自适应处理算法(Adaptive processing).本文采用的FastICA算法就是一种对于胎心提取有效的批处理算法.

FastICA算法，又称固定点(Fixed-Point)算法，是一种快速寻优迭代算法，能大大减少ICA处理的运算量并具有更快的收敛速度.算法的过程大致可以分为两个部分，包括预处理和提取独立分量.

1) 预处理

通过对输入信号进行预处理，可以有效地减少算法的运算量.针对基于FastICA的胎儿心电信号的提取，所涉及的预处理包括信号的去均值和白化处理分别为

$ \mathit{\boldsymbol{\bar X}} = \mathit{\boldsymbol{x}} - \mathit{\boldsymbol{E}}\left[ {{x}} \right], $

(10)

$ \mathit{\boldsymbol{\tilde X}} = \mathit{\boldsymbol{E}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}^{ - 1/2}}{\mathit{\boldsymbol{E}}^{\rm{T}}}\mathit{\boldsymbol{\bar X}}. $

(11)

其中，x表示观测信号，E[x]表示观测信号的数学期望值，$\mathit{\boldsymbol{\bar X}} $表示去均值后的观测信号，Λ表示N×N维的对角矩阵，E表示正交矩阵，$ \mathit{\boldsymbol{\tilde X}}$表示白化处理后的观测信号，即完成预处理后的信号.

2) 提取独立分量

以负熵最大化作为算法的搜索方向，最终迭代找到一个矩阵用以提取出感兴趣的独立成分.根据牛顿迭代理论，FastICA算法的迭代公式表示如下：

$ {\mathit{\boldsymbol{w}}^ * } = \mathit{\boldsymbol{E}}\left\{ {\mathit{\boldsymbol{z}}g\left( {{\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{z}}} \right)} \right\} - \mathit{\boldsymbol{E}}\left\{ {g'\left( {{\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{z}}} \right)} \right\}\mathit{\boldsymbol{w}}, $

(12)

$ \mathit{\boldsymbol{w}} = {\mathit{\boldsymbol{w}}^ * }/{\mathit{\boldsymbol{w}}^ * }. $

(13)

其中，w为迭代的初始权值，g(·)为非线性函数，z为处理后的信号，w^Tz为z在w^T方向的投影.循环迭代，直到收敛.对于提取多个独立分量，可以通过重复上述的步骤实现，直到提取完所有感兴趣的独立分量为止.当假设已经提取到了P-1的独立分量时，对于权值w需要进行处理，如式(14)、(15)所示.

$ {\mathit{\boldsymbol{w}}_p} = {\mathit{\boldsymbol{w}}_p} - \sum\nolimits_{j = 1}^{p - 1} {\mathit{\boldsymbol{w}}_p^{\rm{T}}{\mathit{\boldsymbol{w}}_j}{\mathit{\boldsymbol{w}}_j}} , $

(14)

$ {\mathit{\boldsymbol{w}}_p} = {\mathit{\boldsymbol{w}}_p}/{\mathit{\boldsymbol{w}}_p}. $

(15)

综上所述，本文提出的判别方法可以根据信号间的相关度测量情况，有效地判断观测信号是否满足常规的线性叠加盲分离模型，从而获取到更加适用常规盲分离处理方法的观测信号，大大降低信号处理的难度，使得通过常规的FastICA处理就可以实现胎儿心电的提取.同时，考虑到心电采集过程中会出现噪声干扰，使得信号的相关性被破坏的情况，本文将针对其中最为常见的工频干扰进行相应的滤波处理.整个过程的示意图如图 2所示.

4 仿真分析与讨论 4.1 数据来源

仿真采用De Moor^[20]记录的一组8通道的孕妇心电信号，如图 3所示.将编号为x₁~ x₈，信号对应的采样率为250 Hz，信号长度为2 500，采样时长为10 s.其中x₁~ x₅为孕妇腹部混叠心电，包含母体心电以及胎儿心电成分; x₆~ x₈为孕妇胸部的母体心电.在本文的论述与处理中只涉及腹部混叠心电，因此，下面将针对x₁~x₅这5路信号进行讨论.

4.2 Matlab仿真实验

仿真实验的过程如图 4所示.由于实验中采用的信号本身不存在明显的干扰成分，能较真实地反映心电信号的叠加情况.因此，首先将不考虑实际采集过程中的噪声因素，直接利用本文提出的判别分析方法对5路混叠的观测信号进行相关度的分析，从而得到信号中真实的胎儿心电与母体心电的相关度测量值，并利用判断规则找到其中最符合线性叠加盲分离模型的多组信号.其次，本文将在观测信号中引入工频干扰，以模拟临床采集环境中观测信号的真实情况，分别在忽略噪声以及去除噪声的情况下利用判别分析方法分别进行计算.将所得到的测量值与真实值进行比较，从而确定噪声对观测信号中心电成分的相关性的影响情况.最后，根据判别方法的结果，以两路信号为一组，找到其中多组不同程度地满足线性叠加模型的信号组，并分别对它们进行FastICA处理，试图分离提取出干净的胎儿心电并对比提取的效果，由此来验证判别方法的准确性.

通过现有的R波定位技术，找到混叠信号中母体心电与胎儿心电的R波位置信息，如表 1所示.可以发现心电信号具有一定的周期性且胎儿心电频率要高于母体心电.利用PQRS波段所在位置的估计信息，找到母体心电与胎儿心电相对稀疏即没有发生混叠的区域，如图 5所示.其中，PQRS_F与PQRS_M分别表示观测信号中母体心电与胎儿心电的PQRS波段所在的位置.通过母体心电以及胎儿心电的PQRS波段所在位置进行比较，找到信号相对稀疏区域，分别是806~981、1 360~1 544.接着，对各路信号中对应的稀疏区域进行截取组合，可以得到5路新信号k₁~k₅，如图 6所示.为了进行信号相关性的判定，利用式(1)~(2)对信号k₁~ k₅进行去均值以及归一化的预处理，处理后的信号如图 7所示.最后，利用穷举法找出信号两两的任意组合，再利用式(3)计算信号间的相关性测量值，结果如表 2所示.其中，C(i, j)表示信号x_i与 x_j间的线性相关度.

但是，在实际的心电采集中，常常会遇到各种干扰，而这些干扰可能会影响到相关性判断的结果.下面将模拟存在工频干扰的情况：在原始观测信号x₁~ x₅中，分别引入工频干扰信号，并将它们视为在实际的采集过程中所采集到的带有噪声干扰的心电信号x_n1~ x_n5，如图 8所示.图 9给出了对其中一路含噪信号x_n1作时频域分析的结果，其中可以看到信号中存在明显的工频干扰，即频率为50 Hz及100 Hz的干扰成分.在不对x_n1~ x_n5作任何去噪处理的情况下，同样截取信号的稀疏区域，组成新的信号k_n1~ k_n5，并利用文章中提到的相关性判断方法进行测量，测量结果如表 3所示.从表 3中可以看到，在引入噪声后，如果不作去噪处理直接进行相关性计算，其结果与信号的真实值存在着明显的误差.

考虑到引入的是工频干扰，本文利用梳状滤波器进行去噪处理.通过仿真最终得到滤波器的频率响应曲线如图 10所示，其中可以看到50Hz对应的ω₁=(2π×50Hz)/250Hz=0.4π(radians/sample)以及100Hz对应的ω₂=(2π×100Hz)/250Hz=0.8π(radians/sample)存在明显的抑制.将信号x_n1~ x_n5分别都利用滤波器进行去噪，可以得到滤波后的信号x_f1~x_f5，波形如图 11所示.图 12给出了对其中一路去噪信号x_f1作时频域分析的结果，其中发现滤波处理后的信号中的工频干扰得到了明显的抑制.最后再次利用判断方法计算去噪后的信号的相关性，所得到的相关性测量值如表 4所示.

为了对比去噪前后的结果，本文将按照公式(16)、(17)分别计算表格中有关测量值的均方误差：

$ {\rm{ME}}{{\rm{S}}_n} = \left( {1/C_5^2} \right)\sum\nolimits_{i = 1,j = 1}^5 {{{\left( {C\left( {i,j} \right) - {C_n}\left( {i,j} \right)} \right)}^2}} ,\left( {i < j} \right), $

(16)

$ {\rm{ME}}{{\rm{S}}_f} = \left( {1/C_5^2} \right)\sum\nolimits_{i = 1,j = 1}^5 {{{\left( {C\left( {i,j} \right) - {C_f}\left( {i,j} \right)} \right)}^2}} ,\left( {i < j} \right). $

(17)

其中，MES_n表示C与C_n间的均方误差; MES_f表示C与C_f间的均方误差; C₅²表示5个信号中取出2个信号的所有组合数.计算结果为MES_n=5.863 3%，MES_f=0.428 7%.因此，可以看出噪声环境下极大地破坏了信号真实的相关性，导致测量值无法真实描述信号间的相关度，更无法得知信号是否满足线性叠加的盲分离模型; 而去噪后的测量值与真实信号x₁~x₅的相关性情况基本一致，说明了对采集到信号进行响应的去噪处理是非常有必要的，能够有效地还原信号的真实相关性，为判定观测信号是否满足线性叠加盲分离模型提供有效的保证.

最后，为了验证本判别方法的有效性，本文从中抽取4组相关性测量值由小到大排列的信号组合，并对其进行FastICA处理.4组信号组合分别是组合1：[x₃, x₄]，组合2：[x₁, x₄]，组合：3[x₁, x₃]，组合4：[x₃, x₅]，其相关性测量值分别为C(3, 4)=0.097 0; C(1, 4)=0.313 0;C(1, 3)=0.588 6;C(3, 5)=0.949 6，由表 2可知组合4的相关性最高.4组信号的输出结果如图 13所示.其中，根据胎儿心率要高于成人心率的特点，可以看到输出结果中MECG与FECG分别代表了母体心电以及胎儿心电，而组合4的分离效果显然是最佳的，得到了较为干净的胎儿心电.

4.3 结果分析

从仿真结果可以看到，在观测信号本身不含太多非心电噪声的情况下，通过基于心电信号相关性的判别分析方法，可以准确地反映信号中的胎儿心电与母体心电满足线性叠加盲分离模型的程度.当测量值较高且高于某一阈值时，通过常规的盲分离算法，如文中采用的FastICA算法就可以实现胎儿心电信号的分离与提取，其提取效果如仿真实验中组合4的结果所示; 反之，当组合的相关度测量值较小时，比如组合1，其分离结果没能实现胎儿心电的提取，同时获得的母体心电中也含有较明显的胎儿心电成分，即盲分离的效果很差.由此可见，信号间的线性相关度测量值反映了原始混叠信号是否满足时域稀疏线性叠加盲分离模型的程度，其值越大，满足程度越高，因此越适用于常规的线性混叠盲分离处理.因此，本文提出的判别方法是有效的.

其次，通过表 2~4以及MES_n、MES_f进行比较，本文发现噪声的存在对观测信号中心电信号的相关性存在极大的破坏，使得上述的判别方法可靠性变差，甚至出现判别结果与真实情况相反的情况.因此，在进行相关性分析前，对信号进行相应的去噪处理是非常必要的.

5 结语

基于胎儿心电监护的优势，越来越多的科研人员关注这一领域，而基于盲分离技术的胎儿心电的提取更是关注的其中一个热点.为此，如果能快速、准确地判定观测信号并使之满足盲分离的线性混叠模型，将能为盲分离处理算法应用于胎心电的研究提供有效帮助，具有重要的现实意义.

本文提出的基于心电信号时域稀疏段的相关性分析方法步骤简单且运算量小，可以快速判定观测信号中的胎儿心电与母体心电是否满足线性混叠的盲分离模型.同时，从实验仿真结果看出，相关性测量值较好地反映了信号间的线性混叠程度.在同样通道数量的观测信号中，当测量值较大其高于某一阈值时，观测信号中的心电成分满足线性混叠模型，仅利用常规的盲分离手段，比如FastICA，就可以得到良好的胎心电提取效果; 反之，当测量值较小时，常规的盲分离处理则未能达到提取胎心电的目标.

再者，考虑到实际临床采集环境的复杂性，采集到的观测信号中往往包含了各种非心电的噪声成分，这些噪声都将不同程度地破坏信号中心电信号的相关度.因此，为了使测量值能准确反映心电信号间的相关度，相应的去噪处理将是非常必要的.

最后，通过仿真实验表明，本文提出的噪声环境下基于胎心电检测信号相关性判别的盲分离分析方法满足快速、有效、准确地判断观测信号是否满足线性混叠盲分离模型的要求，可为线性混叠的盲源分离问题提供一种简单有效的验证及判据，或转而寻找更优的胎心电信号处理与盲分离相结合的算法.