為研究T波交替(TWA)的定量檢測方法,探討散點圖外部形態與TWA的關系,研究出適用于散點圖的一種邊緣提取算法"橫向搜索法",然后對TWA散點圖作圖形化處理,并以圖形的短長軸之比Axial_ratio作為定量檢測指標。采用Matlab仿真并與譜分析法比較,找出合適閾值并判斷TWA是否存在,證明由"橫向搜索法"提取的Axial_ratio能夠準確檢測TWA。
引用本文: 陳甜甜, 趙捷, 尹文楓, 張春游, 李大鵬, 安佰京, 張軍建. T波交替散點圖法的定量檢測. 生物醫學工程學雜志, 2014, 31(3): 538-542. doi: 10.7507/1001-5515.20140100 復制
引言
T波交替(T-wave alternans,TWA)是指在規則的心律時,體表心電圖上T波振幅、形態甚至極性的逐搏交替變化,它與惡性室性心律失常的發生密切相關[1]。TWA一般為微伏級,肉眼可分辨的TWA現象非常少[2]。根據TWA檢測統計方法不同,將TWA檢測分為短時傅里葉變換法、符號變換法和非線性法[3]。譜分析法是短時傅里葉變換法的一種,對64~128個ST-T結構中的T波采樣值進行快速傅里葉變換,將心電圖上T波段幅度的變化轉變成能量譜,分析0.5周期/心跳處的頻率成分,來測定是否存在TWA[4]。龐加萊映射[5](Poincare mapping,PM)是非線性法的一種,利用散點圖適于描述周期性變化數據的特性,可以使TWA中T波幅值逐拍變化的特點從形態上直觀表現出來。本文首先對散點圖進行圖形處理并定義散點圖短長軸之比為量化指標,通過與譜分析法的檢測結果進行相關性分析,證明散點圖法是檢測TWA的有效手段,并提出了一種新的TWA檢測有效量化指標。
1 算法介紹
1.1 心電信號預處理
首先通過簡單整系數50 Hz陷波器、零相位數字濾波器、結合小波分解與重構理論和改進的閾值算法去除工頻及其諧波干擾、基線漂移、肌電干擾等心電噪聲;然后采用二階微分Marr小波,應用多孔算法對ECG中QRS波群進行標定[6],在S波峰后150~350 ms的范圍內搜索極大值得到T波峰;隨后,對系列心跳進行復合預處理[7],提高TWA檢測的精確度。
1.2 TWA檢測
1.2.1 T波采樣點的選取
為了更準確的檢測微伏TWA(MTWA),我們采用N個心跳周期的改進T波窗口法。根據RR間期選擇T波窗口起點,通過Bazzet 公式確定T波窗口終點,T波選取按經驗原則,如表 1所示。

以T波波峰對齊,按T波平均長度的1/6為間隔,前后各取3個點構成7點檢測集xi={xi1,xi2,…,xi7},其中i代表所識別的心跳序列序數。本文選取128個T波段,每段提取7個采樣點,構成T波采樣點序列P={x1,x2,…,xi,…}。
1.2.2 散點圖法檢測TWA
對P={x1,x2,…,xi,…}進行一次差分組成新的序列Q={x2-x1,x3-x2,xi-xi-1,…},以X={x2-x1,x3-x2,…,xi+1-xi,…}為橫坐標序列,Y={x3-x2,x4-x3,…,xi+2-xi+1,…}為縱坐標序列,作Poincare散點圖,則該散點圖中點的坐標為(xi+1-xi,xi+2-xi+1)。本文選取128個T波段,每段提取7個采樣點,可作出882個散點。
若存在TWA,相鄰T波的采樣點幅值會以ABAB形式交替變化,即T波采樣點序列P={x1,x2,…,xi,xi+1,xi+2,…}的值會交替變化,則(xi+1-xi)與(xi+2-xi+1)異號且絕對值相當,也就是散點(xi+1-xi,xi+2-xi+1)主要以(+,-)或(-,+)的規律變化,此時散點圖呈長軸在直線y=-x上的傾斜橢圓形狀。若無TWA,相鄰T波的采樣點幅值呈AAAA形式,此時散點以原點為中心均勻分布,散點圖接近圓形。
經Matlab仿真繪圖,TWA散點圖外觀上有較明顯的邊緣[8],如圖 2所示。此邊緣是散點圖的最基本特征。為了方便對散點圖進行特征描述、分析,以及大量縮減待分析數據量,我們對TWA散點圖進行邊緣提取。根據其散點分布特點,研究出“橫向搜索法”來確定散點圖的“邊緣”,其具體步驟如下: ① 散點圖預處理。將散點圖轉換成M×N的二值化的數字圖像矩陣G。根據其橫、縱坐標序列X={x2-x1,x3-x2,…,xi+1-xi,…} 、Y={x3-x2,x4-x3,…,xi+2-xi+1,…}及散點坐標(xi+1-xi,xi+2-xi+1),首先將散點圖坐標序列值取整后上移和右移一定單位,使其剛好完全移至平面直角坐標系第一象限內,然后通過如下方式進行轉換:
$\begin{align} & G=\left[ \begin{matrix} {{g}_{11}} & {{g}_{12}} & \ldots & {{g}_{1N}} \\ {{g}_{21}} & {{g}_{22}} & \cdots & {{g}_{2N}} \\ \vdots & \vdots & \ddots & \vdots \\ {{g}_{M1}} & {{g}_{M2}} & \cdots & {{g}_{MN}} \\ \end{matrix} \right], \\ & (M=max\left[ Y\left( n \right) \right],N=max\left[ X\left( n \right) \right],n=1\,\!:L) \\ \end{align}$ |
其中Y(n)、X(n)分別表示Y、X中的第n個元素值;max[Y(n)],max[X(n)]分別表示X、Y中所有元素的最大值;矩陣元素gi,j表示圖像G中第i行,第j列的像素灰度值,其定義為
${{g}_{i,j}}=\left\{ \begin{align} & 1,\left( i=Y\left( n \right),j=X\left( n \right),n=1\,\!:L \right) \\ & 0,其他 \\ \end{align} \right.$ |
從而得散點圖的二值化數字圖像G。
② 尋找邊緣候選點。對圖像G進行塊操作,初步去除散點圖中心非邊緣點及邊界外較稀疏的非邊緣點。定義5×5的整數格,令整數格以gi,j,(i=1…M,j=1…N)為中心像素開始對散點圖進行行掃描,計算該整數格覆蓋區域A(見下式)的灰度值和a,即A中各元素之和,若2≤a≤3,則將該中心像素標記為邊緣候選點。
$A = \left| {\matrix{ {{g_{i - 2,j - 2}}} & {} & \cdots & {} & {{g_{i - 2,j + 2}}} \cr {} & \ddots & {} & {} & {} \cr \vdots & {} & {{g_{i,j}}} & {} & \vdots \cr {} & {} & {} & \ddots & {} \cr {{g_{i + 2,j - 2}}} & {} & \cdots & {} & {{g_{i + 2,j + 2}}} \cr } } \right|$ |
③ 確定邊緣點。對由上個步驟確定的邊緣候選點進行篩選,去除散點圖中心區域外的一些異常小點集,得到最終邊緣點。首先以邊緣候選點為中心像素構建21×21的整數格,則該整數格覆蓋區域B,即
$B = \left| {\matrix{ {{g_{i - 10,j - 10}}} & {} & \cdots & {} & {{g_{i - 10,j + 10}}} \cr {} & \ddots & {} & {} & {} \cr \vdots & {} & {{g_{i,j}}} & {} & \vdots \cr {} & {} & {} & \ddots & {} \cr {{g_{i + 10,j - 10}}} & {} & \cdots & {} & {{g_{i + 10,j + 10}}} \cr } } \right|,$ |
將矩陣B按B=劃分,得矩陣B1、B2、B3、B4,分別計算其灰度值和b,若所有b值都小閾值6,則刪除此候選點,將其余點保存在邊緣點數組中,得散點圖的邊緣點集。
為了更好地度量TWA的幅值大小,我們需要一種定量的檢測標準[9]。本文提出了短軸、長軸和軸比三個定量指標,并對軸比進行仿真測評。 我們設邊緣點集為Ω,(X,Y)表示Ω中點的坐標。以直線y=-x+z對Ω進行遍歷,當縱截距z(z=x+y)的值取最大和最小時,即分別找到該邊緣點集的兩條切線BC:y=-x+zmax和AD:y=-x+zmin。同理,以直線y=x+b對Ω進行遍歷,分別找到該邊緣點集的兩條切線AB:y=x+bmax和CD:y=x+bmin。由此得到該邊緣點集的最小外接矩形ABCD,散點圖幾何特征如圖 1所示。

① 定義短軸(short_axis)為最小外接矩形ABCD的短邊長,即:
$\begin{align} & Short\_axis=\left| AB \right|=\left| CD \right|=\frac{{{z}_{max}}-{{z}_{min}}}{\sqrt{2}} \\ & (z=X+Y,\left( X,Y \right)\in \Omega ) \\ \end{align}$ |
② 定義長軸(Long_axis)為最小外接矩形ABCD的長邊長,即:
$\begin{align} & Long\_axis=\left| BC \right|=\left| AD \right|=\frac{{{b}_{max}}-{{b}_{min}}}{\sqrt{2}} \\ & (b=X-Y,\left( X,Y \right)\in \Omega ) \\ \end{align}$ |
③ 定義軸比(Axial_ratio)為短軸與長軸的比值,這個指標能夠描述該邊緣點集沿直線y=-x的伸長程度。指標值越小,散點分布越趨近于直線y=-x;反之,散點分布越接近圓形,即:
$\begin{align} & Axial\_ratio=Short\_axisLong\_axis=\frac{{{z}_{max}}-{{z}_{min}}}{{{b}_{max}}-{{b}_{min}}} \\ & \text{ }(z=X+Y,b=Y-X,\left( X,Y \right)\in \Omega ) \\ \end{align}$ |
經Matlab仿真分析,規定Axial_ratio的閾值為0.57,即Axial_ratio≤0.57時,存在TWA;Axial_ratio>0.57時,無TWA。圖 2為TWA檢測的散點圖。

(a)不含TWA的散點圖;(b)含TWA的散點圖
Figure2. Poincare maps of TWA detection(a) the Poincare mapping without TWA; (b) the Poincare mapping with TWA
2 結果分析
采用美國MIT-BIH心律失常數據庫和歐洲ST-T標準心電數據庫中的數據進行算法驗證,首先統一抽樣頻率為200 Hz,進行Matlab仿真并與傳統的譜分析法比較,現將部分樣本數據仿真結果列出,如表 2所示。

對表 2中的Vtwa和Axial_ratio兩組數據進行相關分析[10],得互相關系數γ=-0.9676,說明Axial_ratio值越小,TWA幅值越大,且相關系數絕對值接近1,表明這兩個指標相關性很強,證明了本算法檢測TWA的有效性。再對其進行二次曲線擬合,得兩者之間的關系式為:
$\begin{align} & Axial\_ratio=-1.645\times {{10}^{-6}}\times {{V}_{TWA}}^{2} \\ & -7.233\times {{10}^{-4}}\times {{V}_{TWA}}+0.6361 \\ \end{align}$ |
Vtwa和Axial-ratio的二次曲線擬合如圖 3所示,其二次曲線擬合度為0.9443。

3 討論與結論
本文從圖形處理角度提出邊緣提取算法“橫向搜索法”,特點是將散點圖預處理為二值化的數字圖像矩陣,能夠提取散點圖最有價值的目標邊界信息,有選擇的處理T波采樣值,大量縮減待分析數據的同時剔除異常數據,算法簡單快捷易于計算機實現,明顯提高TWA檢測的實時性和準確性。其定量參數Axial_ratio是在散點圖“邊緣”的最小外接矩形上定義的,能夠恰當地描述散點圖形態與TWA的關系,其計算公式是基于數值分析的,采用線性規劃法尋找最優解,使數學計算最簡化,易于算法實現,大大提高計算機執行效率,進而有利于TWA的實時檢測。相比而言,譜分析法假定TWA分布在所有待檢測T波段上,無法檢測非平穩TWA信號,且檢測出的TWA幅值較真實值有所降低,再加上涉及譜分析,計算量較大[11]。散點圖具有直觀形態且計算簡單,由于是差值作圖,抗干擾性較強,對信號沒有高質量要求,不僅可以測得相鄰TWA幅值,而且包含時域信息,是一種新的可行的檢測方法。
目前,國內外關于散點圖法檢測TWA的論述和研究都比較少,該研究還處于初步階段,需要更多的數據分析來完善算法,其應用前景也有待于進一步臨床測驗。
引言
T波交替(T-wave alternans,TWA)是指在規則的心律時,體表心電圖上T波振幅、形態甚至極性的逐搏交替變化,它與惡性室性心律失常的發生密切相關[1]。TWA一般為微伏級,肉眼可分辨的TWA現象非常少[2]。根據TWA檢測統計方法不同,將TWA檢測分為短時傅里葉變換法、符號變換法和非線性法[3]。譜分析法是短時傅里葉變換法的一種,對64~128個ST-T結構中的T波采樣值進行快速傅里葉變換,將心電圖上T波段幅度的變化轉變成能量譜,分析0.5周期/心跳處的頻率成分,來測定是否存在TWA[4]。龐加萊映射[5](Poincare mapping,PM)是非線性法的一種,利用散點圖適于描述周期性變化數據的特性,可以使TWA中T波幅值逐拍變化的特點從形態上直觀表現出來。本文首先對散點圖進行圖形處理并定義散點圖短長軸之比為量化指標,通過與譜分析法的檢測結果進行相關性分析,證明散點圖法是檢測TWA的有效手段,并提出了一種新的TWA檢測有效量化指標。
1 算法介紹
1.1 心電信號預處理
首先通過簡單整系數50 Hz陷波器、零相位數字濾波器、結合小波分解與重構理論和改進的閾值算法去除工頻及其諧波干擾、基線漂移、肌電干擾等心電噪聲;然后采用二階微分Marr小波,應用多孔算法對ECG中QRS波群進行標定[6],在S波峰后150~350 ms的范圍內搜索極大值得到T波峰;隨后,對系列心跳進行復合預處理[7],提高TWA檢測的精確度。
1.2 TWA檢測
1.2.1 T波采樣點的選取
為了更準確的檢測微伏TWA(MTWA),我們采用N個心跳周期的改進T波窗口法。根據RR間期選擇T波窗口起點,通過Bazzet 公式確定T波窗口終點,T波選取按經驗原則,如表 1所示。

以T波波峰對齊,按T波平均長度的1/6為間隔,前后各取3個點構成7點檢測集xi={xi1,xi2,…,xi7},其中i代表所識別的心跳序列序數。本文選取128個T波段,每段提取7個采樣點,構成T波采樣點序列P={x1,x2,…,xi,…}。
1.2.2 散點圖法檢測TWA
對P={x1,x2,…,xi,…}進行一次差分組成新的序列Q={x2-x1,x3-x2,xi-xi-1,…},以X={x2-x1,x3-x2,…,xi+1-xi,…}為橫坐標序列,Y={x3-x2,x4-x3,…,xi+2-xi+1,…}為縱坐標序列,作Poincare散點圖,則該散點圖中點的坐標為(xi+1-xi,xi+2-xi+1)。本文選取128個T波段,每段提取7個采樣點,可作出882個散點。
若存在TWA,相鄰T波的采樣點幅值會以ABAB形式交替變化,即T波采樣點序列P={x1,x2,…,xi,xi+1,xi+2,…}的值會交替變化,則(xi+1-xi)與(xi+2-xi+1)異號且絕對值相當,也就是散點(xi+1-xi,xi+2-xi+1)主要以(+,-)或(-,+)的規律變化,此時散點圖呈長軸在直線y=-x上的傾斜橢圓形狀。若無TWA,相鄰T波的采樣點幅值呈AAAA形式,此時散點以原點為中心均勻分布,散點圖接近圓形。
經Matlab仿真繪圖,TWA散點圖外觀上有較明顯的邊緣[8],如圖 2所示。此邊緣是散點圖的最基本特征。為了方便對散點圖進行特征描述、分析,以及大量縮減待分析數據量,我們對TWA散點圖進行邊緣提取。根據其散點分布特點,研究出“橫向搜索法”來確定散點圖的“邊緣”,其具體步驟如下: ① 散點圖預處理。將散點圖轉換成M×N的二值化的數字圖像矩陣G。根據其橫、縱坐標序列X={x2-x1,x3-x2,…,xi+1-xi,…} 、Y={x3-x2,x4-x3,…,xi+2-xi+1,…}及散點坐標(xi+1-xi,xi+2-xi+1),首先將散點圖坐標序列值取整后上移和右移一定單位,使其剛好完全移至平面直角坐標系第一象限內,然后通過如下方式進行轉換:
$\begin{align} & G=\left[ \begin{matrix} {{g}_{11}} & {{g}_{12}} & \ldots & {{g}_{1N}} \\ {{g}_{21}} & {{g}_{22}} & \cdots & {{g}_{2N}} \\ \vdots & \vdots & \ddots & \vdots \\ {{g}_{M1}} & {{g}_{M2}} & \cdots & {{g}_{MN}} \\ \end{matrix} \right], \\ & (M=max\left[ Y\left( n \right) \right],N=max\left[ X\left( n \right) \right],n=1\,\!:L) \\ \end{align}$ |
其中Y(n)、X(n)分別表示Y、X中的第n個元素值;max[Y(n)],max[X(n)]分別表示X、Y中所有元素的最大值;矩陣元素gi,j表示圖像G中第i行,第j列的像素灰度值,其定義為
${{g}_{i,j}}=\left\{ \begin{align} & 1,\left( i=Y\left( n \right),j=X\left( n \right),n=1\,\!:L \right) \\ & 0,其他 \\ \end{align} \right.$ |
從而得散點圖的二值化數字圖像G。
② 尋找邊緣候選點。對圖像G進行塊操作,初步去除散點圖中心非邊緣點及邊界外較稀疏的非邊緣點。定義5×5的整數格,令整數格以gi,j,(i=1…M,j=1…N)為中心像素開始對散點圖進行行掃描,計算該整數格覆蓋區域A(見下式)的灰度值和a,即A中各元素之和,若2≤a≤3,則將該中心像素標記為邊緣候選點。
$A = \left| {\matrix{ {{g_{i - 2,j - 2}}} & {} & \cdots & {} & {{g_{i - 2,j + 2}}} \cr {} & \ddots & {} & {} & {} \cr \vdots & {} & {{g_{i,j}}} & {} & \vdots \cr {} & {} & {} & \ddots & {} \cr {{g_{i + 2,j - 2}}} & {} & \cdots & {} & {{g_{i + 2,j + 2}}} \cr } } \right|$ |
③ 確定邊緣點。對由上個步驟確定的邊緣候選點進行篩選,去除散點圖中心區域外的一些異常小點集,得到最終邊緣點。首先以邊緣候選點為中心像素構建21×21的整數格,則該整數格覆蓋區域B,即
$B = \left| {\matrix{ {{g_{i - 10,j - 10}}} & {} & \cdots & {} & {{g_{i - 10,j + 10}}} \cr {} & \ddots & {} & {} & {} \cr \vdots & {} & {{g_{i,j}}} & {} & \vdots \cr {} & {} & {} & \ddots & {} \cr {{g_{i + 10,j - 10}}} & {} & \cdots & {} & {{g_{i + 10,j + 10}}} \cr } } \right|,$ |
將矩陣B按B=劃分,得矩陣B1、B2、B3、B4,分別計算其灰度值和b,若所有b值都小閾值6,則刪除此候選點,將其余點保存在邊緣點數組中,得散點圖的邊緣點集。
為了更好地度量TWA的幅值大小,我們需要一種定量的檢測標準[9]。本文提出了短軸、長軸和軸比三個定量指標,并對軸比進行仿真測評。 我們設邊緣點集為Ω,(X,Y)表示Ω中點的坐標。以直線y=-x+z對Ω進行遍歷,當縱截距z(z=x+y)的值取最大和最小時,即分別找到該邊緣點集的兩條切線BC:y=-x+zmax和AD:y=-x+zmin。同理,以直線y=x+b對Ω進行遍歷,分別找到該邊緣點集的兩條切線AB:y=x+bmax和CD:y=x+bmin。由此得到該邊緣點集的最小外接矩形ABCD,散點圖幾何特征如圖 1所示。

① 定義短軸(short_axis)為最小外接矩形ABCD的短邊長,即:
$\begin{align} & Short\_axis=\left| AB \right|=\left| CD \right|=\frac{{{z}_{max}}-{{z}_{min}}}{\sqrt{2}} \\ & (z=X+Y,\left( X,Y \right)\in \Omega ) \\ \end{align}$ |
② 定義長軸(Long_axis)為最小外接矩形ABCD的長邊長,即:
$\begin{align} & Long\_axis=\left| BC \right|=\left| AD \right|=\frac{{{b}_{max}}-{{b}_{min}}}{\sqrt{2}} \\ & (b=X-Y,\left( X,Y \right)\in \Omega ) \\ \end{align}$ |
③ 定義軸比(Axial_ratio)為短軸與長軸的比值,這個指標能夠描述該邊緣點集沿直線y=-x的伸長程度。指標值越小,散點分布越趨近于直線y=-x;反之,散點分布越接近圓形,即:
$\begin{align} & Axial\_ratio=Short\_axisLong\_axis=\frac{{{z}_{max}}-{{z}_{min}}}{{{b}_{max}}-{{b}_{min}}} \\ & \text{ }(z=X+Y,b=Y-X,\left( X,Y \right)\in \Omega ) \\ \end{align}$ |
經Matlab仿真分析,規定Axial_ratio的閾值為0.57,即Axial_ratio≤0.57時,存在TWA;Axial_ratio>0.57時,無TWA。圖 2為TWA檢測的散點圖。

(a)不含TWA的散點圖;(b)含TWA的散點圖
Figure2. Poincare maps of TWA detection(a) the Poincare mapping without TWA; (b) the Poincare mapping with TWA
2 結果分析
采用美國MIT-BIH心律失常數據庫和歐洲ST-T標準心電數據庫中的數據進行算法驗證,首先統一抽樣頻率為200 Hz,進行Matlab仿真并與傳統的譜分析法比較,現將部分樣本數據仿真結果列出,如表 2所示。

對表 2中的Vtwa和Axial_ratio兩組數據進行相關分析[10],得互相關系數γ=-0.9676,說明Axial_ratio值越小,TWA幅值越大,且相關系數絕對值接近1,表明這兩個指標相關性很強,證明了本算法檢測TWA的有效性。再對其進行二次曲線擬合,得兩者之間的關系式為:
$\begin{align} & Axial\_ratio=-1.645\times {{10}^{-6}}\times {{V}_{TWA}}^{2} \\ & -7.233\times {{10}^{-4}}\times {{V}_{TWA}}+0.6361 \\ \end{align}$ |
Vtwa和Axial-ratio的二次曲線擬合如圖 3所示,其二次曲線擬合度為0.9443。

3 討論與結論
本文從圖形處理角度提出邊緣提取算法“橫向搜索法”,特點是將散點圖預處理為二值化的數字圖像矩陣,能夠提取散點圖最有價值的目標邊界信息,有選擇的處理T波采樣值,大量縮減待分析數據的同時剔除異常數據,算法簡單快捷易于計算機實現,明顯提高TWA檢測的實時性和準確性。其定量參數Axial_ratio是在散點圖“邊緣”的最小外接矩形上定義的,能夠恰當地描述散點圖形態與TWA的關系,其計算公式是基于數值分析的,采用線性規劃法尋找最優解,使數學計算最簡化,易于算法實現,大大提高計算機執行效率,進而有利于TWA的實時檢測。相比而言,譜分析法假定TWA分布在所有待檢測T波段上,無法檢測非平穩TWA信號,且檢測出的TWA幅值較真實值有所降低,再加上涉及譜分析,計算量較大[11]。散點圖具有直觀形態且計算簡單,由于是差值作圖,抗干擾性較強,對信號沒有高質量要求,不僅可以測得相鄰TWA幅值,而且包含時域信息,是一種新的可行的檢測方法。
目前,國內外關于散點圖法檢測TWA的論述和研究都比較少,該研究還處于初步階段,需要更多的數據分析來完善算法,其應用前景也有待于進一步臨床測驗。