基于超聲造影的頸部淋巴結灌注流向可視化①

李明豪, 尹 皓, 趙海娜, 劉東權

1(四川大學 計算機學院, 成都 610065)

2(四川大學華西醫院 超聲醫學科, 成都 610041)

頸部淋巴結會受甲狀腺病變細胞轉移的影響而發生病變. 發生腫瘤轉移的淋巴結, 腫瘤細胞會通過淋巴管首先聚集于邊緣竇, 繼而逐漸侵入淋巴結實質; 其可破壞淋巴結正常的動靜脈結構, 誘導新生血管的形成. 腫瘤細胞的轉移會影響頸部淋巴結的血流情況. 基于超聲造影(contrast-enhanced ultrasound, CEUS)可以通過造影劑判斷血流灌注的信息, 輔助醫生進行病情診斷[1].

傳統的灌注曲線(時間強度曲線, TIC)是基于病灶區域每幀取強度值的平均值而沿時間軸成時間序列而得. 然而, 由于感興趣區域(region of interest, ROI)的灌注強度可能并不均勻, 這樣將使提取的特征不準確[2].因此, 本文對感興趣區域進行像素級分析, 即視頻的一個坐標將沿時間軸擬合一條時間強度曲線. 所有感興趣區域內的曲線對一個灌注特征(達峰時間、峰值等)進行提取可以形成一個二維矩陣. 二維矩陣可以進行編碼顯示成圖, 也可以進一步提取更深層次的信息來進行可視化.

本文主要的創新實驗是, 基于達峰時間(time to peak,TTP)二維矩陣進行研究, 猜測感興趣區域中相距不遠的達峰時間接近的兩坐標可能存在流向先后的關聯性并因此尋找血流流向的信息, 來進行灌注流向的成像可視化.

1 相關技術與研究現狀

1.1 研究現狀

探頭固定時, 頸部淋巴結采集的造影視頻, 目標跟蹤過程是不考慮尺度變化的. 比如, Zhang等人在2011年認為頸部淋巴結在探頭不動的情況主要會受呼吸作用影響而做類似周期性的運動, 而直接利用基于塊匹配的方法來進行目標跟蹤[3]. 其利用了最小化絕對差和(sum of absolute difference, SAD). 基于光流法也能進行目標跟蹤, 但是對光照要求太過嚴格[4]. Ta等人在2012年利用了基于互信息的方法來進行目標跟蹤[5].然而互信息始終計算量太大, 并且也不能很好地處理異常幀. 而后, 林細林等人2018年利用了感知壓縮跟蹤[6], 該方法通過改進可以處理異常幀的情況[7], 但壓縮感知跟蹤的跟蹤結果位置不夠精確, 甚至位置還會跑得很遠[8]. 由此本文實驗了KCF (核相關濾波)算法來進行目標跟蹤, 其利用了相關濾波理論與核函數技術, 使得目標跟蹤很魯邦而且計算量不會很大[9].

時間強度曲線可以反映血液中的造影劑濃度變化.Zhang等人早期就對感興趣區域提取亮度灰度值的平均值, 隨幀數生成時間強度曲線(time intensity curves,TIC), 提取出達峰時間、峰值(peak value, PI)等灌注參數來分析視頻[3]. Gaia分析類風濕性關節炎時, 指出可以利用像素級的時間強度曲線, 即每個坐標點擬合一條曲線. 然后對每一個灌注參數可以生成一個二維矩陣, 可以編碼成像來幫助醫生診斷病情. 其擬合曲線時利用的是Gamma公式[10]. 王本剛等人分析肝臟病例也進行了基于像素的時間強度曲線分析, 不過只是利用了savgol濾波器進行曲線擬合[2]. 這些文獻沒有研究怎么樣將沒有很好灌注趨勢的曲線篩除, 本文將對此進行研究.

針對灌注流向問題, 本文嘗試利用流線與向量場來進行實驗. 自Brian 等提出線積分卷積 (line integral convolution, LIC)的方法后, 對向量場進行可視化大部分是利用了LIC方法及其變體[11]. 雖然Wegenkittl 等人提出 Oriented LIC 方法來使結果紋理含有方向指向信息[12]. 但本文病例感興趣區域太小, 因此不能用該方法表示方向. 本文使用了Brain提出的基本的LIC方法, 同時利用顏色編碼來表示坐標的先后關系.

1.2 相關技術

第1步. 對淋巴結進行目標跟蹤.

實驗中所采集的超聲造影視頻都為雙模式并排顯示, 如圖1所示, 左為超聲灰階B模式, 右為超聲造影模式. 目標跟蹤在灰階B模式上面進行, 然后再將位置信息應用到彩色造影模式中以分析灌注情況.

圖1 雙模式視頻幀

由于淋巴結主要是因為呼吸會產生運動, 但沒有明顯尺度變化. 所以可以由醫生進行感興趣區域的標記, 然后提取出掩膜mask以便在視頻中進行目標跟蹤, 如圖2所示.

圖2 提取mask掩膜示意圖

這里手工標定的幀由醫生選取, 被用來做目標跟蹤的起始幀, 被稱為“參考幀”. 參考幀不一定是原視頻的第1幀.

本文選擇核相關濾波(kernelized correlation filters,KCF)[9]和Kalman濾波結合的方式進行目標跟蹤[13].

第2步. 進行像素級的時間強度曲線擬合.

基本思路是: 在B模式對淋巴結進行目標跟蹤后,將位置信息用到造影模式. 造影模式ROI每一個坐標沿著時間軸都能擬合一條時間強度曲線. 本文利用Gamma公式模型, 同時也研究怎么篩選符合灌注特征的曲線.

第3步. 提取灌注特征參數. 因為是基于像素的分析, 所以對符合標準的曲線提取灌注參數(達峰時間、峰值等), 形成二維矩陣. 對二維矩陣可以編碼成像顯示, 輔助醫生診斷病情.

第4步. 利用TTP參數二維矩陣可視化灌注流向情況. 設種子點為TTP值低于設定閾值的坐標點. 先找種子坐標點, 每個種子點會找一條流線, 這些流線可以進行可視化. 實驗也嘗試了利用流線信息生成向量場,并利用線積分卷積的方式可視化.

2 病灶跟蹤

2.1 淋巴結關于目標跟蹤方面的特點

探頭固定的情況下, 視頻灰階B模式相鄰兩幀之間的運動變化很小可以忽略, 整體運動以呼吸運動的影響為主, 而呼吸運動可以看成周期運動[3]. 在視頻中淋巴結很少發生巨大的尺度變化(即縮放變化), 也不會有大的形變(指形狀改變成另外的樣子).

灰階B模式采用的是基波, 造影模式采用的是諧波. 造影劑的濃度變化會使彩色造影部分的亮度跟著改變. 這也會對灰階B模式產生影響, 但亮度不會發生很大變化.

視頻中會出現短暫的異常幀, 持續時間很短, 并且之后目標形態會恢復正常. 異常幀常是由于信號噪聲、病人的較強的呼吸引來的其他平面信息等因素而導致的, 主要表現是目標細節暫時丟失.

總結下來, 沒有異常幀時可以利用簡單目標跟蹤方法(如基于灰度的模板匹配)來處理. 有異常幀時就要考慮如何檢測與處理異常幀.

2.2 目標跟蹤模型

傳統的基于灰度的模板匹配(如最小化絕對差和SAD、最小化平均絕對差MAD、最大化歸一化互相關NCC等)進行目標跟蹤, 一幀一幀更換參考模板的可以漸變地跟蹤目標. 但是, 其并沒有提出好的檢測異常幀的方法.

基于要處理異常幀的思路, 本文選擇了相關濾波的方法進行目標跟蹤.

相關濾波的思想是利用信號的相關性, 設計一個濾波模板, 利用模板與目標候選區做相關運算生成響應圖, 響應圖上值越大的位置越可能是待搜索幀的目標位置. 相關濾波跟蹤的優點是可以利用卷積計算在傅里葉域的性質減少計算量.

用KCF (核相關濾波)進行目標跟蹤, 利用HOG特征會對光照變化有魯邦性, 異常幀檢測可以利用響應圖的峰值旁瓣比(peak to sidelobe ratio,PSR)[9].

峰值旁瓣比公式為[14]:

其中,gmax為相關響應圖最大峰值, 峰值周圍半徑為11像素大小的圓形區域之外為旁瓣, μs與δs分別為旁瓣的均值與方差.

2.3 卡爾曼濾波

卡爾曼濾波器是一種高效的遞歸濾波器, 能夠從一系列不完全包含噪聲的測量中, 估計出動態系統的狀態[13]. 其常用在軌跡預測與目標定位中. 對于一般的線性隨機系統, 可表示為:

其中,xk∈Rn表示狀態向量,zk∈Rm表示測量向量,與Gk表示狀態轉移矩陣,Hk表示測量矩陣,wk是過程噪聲,vk是 測量噪聲.wk的協方差矩陣為Qk,vk的協方差矩陣為Rk. 實驗中假定Qk與Rk均為零均值高斯白噪聲, 且兩者互不相關.

上述模型是線性的Kalman濾波模型, 常有5個基本公式. 用Kalman濾波進行數據的濾波, 可以從后驗結果得到最終的數據濾波結果, 比如濾波后的曲線數據; 用Kalman濾波進行預測, 可以從先驗結果里得到預測結果, 比如預測的位置信息.

2.4 總的目標跟蹤思路

視頻中大部分時間淋巴結的運動速度都是較慢的,但突發情況時會發生較快的運動. 較快的運動就會有一些幀中目標變得模糊, 或產生偽影. 因此建立線性的卡爾曼模型來檢測目標是否運動過快, 一旦速度超過閾值, 就不對KCF濾波模板進行更新[7].

最終, KCF的模板大小為長寬96 像素. 如果圖像太小, 模板邊長要適當變小. PSR的最小閾值為40,Kalman預測位置與相關濾波檢測位置的歐式距離最大閾值為60 像素.

如果上述兩閾值條件都滿足, 才更新KCF的濾波模板. 對于Kalman模型的參數更新, 如果某幀不滿足PSR最小閾值就用Kalman預測位置更新; 如果滿足PSR最小閾值, 就用KCF檢測位置更新. 目標跟蹤中的更新操作流程如圖3所示.

圖3 目標跟蹤中的更新操作流程

3 進行像素級的時間強度曲線擬合

3.1 二維平滑窗口的選擇

時間強度曲線是每個坐標都沿時間軸提取一條曲線. 但是會出現大量坐標的數據不能曲線擬合成功.

針對該問題, 猜測淋巴結內部各結構會隨著呼吸相對于淋巴結整體發生很細微的運動. 考慮利用一個二維平滑窗口, 將呼吸引起的細微運動的軌跡籠罩住,窗口內提取平均值. 平均值的變化就能符合灌注趨勢了. 這里平滑窗口高寬是“29×29”.

3.2 TIC曲線數據的提取

TIC (time-intensity curve)曲線即時間強度曲線. 將造影圖片從RGB顏色空間轉換為“HSV”顏色空間[15].生成H分量的直方圖之后發現數據分布集中在很窄的一段, 如圖4, 即造影部分色調都比較接近. 因此, 圖像的色調不存在太多信息, 飽和度又不能說明亮度變化.最終選擇了“HSV”的“V”分量(縮放到[0, 255])來表示“強度”.

圖4 一幀造影與其HSV空間的H分量的直方圖

每一坐標點可以沿時間軸得到一條TIC曲線的原始數據. 當TIC的原始數據得到之后, 會首先采用S-G平滑濾波(Savitsky-Golay濾波器)[2], 多項式階數是2,窗口大小是81 (原視頻幀率是10 fps, 所以這里代表8.1 s). 濾波目的是減弱噪聲.

3.3 基于Gamma公式的曲線擬合

Gamma模型公式[10]為:

其中,A≥0,C＞0,k＞0 ,B是初始值即t=0時刻的y值.B為固定值, 不屬于可變參數,AT(arrival time)為造影劑到達時間.

利用Gamma公式進行曲線擬合, 可以描述造影劑在血管中的洗入(wash-in, 表現為曲線上升)與洗出(wash-out, 表現為曲線下降)的時間段情況, 便于提取相關的灌注特征參數, 如圖5所示.

圖5 原數據和Gamma擬合圖

3.4 曲線數據的篩選之粗篩選與細篩選

并非ROI內所有的點都能很好地擬合曲線, 考慮原因是并非所有區域都是液態的有血流經過的區域.因此, 需要將不符合的曲線數據篩除.

這里考慮用“粗篩選-細篩選”的方式, 只有兩種篩選都通過的曲線才能用于以后的實驗.

第1步, 粗篩選.

粗篩選大致分兩個思路, 兩個子篩選條件都滿足了才算通過粗篩選. 一個是真實值的判斷, 一個是縮放值的判斷.

(1)真實值判斷, 就是對曲線原數據的峰值、起始值判斷大小.

比如在實驗中, 峰值不能小于30, 曲線起始值不大于100. 這里注意由于利用HSV的方法, 所有曲線中強度值范圍可統一為[0, 255], 如圖6所示.

圖6 真實值判斷示意圖

(2)縮放值判斷, 將曲線原數據值縮放到[0,n],n為正整數, 這里為8. 對數據中的第2峰值、后面(1/2時間之后)最大值、起伏過程中的上升幅度、起伏過程中的異常上升總次數等進行判斷.

比如實驗中, 第2峰值在數據后面不能大于5, 在前面不能相比之前的波谷上升超過3. 后端1/5時間段最大值不能是6. 上升幅度大于3為異常上升, 異常上升次數不能多于3個, 如圖7所示.

圖7 縮放值判斷示意圖

第2步, 細篩選.

滿足粗篩選的曲線再次用細篩選進行處理. 細篩選是利用Gamma公式去擬合數據, 根據擬合情況對數據進行最后篩選.

因為“R-squared”一般應用于線性擬合中的擬合度評價準則, 所以沒有使用在這里[16]. 本文的擬合度的評價準則如下:

(1)NRMSE, 即歸一化的RMSE(均方根誤差).

這里RMSE公式為:

其中,ri是 原始數據,fi是擬合后的數據.

而NRMSE的公式為:

NRMSE的最大閾值為0.15.

(2) Adjusted Cosine即修正的余弦相似性[17]. 公式為:

其中,u即原數據r與擬合數據f合在一起共同的平均值. 使用修正余弦相似性有一個重要條件, 需要兩個時間序列一樣長.

考慮到擬合的曲線可能會直接平滑掉起伏頻率高的數據段, 可能出現NRMSE小但是跟原數據根本就不像的情況. 因此, 利用余弦相似度進行衡量, 這里設置的最小閾值為0.75.

4 提取灌注特征參數

4.1 提取2-D灌注參數特征

對于一條TIC曲線可以提取灌注參數, 而參數大致分兩類[5]: (1)基于時間的, 如TTP (達峰時間)、MTT(平均渡越時間); (2)基于數值的, 如AUC (曲線下面積)、PI (峰值).

在視頻中, 像素級提取時間強度曲線, 則對于一個灌注特征參數可得到一個二維矩陣, 可稱為“灌注參數二維矩陣”.

另外, 曲線擬合成功與失敗的坐標需要記錄下來,下文需要使用.

4.2 灌注參數二維矩陣的可視化

曲線擬合成功的點才能進行正常的顏色編碼來可視化, 失敗的坐標點將會設置為黑色.

首先會將參數二維矩陣中成功的坐標的值統一縮放到[0, 240]之間. 因為在HSV顏色空間中, 色調分量H從“0到240”會對應“紅到藍”. 將S和V統一設置為1 (或滿值), 不同坐標點只需設置各自的H分量, 便可得到一張彩色圖.

然后, 對縮放后的二維矩陣進行編碼可顯示成圖片. 不過, 顏色編碼方式需結合實際情況指定.

比如, 對于TTP (達峰時間)來說, 達峰時間越小就說明達峰時間早, 越早越是醫生感興趣的越偏向暖色(紅色). 此這一類數據采取“藍到紅, 大到小”的編碼方法, 如圖8.

圖8 TTP二維矩陣成像圖

而對于AUC (曲線下面積), 值越大說明灌注總量越大, 越表示醫生感興趣的而偏暖色(紅色). 因此這一類數據采取“藍到紅, 小到大”的編碼方法, 如圖9.

圖9 AUC二維矩陣成像圖

5 基于達峰時間的灌注流向估計與可視化

5.1 灌注方向的求法

灌注大致方向對于病情的判斷是有幫助的. 所以,為了可視化灌注方向, 本課題利用了達峰時間二維矩陣進行估計. 總體思路是, 造影會從TTP小的地方流向TTP大一點的地方.

5.2 流線種子點的選取

針對達峰時間二維矩陣, 設置一個TTP最大閾值,將低于最大閾值的區域提取出來作為種子點.

對于TTP二維矩陣查找最大值maxValue和最小值minValue, 然后閾值Th為:

其中,frac為“0”到“1”之間的小數, 這里frac設置為0.3.

5.3 每條流線的生成

流線是一條其上任意一點的矢量方向都與流線相切的曲線. 流線方程定義如下[18]:

其中, μ(τ) 代 表空間中點的位置,v(μ(τ))為流線上該點的切線方向,v是τ 的函數, τ可以是時間或弧長等參量.

本文一條流線可以利用一組具有先后順序的采樣點來確定. 這里利用鏈表的形式給出, 為了理解, 這里借用C++的容器表示.

class Point{

int x;

int y;

};

std::vector＜Point＞ streamlineList; //即一條流線的點組成的list.

而這里尋找流線的思想: 每個點以自己為起點找下一個符合的點, 并以找到的點為起點再找下一個點.如圖10(a), 一條流線的例子, 從A點開始, 直到F點結束就形成了一條流線的采樣點.

圖10 尋找流線示意圖

每次一個點為起點尋找下一個點時, 找大于該點且時間差為“0.2”或“0.3”秒的最小值(兩個點不能同時發生, 所以時間差要有最小閾值), 搜索范圍并不是采用固定的8鄰域, 而是以當前點為中心, 一層一層(1、2到3的順序)地搜尋外圍的鄰域點, 直到找到一個TTP大于該點且超過一定數值的最小值.

如圖10(b), 以當前點P搜索下一個候選點的過程如下:

(1)設當前參考點為P, 并設置預定標準為下一點是該點值是比當前參考點大0.2的點的最小值.

(2)搜索P點的8鄰域點, 如圖中標記為“1”的那些點. 如果找到預定標準的候選點就結束; 否則將當前點和其8鄰域點加入一個點群pointGroup, 然后進行第(3)步;

(3)繼續搜索pointGroup點群的4鄰域點(比如上次搜索標記為“1”, 這次搜索標記為“2”的那些點). 當前參考點依舊為中心點P (即第(1)步的點). 如果找到預定標準的點就結束; 否則就將其周圍點加入pointGroup點群, 重復本步驟. 另外, 此步驟中pointGroup點群的長度大于限制值或ROI內所有點都已經遍歷也會終止.

5.4 流線顯示

找到每條流線采樣點中, 一條流線上每相鄰的兩個點可能在圖上空間不相鄰, 可以利用Bresenham直線算法用一條近似直線的像素點連接兩點[19].

可問題是, 方向怎么表示. 圖像小, 不考慮用箭頭表示流向, 而利用TTP參數成像圖里的顏色(色調)不同進行標記.

另外, 對區域中每一個坐標點統計有多少條流線經過(遍歷所有曲線得出), 可以生成一個流線計數二維矩陣, 流線經過比較多的區域有很大可能是有血管的區域. 由于流線計數圖亮度低的區域的細節看不清.可以利用log增強技術來提高低灰度區域的對比度.

流線圖、流線計數圖以及流線計數圖log增強圖都顯示在圖11中.

圖11 流線結果顯示

6 生成向量場并可視化

6.1 生成向量場

不同的種子點會有不同的流線, 可生成大小與原參數成像圖相同的向量場, 每個坐標點有一個向量, 為了表示每一個坐標點的造影劑在下一時刻的可能去向[18].

一條流線可以利用Bresenham直線算法填充空隙,然后從前到后依次指向, 每個非終止點可以指向流線上下一點而成向量. 而每一條流線的終點會設置成零向量.

可以用一張mask來表明某坐標是否已經因為一條流線經過而被賦值過. 遍歷每一條流線, 然后在向量場相應位置修改向量值. 只有流線經過的點有非零的向量, 其他點向量都為零向量. 當出現某個點有多條流線經過的時候, 選取流線計數大的下一個點來作為下一個最優指向.

6.2 線積分卷積

LIC技術根據向量場, 用一維濾波器沿流線卷積白噪聲圖像, 合成紋理. 該紋理與向量場的向量方向高度相關.

其生成過程需要輸入向量場和一張噪聲圖(通常是白噪聲)[11], 如圖12所示.

圖12 LIC圖生成示意圖

其基本思想是, LIC會選擇噪聲圖像作為輸入背景紋理, 然后對每一個坐標點正反向去各追蹤一條流線, 流線上的每一個像素點根據卷積核算出的權重來進行卷積計算, 得出紋理值. 而采用不同的核函數, 加之不同的背景噪聲, 最終得到的效果也不一樣[20], 本實驗選擇的是盒型核與白噪聲.

6.3 灌注流向的向量場成像

生成LIC圖時要對零向量區域進行一步抑制, 縮小其亮度范圍, 利用了TTP成像圖的顏色來著色表示坐標點的先后關系(越接近暖色即紅色, 則達峰時間越小).

如圖13, 是對上面流線圖的例子進行計算的結果.其紋理顯示了向量場的大致狀態.

圖13 LIC結果圖

7 實驗結果分析

7.1 數據

本課題所選擇的病例視頻由四川成都華西醫院提供, 造影劑采用意大利博萊科公司生產的第二代造影劑聲諾維. 所有病例來自于同一個邁瑞的機器. 采集期間, 探頭固定, 囑患者固定體位, 平穩呼吸, 勿做吞咽動作. 目前實驗數據共18例病例, 12例良、6例惡.

7.2 曲線擬合篩選后的有效點面積占比

如表1所示, 18例視頻都進行了實驗. 有效坐標點是上文曲線擬合成功坐標點. 有效坐標點占比是有效坐標點數占ROI內總點數的比重.

表1 有效坐標點占比結果

對于造影視頻, 符合灌注趨勢的坐標應是有血流經過的坐標. 有效坐標點占比太小則合格血管的覆蓋面積太少.

因為本文的重點是探尋灌注流向信息可視方法,所以實驗的數據不是很多. 但在已有的數據中, 大部分良性視頻的有效成功點數占比為“50%–100%”, 而有一半的惡性視頻的有效成功點數占比為“0–50%”. 如果發生了有效成功點占比低于50%的事情, 樣本極有可能是惡性的, 也可能是采集異常, 可標記為“可疑的”.

7.3 4種不同面積占比的視頻的灌注可視化結果

如圖14, 觀察不同的擬合成功點占比下灌注情況可視化結果. 其中, 每一行圖片為一個視頻的結果, 從上而下有效坐標占比依次為“0–20%”“20%–50%”“50%–80%”“80%–100%”. 每一行從左到右分別為TTP成像圖、流線圖、流線技術圖(log增強)、LIC結果圖.

圖14 4種不同有效點面積占比的結果

當有效坐標點占比過于小時(低于50%時), 有關灌注情況的效果會很差. 因此, 建議觀察流向信息時,有效坐標占比要超過50%.

7.4 總的結果分析

當確定ROI內基本是血管時, 可以將有效坐標點占比小(比如低于50%)的視頻標記為“可疑”視頻供醫生分析. 有效坐標點(擬合曲線成功的坐標)占比太小,不建議進行灌注可視化. 利用上述方法尋找到的流線,流線圖與流線計數圖非常類似于血管的分布, 應具有一定的研究價值.

8 結論與展望

本文研究了對造影視頻的灌注特征的提取與可視化. 提取灌注參數(達峰時間、峰值等)來編碼成像可以協助醫生分析病情, 使用流線等方法可視化造影劑在ROI內的流向信息. 結果中流線圖與流線計數圖中的線條類似血管, 應具有一定的研究意義.

經歷了目標跟蹤、時間強度曲線擬合、灌注參數提取、根據達峰時間二維矩陣進行灌注情況可視化共4個步驟. 目標跟蹤根據視頻的特點而選擇了KCF與Kalman濾波結合的方式. 時間強度曲線擬合是像素級的, 經歷了粗細篩選的方式, 將符合灌注特征的曲線留下來. 灌注參數提取主要是提取達峰時間、峰值等. 像素級的曲線擬合與提取參數可以形成一個二維矩陣,可顏色編碼成像顯示, 也可做一些處理提取特征. 根據達峰時間二維矩陣, 做了流向情況可視化, 最后可得到流向可視化.

做流向情況可視化時, 有效坐標占比要大于50%才能有較好的可視化效果.

未來更進一步的實驗將是: 更好的曲線擬合公式與篩選方式, 因為目前Gamma公式可能并非對灌注情況描述最佳的公式, 雖然達峰時間誤差小但峰值處并不能很好地貼合原數據; 根據參數二維矩陣更好的尋找流線的方式, 目前尋找流線的方式還需要進行更進一步的標準; 使用更多的數據, 只有數據量足夠才能找到更精確的標準.