《1 引言》

1 引言

人体信号可分为两类:一类是主动信号,它是人体自发地发出的电、磁、光、热、声等信号;古一类是被动信号,它是外源作用于人体后所产生的信号,外源也是载体,运载着人体的信息。临床上常用的外源有X射线、激光、超声、红外等等。

人体既可以作为一种有源网络产生主动信号,还可以作为一种无源网络在外源的激励下产生响应,即产生被动信号。

《图1》

图1 人体信号的分类

Fig.1 The classification of the signals from human body

人体信号具有强随机性和高背景噪声的特点:

1) 随机性    影响人体信号的因素很多,很难用确定的数学函数来描述,现已能获得其统计特性。另外,人体信号作为随机信号还具有非平稳特性,即借号的统计特性随时空而变。

2) 背景噪声强     人体信号往往处于强的背景噪声之中,强噪声弱信号的人体信号给检测带来图难。噪声n(t)与信号s(t)绪合的方式有三种形式:加法性噪声,x(t)=s(t)+n(t);乘法性噪声,x(t)=s(t)· n(t);卷积性噪声,x(t)=s(t)*n(t)。对于这三类噪声,有不同的借号处理方法可将信号和噪声分开。

在实际人体信号研究过程中,还有正问题和逆问题之分。以心电为例,从人体的体表检测心电为正问题,而从检测到的心电图来推断心脏的电活动则是逆问题。广义来说,医生应用心电图进行诊断,其实质就是在求解心电逆问题。

本文着重闸述超声血流的无创伤检测和医学信号的特征提取,其中既有人体血流信息检测的正问题,又包括了经信号处理和特征提取后为临床诊断提供参考或依据的逆问题。

《2 超声血流信息检测的准确性研究》

2 超声血流信息检测的准确性研究

《2.1 血流Doppler效应的机理》

2.1 血流Doppler效应的机理[1]

1960年日本的里村茂夫首次将Doppler效应应用于医学中,开创了无创伤检测人体内血流速度的历史。但在以往的模型中,假设介质是静止的,而忽略了介质(血浆)也是运动的事实。因此,以往的理论很隼解释运动的没有红细胞的血浆中也存在Doppler效应现象,以及在流途较大时,常出现的计算与实验间的偏差。

针对上述问题,对于血液中的Doppler效应,提出了一种新的解释。

在运动介质中传播声波,其传播速度是介质运动速度和声波在此介质静止时的传播速度的矢量合成。血浆作为介质可认为是一种连续流体(悬浮液),红细胞则是悬浮体,它相对于悬浮液有速度差△ν。△ν的大小与血液的粘度有关,△ν取负号是因为红细胞的旋转和弯曲,使得它的速度总慢于血浆速度。

计人红细胞速度与血浆途度的差异后,可以得出用超声Doppler效应测量血流速度的新公式。

此外,还可以解释测量高速时实验值偏离理论计算的原因。在血流速度大于80cm/s时,计算值与实验值之间的偏差在于红细胞潜后于血浆的速度。用Doppler效应的检测仪器所测得的速度是红细胞的速度,而并非血挤途度。在计人血流量时要考虑到红细胞湍后校正。为此,引入了校正电路。

《2.2 起声定量测量血流速度的方法》

2.2 起声定量测量血流速度的方法[2~8]

用超声Doppler效应测量血流速度,具有无创伤的优点,但也存在不能定量测量的缺点。其原因:待测速度o和(声束与血流间)夹角0是描述Doppler效应的方程式(1)中存在两个未知数:

\(f_{\mathrm{D} }=f\left(\theta, v\right)\)   （1）

未知数多于方程数,不存在唯一解。这里通过人工设计增加方程式,使得方程数和未知数一致。求得待定未知数的唯一解,使得定量测速成为可能。

首先,采用两路独立的超声束的Doppler效应,它们的接收借号U1,U2与θ1,θ2间存在关系

\(\begin{array}{l} U_{1}=K_{1} f_{\mathrm{D} 1}=f_{1}\left(\theta_{1}, v\right) \\ U_{2}=K_{2} f_{\mathrm{D} 2}=f_{2}\left(\theta_{2}, v\right) \end{array}\)    （2）

式(2)中K1,K2为频率一电压转换常数。其次,设计一种探头,使两个超声东信号存在约束关系。

\(g\left(\theta_{1}, \theta_{2}\right) = 0 \)    (3)

籍此,三个方程中有三个未知数θ1、θ2、ν,就可以得到ν的唯一解。

在技术上,提出了三种不同的实现方法:

1) 求得两道输出U,U的差可以作为最佳位置指示值

\( \left|U_{1}-U_{2}\right|=K v \sin \theta_{0} \sin (\Delta \theta) \)    (4)

式中:0\(\theta_{0}=\frac{\theta_{1}+\theta_{2}}{2} \);\(\Delta \theta=\frac{\theta_{1}-\theta_{2}}{2}\);   K为与发射频率f0、声速c有关的一个常数。

将两道输出U1,U2的线性组合组成补偿

\(\begin{array}{l} U_{0}=\left(U_{1}+U_{2}\right) / \cos \theta_{0}+n\left|U_{1}-U_{2}\right| / \sin \theta_{0}= \\ K v[\cos (\Delta \theta)+n \sin (\Delta \theta)] \end{array}\)    (5)

通过调节常数,可以在较大的偏角范围内,输出U随角度的误差小于规定值。这就是在一定精度下的定量测速。

2) 在最佳位置的附近由三个方程解出

\(\begin{array}{c} \sec \theta_{1}=\sec \left\{\operatorname { a r c t a n } \left[\left(U_{1} \cos 2 \theta_{0}-U_{2}\right) /\right.\right. \\ \left.\left.\left(U_{1} \sin 2 \theta_{2}\right)\right]\right\} \end{array}\)    (6)

将此关系式事先存人只读存储器(EPROM)中,当测量时,用两个声束的输出U1,U2作为地址码去查表得secθ1,然后将此secθ1,乘以U1,即可得到与θ1无关的血流速度值。

第三种实现方法较之前两种有三点改进之处:

a. 探头处于任意位置都能定量测定血流途度,不需要像第一、二种方法那样寻找探头最佳位置。直接通过求解方程组计算流速值ν:

\(\begin{array}{c} v=\frac{c}{4}\left[\left(\frac{f_{\mathrm{D} 1}}{f_{10}}+\frac{f_{\mathrm{D} 2}}{f_{20}}\right)^{2} / \cos ^{2} \theta_{0}+\right. \\ \left.\left(\frac{f_{\mathrm{D} 1}}{f_{10}}-\frac{f_{\mathrm{D} 2}}{f_{20}}\right)^{2} / \sin ^{2} \theta_{0}\right]^{-1} \end{array}\)    (7)

式中,f10,f20为两超声源频率。

b. 可以计算出体内血管的走向△θ

\(\Delta \theta=\arctan \left[\cot \theta_{0}\left(\frac{f_{\mathrm{D} 2}}{f_{20}}-\frac{f_{\mathrm{D} 1}}{f_{10}}\right) /\left(\frac{f_{\mathrm{D} 1}}{f_{10}}+\frac{f_{\mathrm{D} 2}}{f_{20}}\right)\right]\)    (8)

c. 按实际情况进行误差再校正,可进一步提高测量流速的精度。误差再校正式为:

\(v_{\mathrm{o}}(\Delta \theta)=h(\Delta \theta) \dot{v}_{\mathrm{I}}(\Delta \theta)\)    (9)

式中ν1(△θ)为校正输入的流速;ν0(△θ)为校正输出流速;校正函数

\(h(\Delta \theta)=1 /\left(a_{1} \Delta \theta^{2}+a_{2} \Delta \theta+a_{3}\right)=1 / v_{0}(\Delta \theta)^{\prime} \)    (10)

通过流速实际值与△θ二次曲线拟合可求得a1,a2,a3系数。

《3 医学信号特征提取的新方法》

3 医学信号特征提取的新方法

这里将现代理论中的一些新概念引人医学超声学,形成医学信号特征提取的新方法。

《3.1 基于分形(Fractal)的Doppler血流信号分析》

3.1 基于分形(Fractal)的Doppler血流信号分析[9~17]

Mandelbrot于20世纪70年代未首先提出了分形概念。80年代末,分形开始应用于地学、宇宙、生物等方面。我们提出超声音频Doppler血流信息是具有自相似和没有囹定特征长度的图形、构造和现象,并将分形概念引人到超声血流信息处理中。根据这个新观点,我们在国内外首次报导了短时(10ms)的和整个心动周期的音频Doppler血流信号的分形维数分析、复杂度分析和近似熵分析;寻找分形新参数和它的提取技术,并对临床应用的前景进行了分析。这里仅给出分形维数的计算方法及其临床应用结果。关于复杂度、近似熵的讨论参阅文献[9一16]。

3.1.1 折线和网格分形维数计算方法

1) 折线分形维数计算:将待分析的时间序列

作式(11)计算:

\(\lg [L(r)]=\lg \left(\Sigma \sqrt{\left.\left(x_{i}-x_{i-1}\right)^{2}+r_{i}^{2}\right)}\right.\)   (11)

式中:xi为第i点归一化信号幅度,改变网格尺度r,得到在不同观察尺度下信号波形的不同长度Ig[L(r)]。设最小二乘法拟合直线的斜率为a,可得到该时间序列的折线分形维数估计值D为:

\(D=1+\alpha\)    (12)

2)网格分形维数的计算:对经过归一化处理的时间序列取间隔r作为网格尺度。则此时间床列的函数图像可分成m.个网格。然后统计覆盐住函数图像的网格格子数N(r)。改变网格尺度.作\(\lg N(r) \sim \lg r\)的双对数图。在图中确定出一段直线区,通过最小二乘法拟合该段直线,其斜率的绝对值就是时间序列的分形维数估计值。

在短时分形维数估计中,由于受时间序列数据长度的制约,只取采样间隔△以及2△为测量的网格尺度,可测得:

\(N(\Delta)=\sum_{j=1}^{n}\left|x_{j}-x_{j+1}\right| / \Delta \),    (13)

\(N(2 \Delta)=\sum_{j=1}^{n / 2}\left(\max \left\{x_{2 j-1}, x_{2 j}, x_{2 j+1}\right\}-\right. \\ \left.\min \left\{x_{2 j-1}, x_{2 j}, x_{2 j+1}\right\}\right) / 2 \Delta\)    (14)

该段时间序列的短时分形维数D为:

\(D=\frac{\lg N(\Delta)-\lg N(2 \Delta)}{\lg 2}\)    (15)

3.1.2  对102例胎儿的脐动脉超声血流信号进行非线性分析    102例胎儿脐动脉的超声多普勒音频血流信号取自上海医科大学附属妇产科医院,孕周数为20~42周,年龄为20~46岁,其中临床诊断为正常的57例,临床诊断存在疑问或异常的45例,血流信号经12bit A/D板椿样送人计算机。对每个病人采集10s左右的数据(大约有12~13个心动周期),然后拍心动周期分段进行整个心动周期的分形维数计算。

通过对57例正常病例的分数维值计算和统计分析,初步结果为: a. 正常的分数维值为1.847±0.018; b.正常分数维值与妊娜周数存在相关性,正常分维值随着妊娠周数增大而变大; c. 正常分数维值与孕妇年龄之间无显著相关性。

在25例病例中,分数维值在正常范围内的共有10例,大于正常范围的有6例,小于正常范的有9例,总计异常检出率为60%。它与常规声谱图参数S/D进行了比较:对阳性检出率、正确率等指标。就现有的样本而言,分形分析方法要优于常规的声谱参数法。

血流信号的非线性的分形特征同时反映了信史的频率信息和相位信息,可作为传统的血流信号频谱分析方法的补充和扩展。

《3.2 在医学超声中引入数学形态学(Morphology) 方法》

3.2 在医学超声中引入数学形态学(Morphology) 方法[17，18]

我们率先在医学超声研究中引人了数学形态学方法,主要是在传统的超声多普勒声谱参数提取技术中,提出了基于数学形态学边缘提取和形态滤波的声谱包络及包络特征点提取的新方法,并应用于超声多普勒声谱参数的自动计算中。

3.2.1 超声Doppler声谱色络的数学形态学提取方法    用a-裁剪均值形态边缘检测算子(a 一 trimmed morphological edge detector)来提取声谱包络是数学形态学在医学超声多普勒上的首次应用尝试。该方法阐述如下:

设输入信号为\(f_{1}(z) \subseteq Z^{n}\),中间结果为\(f(z) \subseteq Z^{n} \),输出信号为\(f_{\mathrm{o}}(z) \subseteq Z^{n}\)。

首先,将输入信号\(f_{\mathrm{I}}(z)\)进行a-裁剪平均模糊处理,得到中间结果 f(z):

\(f(z)=\sum_{l=\alpha+1}^{k-\alpha} \frac{f_{i l}(z)}{k-2 \alpha}\)   （16）

式中,为结构元素B包含的总点数,f(z)是用移位至z处的结构元素Bz.度量的输入信号中第l个最小的样本。

然后对中间结果进行形态边缘增强,得到输出借号fo(z):

\(\begin{array}{c} f_{\mathrm{o}}(z)=\min \left\{\left[f_{\mathrm{B}}(z)-(f \ominus B)(z)\right]\right. ,\\ \left.\left[(f \oplus B)(z)-f^{\mathrm{B}}(z)\right]\right\} \end{array}\)。(17)

最后,对输出信号进行阈值化处理,得到声谱包络。

传统的模板匹配和增强/阈值化边缘检测方法设计简单,但鲁棒性差,即随着信噪比的下降,检测性能也显著下降;而数学形态学的图像边缘检测,具有很好的统计鲁棒性,尤其是具有较高的脉冲性噪声抑制能力,在信噪比较差的情况下,仍能获得较好的图像边绍检测效果。

3.2.2 超声Doppler声谱包络特征点的数学形态学提取方法    声谱图包络特征点主要是指最大频率曲线上收缩期最大频率(速度)点S和舒张末期最大频率(速度)点D,在声谱图包络上一般表现为极大值点和极小值点。对于不同的血流状况,在一个心动周期内,S点应是最大值点,但是D点不一定是最小值点。

这里提出一种基于数学形态学的定征体系,其结构框图见图2。其中:形态滤波是对声谱包络进行光滑,降低噪声对包络提取的影响,它可以清除信号中脉冲性噪声引起的伪峰谷,更有效地提取形态峰谷。形态滤波法的结构框图见图3。它包指两个并行处理的分支,两分支的平均值作为最终的输出。滁波效果取决于结构子的形状。

《图2》

图2 数学形态学定征体系的结构框图

Fig.2 The outline of the morphology

《图3 》

图3 数学形态滤波法的结构框图

Fig.3 The outline of the method of morphology filter

峰谷提取运算的基本思想是通过使用一种算子产生一个新的波形。在新生成的波形中,收缩期的波峰被突出出来,而其它的小波峰(局部极大值点)则被映射成平坦信号区域;同样,舒张末期的波谷被映射成负波峰,而其它的小波谷(局部极小值点)也被映射成平坦信号区域。峰谷提取质量主要依赖于结构子的选取,这里选零振幅水平线作为结构子。

声谱包络经过形态转换处理后,有用的峰谷已被突出。寻找极大幅值点时,若只有一点,则定为特征点S,若不止一点,取原始波形中幅值最大点作为特征点S;对于特征点D,则由特征点S向前查找,避免舒张末期以前可能出现的反向血流造成S点误判。这样更方便地取出特征点。

3.2.3 基于数学形态学的声谱参数自动计算系统

综合数学形态学的声谱图包络自动提取法和声谱图包络特征点自动定位的方法,提出超声多普勒声谱参数自动计算方法,构成整套的基于数学形态学的声谱自动计算系统,见图4。

《图4》

图4 基于数学形态学的声谱参数自动计算系统框图

Fig.4  The outline of the spectrogram indices automatic calculation algorithm based on morphology

在小儿肺动脉、成人颈动脉以及胎儿脐动脉多普勒血流信号上,用数学形态学方法提取声谱包络和自动提取声谱包络特征点,取得了良好的效果。

《3.3 数量化理论在孕周和胎重无损估计中的应用》

3.3 数量化理论在孕周和胎重无损估计中的应用[19~24]

1983年我们首次将数量化理论(theory of quantification)应用于主动脉关闭不全的无损诊断上,近年来又研究用于无损估计孕周和胎重。对孕妇和胎儿的监护(孕周和胎重的估计是其中重要的两个参数),可以为医护人员诊断、治疗提供重要的参考依据。对孕周和胎重的传统估计是通过B型超声测量胎儿器官的有关参数,辅以经验公式而得到。估计孕周的传统参数一般为双顶径、腹径和股骨长,估计胎重的传统参数一般为肝径和股骨长,但它们在不同的阶段对孕周和胎重变化的敏感程度不同,因而影响了估计的精度;经验公式无统一形式,昆参数皆单独使用,估计结果离散度大。因此,我们研究用数量化理论来无损估计孕周和胎重。

具体估计方法为:参照数量化理论中将定性变量(称作项目)的不同取值划分为类目的方法,将估计用的参数作为项目,等间隔地划分成不同的类目,对定量的基准变量(即为孕周或胎重)进行估计。用短阵形式表示:

\(\boldsymbol{Y}=\boldsymbol{X} \cdot \boldsymbol{b}+\boldsymbol{E}\),(18)

式中,X即为反应矩阵, Y为样本矢量, b为系数矩阵,E为误差矩阵。

因变量估计值\(\hat{\boldsymbol{Y}}\)的表示式为:

\(\hat{\boldsymbol{Y}}=\boldsymbol{X} \cdot \hat{\boldsymbol{b}}\)。(19)利用广义逆短阵法求解

\(\hat{\boldsymbol{b}}=\left[\boldsymbol{X}^{\mathrm{T}} \cdot \boldsymbol{X}\right]^{+} \cdot \boldsymbol{X}^{\mathrm{T}} \cdot \boldsymbol{Y}\)。(20)

这里A+为短阵A的广义逆短阵,求出短阵的奇异值A分解,便可得到广义逆短阵A+。可以分两阶段计算A的奇异值分解,即第一阶段利用Householder变换把短阵A化为双对角短阵形式;第二阶段利用Givens反射变换的QR算法求双对角短阵的奇异值。

为了建立孕周的估计公式,采集了390例样本作为胎儿孕周佼计的原始数据。采用的参数除了传统的双顶径、腹径和股骨长之外,还选用了可能对孕周变化敏感的胎儿眼球纵径、眼球横径和耳廓长这3个参数。把这6个参数作为6个项目进行划分。TQ中的项目及类目划分见表1。

《表1 》

表1 项目中类目的划分

Table1 ltem division

项目

最小值／mm

最大值／mm

间隔／mm

类目数

眼球纵径

4.0

16.4

0.5

27

眼球横径

4.8

16.8

0.5

26

耳廓长

9.5

34.0

1.0

27

双顶径

34.0

96.0

1.0

64

腹径

30.0

98.0

1.0

70

股骨长

16.0

76.0

1.0

62

总计

 

 

 

276

用另外60多例临床数据对孕周估计公式(19)进行了检验,并将该方法与以往方法及线性回归法等其他几种方法作了比较,见表2。其中数量化I采用删除法求解短阵8,数量化II采用广义逆矩阵法求解矩阵。可见,改进后的采用6参数的数量化方法效果最好,优于以往的方法(包括采用6参数线性回归方法),是临床上应用医学超声无损估计孕周的实用方法。

表2 各种孕周估计方法的比较

Table 2    Comparisons of the various gestational age estimation methods

方法

胎儿发育

阶段／周

参数个数

／个

例数

／个

绝对误差

平均值／周

标准差

／周

数量化方法II

14-40

6参数

63

0.6479

1.033

数量化方法1

14-40

6参数

63

0.9862

1.347

线性回归

14-40

6参数

63

1.154

1.510

线性回归

14-40

双顶径

63

1.582

2.037

线性回归

14-40

股骨长

63

1.474

1.830

线性回归

14-40

腹径

63

1.402

1.791

经验公式

14-40

双顶径

63

2.969

2.704

经验公式

14-40

股骨长

63

2.222

2.486

同样方法可对胎重进行估计。采用170例作为训练集,18例作为测试,所用的参数为肝径(类目数取为55个)和股骨长(类目数取为21个)。结果如表3所示。

《表3》

表3 胎重数量化理论的估计

Table 3 Gestational weight estimation based on the quantification theory

 

训

练

测

试

平

均误差/g

标准误差／g

平均误差/g

标准误差/g

数量化理论法

141

191

346

615

线性回归法

200

278

209

300

《3.4 小波变换软阈值降噪方法的研究和应用》

3.4 小波变换软阈值降噪方法的研究和应用[22~34]

许多具有临床价值的血流参数都是从声谱图中获得的,如S/D(收缩期与舒张未期最大流途比)、PI(脉动指数)等。为了避免因声谱图畸变给参数提取造成的困难,对声谱图进行增强。由于音频多普勒血流信号的功率谱是时变的,传统的仅以优化均方误差为目标的降噪算法并不适合。除了优化均方误差外,对音频多普勒血流信号在降噪中尽量减少由降噪本身所带来的不必要的频率成份。为此,提出基于小波变换的软闻值降噪算法,对音频多普勒信号进行降噪从而增强声谱图。

从含噪声的信号Q,中恢复未知函数广的降噪步骤为:

1) 利用预先定义的塔式算法对归一化的测量信号\(\left(d_{i} / \sqrt{n}\right)\)作多分辨率分析,得到小波变换系数\(\omega_{j, k}\)其中,\( n=2^{J+1}, j=j_{0}, \cdots, J, k= 0, \cdots, 2^{j-1}\), j0和J分别为预先定义的分析层数;

2) 对所得小波变换系数进行收缉,引入软阈值法,即加阈值\(\beta=\sqrt{2 \lg (n)} \sigma / \sqrt{n}\),得到估计的小波变换系数\(\hat{\alpha}_{j, k}\),即

\(\hat{\alpha}_{j, k}=\left\{\begin{array}{ll} \omega_{j, k}-\beta & \omega_{j, k} \geqslant \beta \\ 0 & -\beta<\omega_{j, k}<\beta \\ \omega_{j, k}+\beta & \omega_{j, k} \leqslant-\beta \end{array}\right.\)

与硬闻值法不同,软闻值方法并不是简单地将绝刑值小于阅值的系数设为零,而是将大于闻值的系数收缩至零点附近,避免了硬阅值方法中产生不连缎点而引起不必要的频率成份;

3) 对所有j>J,设\(\hat{\alpha}_{j, k} = 0\),利用;\(\hat{\alpha}_{j, k}\)进行小波重构得到一的估计\(\hat{f}\)。

利用模拟多普勒信号和临床多普勒信号进行降喋实验。实验表明,小波结合软闻值的降噪新方法比以往的FFT闻值法以及硬阎值法有更好的谱线增强效果(见图5、图6)。该技术已被应用到肺循环血液动力学参数无损检测系统中,临床试用效果良好。

《3.5 极点轨迹》

3.5 极点轨迹[26~28]

英国Bristol总医院的Skidmore和Woodcock,首次提出了利用三防极点模型对一个心动周期血洗平均频率(速度)曲线进行建模[29]。但血管内血流速度具有一定的分布,而且流速分布在一个心动周期内是不断变化的,反映为Doppler信号的非平稳性。因此,在三阶极点的基础上,我们直接对含有比平均频率曲线更丰富信息的音频Doppler信号进行建模。通过对短时间间隔内的音频Doppler信号进行极点分析,研究整个心动周期不同时相的动态极点轨迹。

《图5》

图5 模拟多普勒信号声谱图增强结果
Fig. 5    Results of the spectrogram enhancement of the simulated Doppler signal 

《图6》

图6 临床多普勒信号声谱图增强结果
Fig. 6    Results of the spectrogram enhancement of the clinical Doppler signal 

音频Doppler信号在短时间T内可以看作是准广义平稳的随机信号。我们的研究表明,随心动时相而变,收缩期T约为8ms,舒张期约为12ms[30]。在实际应用中,通常取工=10ms[31]。在时间段工内的音频Doppler用一个三阶极点模型表示:

\(H(Z)=\frac{b_{1} Z^{-2}+b_{2} Z^{-1}+b_{3}}{Z^{-3}+a_{1} Z^{-2}+a_{2} Z^{-1}+a_{3}}\)(21)

稳定网络的三个极点在Z平面上的分布一舫为:一个实轻极点,两个共轭极点(有时会退化为两个实极点)。实极点与该系统的能量消耗特性有关,复极点与该系统的振荡、阻尼特性有关。

设信号的采样频率为人,则在时间段T中采集得N=fsT点数据\(x(0), x(1), \cdots, x(N-1) \),于是H(z)又可以写成

\(\begin{aligned} H(z) &=x(0)+x(1) z^{-1}+\cdots+\\ & x(N-1) z^{-(N-1)} \end{aligned}\),(22)式中,x(n)为信号序列。由式(21)与式(22)

可得

\(\begin{array}{l} \left\{\begin{array}{c} -x(0)=a_{3} x(3)+a_{2} x(2)+a_{1} x(1) \\ -x(1)=a_{3} x(4)+a_{2} x(3)+a_{1} x(2) \\ \vdots \\ -x(N-4)=a_{3} x(N-1)+ \\ \quad a_{2} x(N-2)+a_{1} x(N-3) \end{array}\right. \\ \end{array}\)                     (23)

\(\begin{array}{l} \left\{\begin{array}{c} b_{3}=a_{3} x(0), \\ b_{2}=a_{3} x(1)+a_{2} x(0) \\ b_{1}=a_{3} x(2)+a_{2} x(1)+a_{1} x(0) \end{array}\right. \end{array}\)                           （24）

式(23)是一个超定方程组,可以用伪逆法由式(23)解出al,a2,as,再由式(24)求得系数,多,53。根据al,aa和a又可得出该模型的三个极点,其中两个为共轮的复极点。由于多普勤信号是随机过程,再考虑到模型的因果性和稳定性,所有的极点均应在单位圆内。

利用上述方法分析一个心动周排或更长时间内的各段音频多普勒信号,就可以得到在这段时间内的极点分布情况,从中可以提取出与疾病症状、生理特征有关的参数。该模型曾用于颈动脉血流的分析(上海第一人民医院)21和脐动脉血流的分析(上海医科大学妇产科医院)7291。

脐动脉是胎儿与母体显养、物质交换的唯一通道,一端源自胎儿心胃,一端终止于胎盘。其终端负载一胎盘的营养物质交换状况及胎盘本身发育状况,直接反映、影响胎儿的生长发育状况,胎盘的功能状态更是临床上关心的重要信息。

分析了47例的孕妇脐动脉血流,其中正常孕妇30例,临床诊断为IUGR或其它妊娜期疾病的孕妇17例。图7a给出了一例正常升妇脐血流信号极点随时间变化的分布情况,而图7b给出了一例临床诊断为IUGR的孕妇,其脐血流信号极点随时间变化的分布情况。

《图7》

图7 正常孕妇、IUGR孕妇脐动脉血流的极点分布

Fig.7 The poles distribution of the umbilical artery blood from normal and IUGR pregnant woman

由图7可见,正常升妇的极点分布较为集中且靠近单位圆,而IUGR孕妇的极点分布较为分散且远离单位圆。极点的分布可以在一定程度上反映胎儿的发育情况。

用上半平面内所有复极点和实极点的重心来定量描述极点分布, 分别记为, \( P_{\mathrm{G}}=\sum_{i=1}^{N} P_{i} / N_{\mathrm{P}}, R_{\mathrm{G}}   =\sum_{i=1}^{N_{\mathrm{R}}} R_{i} / N_{\mathrm{R}}\) , 式中, \( P_{i}\left(i=1,2, \cdots, N_{\mathrm{P}}\right)  \)是上 半平面内的复极点,\(  R_{i}\left(i=1,2, \cdots, N_{\mathrm{R}}\right)  \)是实 极点。接着进一步定义有关重心的 5 个参数: \( P_{\mathrm{G}}  \)的幅长\(  r_{\mathrm{P}}\) 、  幅角\(  \phi_{\mathrm{P}} \)、  上半平面内所有复极点相对\(  P_{\mathrm{G}}  \)的离散程度\(  v_{\mathrm{P}} 、 R_{\mathrm{G}}  \)的幅长\(  r_{\mathrm{R}}  \)以及所有实极点 相对  \(R_{\mathrm{G}}  \)的离散程度\(  v_{\mathrm{R}}  \)。

\(\begin{array}{c} r_{\mathrm{P}}=\left\|P_{\mathrm{G}}\right\|, \phi_{\mathrm{P}}=\arctan \frac{\operatorname{imag}\left(P_{\mathrm{G}}\right)}{\operatorname{real}\left(P_{\mathrm{G}}\right)} \\ r_{\mathrm{R}}=R_{\mathrm{G}}, \\ v_{\mathrm{P}}=\sqrt{\frac{\sum_{i=1}^{N_{\mathrm{p}}}\left\|P_{i}-P_{\mathrm{G}}\right\|^{2}}{N_{\mathrm{P}}}} \end{array}\)
 

\(v_{\mathrm{R}}=\sqrt{\frac{\sum_{i=1}^{N_{\mathrm{R}}}\left\|R_{i}-R_{\mathrm{G}}\right\|^{2}}{N_{\mathrm{R}}}}\)               (25)

图8所示是这5个特征参数相对于孕周数(Weeks)的分布情况。由图8可见,IUGR孕妇的rp与正常人相比偏小,但差异并不明显。对于参数加、op,正常孕妇与IUGR孕妇之间没有明显的差异。而对于表征实极点分布情况的参数r和og,IUCR孕妇的参数r较正常人为小,而ox参数较正常人为大,其中以og参数的敏感性尤为明显,30例正常孕妇中只有1例明显偏大,而17例IUGR或其它妊娜期疾病孕妇中只有4例明显偏小。

脐动脉的终端一胎盘可以看成是它的外周血管床。正常孕妇的胎盘毛细血管发育完全、血流量大、营养物质交挨充分、胎儿发育良好,表现为胎盘阻力较小。而IUGR孕妇的胎盘毛细血管发育不好,血液灌注量小,表现为胎盘阻力较大,胎儿发育不良。这样,实极点的分布与胎盘发育状况密切相关,使实极点指数RI的ra和ox可能成为诊断IUGR的判据,如图9所示。可见,仅有两例正“谱参数RI相对孕周的分布情况。由图10可见,正常孕妇的特征参数画常,仅有三例妊娜期疾病孕妇常孕妇的R参数与异常孕妇相比较小,但两者的的特征参数落人正常区,与其它手段得到的判断相“差异不显著。对声谱参数S/D和PI的研究也发现比,有很好的相关性。图10给出了这些孕妇的声谱谱参数RI相对孕周的分布情况。由图10可见，正常孕妇的RI参数与异常孕妇相比较小，但两者的差异不显著。对声谱参数S／D和PI的研究也发现类似的情况。

《图8》

图8 各特征参数相对于孕周数的分布情况

Fig.8 The relationship between various indices and gestational ages

《图9》

图9 rR、VR参数的两维分布图

Fig.9 Two dimension distribution of rR、UR indices

《图10 》

图10 声谱参数RI相对孕周数的分布情况

Fig.10 The relationship between the spectrogram index RI and the gestational ages

通过建立脐血流的极点模型,分析极点的分布情况,可提取一定的特征参数作为妊娓期疾病的判据。用极点模型提取出的两个特征参数rk、oR可能比传统超声多普勒声谱参数更能反映病理信息,可为医学临床的诊断提供一种新的有效指标。这种技术的临床初步试用表明,它具有有效性、可重复性及无损伤性等优点1,可将它作为新的超声仪器设计恺路。系统结构可以是将检测和建模集成为一体化的仪器,也可以将建模这部分的软硬件设计成后置式,医院现有的超声仪器(有音频多普勒信号输出)为前级,新增的后级是带采集装置的PC机,在PC机上实现极点的计算和特征提取。

如果把胎盘看作一个系统,把脐动脉血流信号作为输入,脐静脉血流信号作为输出,则是个非线性系统。而极点模型是个线性模型,可能会丢失棠些具有生理和病理意义的非线性信息。如果进一步对脐血流信号建立非线性模型,可能会取得更好的效果。

《3.6 脐血流阻抗的血管传输线模型估计》

3.6 脐血流阻抗的血管传输线模型估计［31,32］

从宏观传播方式上看，血流和电流具有相似性。将血流（波）在血管中的传输类比为电磁波在传输线中的传播，这就是建立血管传输线模型的基本思路。

脐动脉起自胎儿心脏，止于胎盘。脐动脉血管粗细均匀，且无分支。可以将整段脐动脉看作由一系列无限小单位长度的小段依次串联组成。脐动脉前端看作为传输线的信号源，终端一胎盘看作传输线的负载，如图11所示。初步临床应用如图12所示。随着妊娠周数增加，胎儿所需血液的灌注量也增大，血管管径增大，血流量增加，表现为阻力减小，图12中显示出脐血流阻抗下降。在22～30周下降趋势更为明显。而至32～33周，血管床发育日渐成熟，胎盘灌注量变化不大，脐血流阻抗变化平缓（图12）。正常胎儿与异常胎儿的阻力有显著差异，相差30％～50％。在小孕周，差异更为明显。临床数十例的数据初步表明，用传输线模型得到的脐血流阻抗是胎儿生长发育状况的一种估计，有一定应用价值。

《4 应用系统的研制》

4 应用系统的研制

依据上述的新概念、新方法研制了几套系统，下面介绍其中主要的三项。

《图11》

图11 血管传输线模型示意图

Fig.11 The illustration of the transmission line model of the vessel

《图12 》

图12 不同孕周正常胎儿血流阻抗分布图

Fig.12 The distribution of the umbilical blood impedance of the fetal with various gestational ages

《4.1 肺动脉血液动力学参数的无创伤估测系统》

4.1 肺动脉血液动力学参数的无创伤估测系统［33］

肺动脉血液动力学参数是先心病、肺心病等疾病的临床诊断、手术指征、预后判断的重要依据。血液动力学参数中以肺动脉压的估测最为重要。估测肺动脉压，并借助超声心动图可估测全肺阻力、肺血流量［35］。肺动脉压常采用有创伤的心导管术测量。我们提出一种基于超声和心电的无创伤估测方法，研制了肺动脉压无创伤估测系统。

肺动脉压PPA无创伤估测公式为：

\(\begin{array}{l} p_{\mathrm{PA}}=\left(F_{\mathrm{PA}} / F_{\mathrm{AO}}\right) \cdot p_{\mathrm{BA}} \\ F_{\mathrm{PA}}=t_{\mathrm{PEP}} /\left(t_{\mathrm{PAP}} \cdot t_{\mathrm{ET}}\right) \\ F_{\mathrm{AO}}=t_{\mathrm{PEP}} /\left(t_{\mathrm{AOP}} \cdot t_{\mathrm{ET}}\right) \end{array}\)    (26)

式中：PBA为体外测得的血压，tpEp为收缩前期间期，tPAP、tAOP为肺动脉和主动脉血流速度的上升时间（起始点至峰值），tET为心射血时间。肺动脉压无创伤估测系统的框图见图13。

《图13 》

图13 肺动脉压无创伤估测系统框图

Fig.13 The architecture of the system of the non-invasive estimation of the pulmonary pressure

系统采用了如下的关键技术：a.音频Doppler信号小波变换软阈值降噪［24,25］；b．改进百分比法提取声谱图的包络线［24,36］；c．基于Mallat算法的声谱包络和心电信号特征点的识别和定位［24,37］。

系统正在上海第二医科大学新华医院和上海市儿童医学中心试用。经动物实验（10条狗）和临床应用（117例）结果表明：效果良好，经与心导管方法对照，取得了较好的相关性。

《4.2 彩色编码声谱仪》

4.2 彩色编码声谱仪［37,38］

国内的Doppler血流（非实时）声谱分析，80年代湖北物理研究所、清华大学最早报导。复旦大学于1988年在国内首次研制成实时声谱仪；1992年首次研制成彩色编码实时声谱仪。嗣后，由江西、江苏的两个工厂进行科研成果产业化。产品在医院应用病例近万例，取得良好的效果［39~41］。
该仪器的硬件由一个 Doppler 超声检测单元和多CPU系统组成。Doppler单元完成超声信号的发射、接收和解调，然后由多CPU系统处理。3个CPU独立处理各自任务，CPU-1为控制单元，CPU-2为谱分析单元，双CPU完成实时谱分析，并以扩展卡方式进入CPU-3的主机单元（PC机），实现声谱图的彩色编码实时显示、图形冻结、参数测量、文件管理及打印输出等功能。

《4.3 超声血流速度定量检测卡和脑血管参数无损测量系统》

4.3 超声血流速度定量检测卡和脑血管参数无损测量系统

血流速度的定量检测从模拟［2,4,5］、数字［6］进入计算机化［7,8］。将超声血流速度定量检测设计成计算机的扩展卡，可组成各种不同用途的系统。脑血管参数无损测量系统就是其中之一。

脑血管参数无损测量系统包括 Doppler 测速、压力传感、接口、PC机、打印等部分。血流速度和压力信号由检测部件获得，通过接口送入PC机，结合有关的病例资料，在PC机内算出指定的参数并显示，最终结果送硬盘存贮和打印机输出。
系统的参数有四大类：反映脑血管供血状况的运动学指标（颈动脉血流速度）；反映脑血管外围阻力的指标（外周阻力和动态阻力）；反映脑血管弹性的指标（特性阻抗，脉搏波速）；反映脑血管闭锁状态的指标（临界压、舒张压与临界压差）。
超声血流速度定量检测卡还应用在数种计算机化的医学仪器中。