第四節(jié)　有限元法心臟運(yùn)動(dòng)建模

基于有限元方法的模型更側(cè)重于分析心室各個(gè)內(nèi)部單元的應(yīng)力應(yīng)變分布情況，在局部運(yùn)動(dòng)分析方面有限元模型比可形變模型更細(xì)致精確。之所以采用四面體單元?jiǎng)澐中呐K表面實(shí)體模型，是因?yàn)橛盟拿骟w單元?jiǎng)澐秩S實(shí)體很靈活，可以逼近較復(fù)雜的幾何形狀，并按需要在局部區(qū)域隨意加密或變稀網(wǎng)格，而心臟的幾何形狀比較復(fù)雜、不規(guī)則，適合用四面體單元?jiǎng)澐帧?/p>

一、心臟有限元模型的應(yīng)變分析

心肌是近似不可壓縮的，心肌被視為是一種非剛性的會發(fā)生形變的彈性體材料，心臟在收縮、舒張時(shí)必然會產(chǎn)生體力或表面力；而應(yīng)變分析是描述連續(xù)體內(nèi)部形變的方法，是目前用于量化分析心臟形變的常用工具。許多心臟運(yùn)動(dòng)分析算法都使用應(yīng)變張量來表示心肌的形變，應(yīng)變張量是彈性力學(xué)中的概念，能包括形變的全部信息。該方法能對心肌組織的形變進(jìn)行精確估計(jì)，因此被廣泛應(yīng)用于描述彈性塑性體的形變。在有限元模型中，重建心臟在整個(gè)心動(dòng)周期中的3D形狀以及任意心肌點(diǎn)的運(yùn)動(dòng)軌跡，并能對心肌的運(yùn)動(dòng)情況進(jìn)行分析，是目前心臟圖像分析領(lǐng)域中重要的研究方向。

圖102為初始時(shí)刻的心臟表面有限元模型，全部由四面體單元組成，該模型由619個(gè)四面體單元組成，共771個(gè)節(jié)點(diǎn)。運(yùn)用四面體單元剖分三維實(shí)體時(shí)，可以更靈活地進(jìn)行剖分，逼近復(fù)雜的幾何形狀，也可以按需要在局部區(qū)域隨意加密或變稀網(wǎng)格。通常在應(yīng)力集中區(qū)域或應(yīng)力梯度高的區(qū)域布置較密的網(wǎng)格，在應(yīng)力變化平緩的區(qū)域布置較稀疏的網(wǎng)格。這樣做可以同時(shí)滿足精度和效率兩方面的要求。

根據(jù)所建立的基于3D表面重建的動(dòng)態(tài)有限元方程的模型以及基于形變模型的約束條件，就可得到整個(gè)結(jié)構(gòu)的所有節(jié)點(diǎn)的位置矢量，由形態(tài)函數(shù)可得到表面上任一點(diǎn)的位置?？紤]到心肌運(yùn)動(dòng)應(yīng)該是連續(xù)的彈性運(yùn)動(dòng)，本書運(yùn)用薄板樣條對模型受力進(jìn)行了優(yōu)化，針對本書建立的四面體元心臟表面三維有限元模型，關(guān)鍵是導(dǎo)入邊界條件求得作用在各單元節(jié)點(diǎn)上的節(jié)點(diǎn)力，計(jì)算模型中各單元節(jié)點(diǎn)的位移。已知位移，計(jì)算單元節(jié)點(diǎn)的應(yīng)變，其中某一點(diǎn)的應(yīng)變狀態(tài)可用六個(gè)獨(dú)立的應(yīng)變分量表示：

其中：ε_x，ε_y和ε_z是正應(yīng)變，γ_xy，γ_yz和γ_zx是切應(yīng)變。用u，v，w表示四面體各節(jié)點(diǎn)沿x，y，z方向的位移，這些位移一般應(yīng)為各節(jié)點(diǎn)坐標(biāo)的函數(shù)，即

為了反映某節(jié)點(diǎn)總的應(yīng)變狀態(tài)，把公式（1016）中6個(gè)應(yīng)變分量用米澤斯屈服條件轉(zhuǎn)換為米澤斯應(yīng)變，即：

所以，如何求左心室表面各單元節(jié)點(diǎn)的位移是求解心臟表面應(yīng)變分布的關(guān)鍵。

心臟表面的幾何形狀比較復(fù)雜、不規(guī)則，適于用四面體單元?jiǎng)澐?，在所建立的有限元模型中，目前非剛體運(yùn)動(dòng)對應(yīng)點(diǎn)估計(jì)方法已成功地運(yùn)用于心臟心室的表面運(yùn)動(dòng)跟蹤以及大腦皮質(zhì)表面點(diǎn)匹配。

UNGCK算法是一種結(jié)合了單位法線（unit normal，UN）算法和高斯曲率（Gaussian curvature，GC）算法各自優(yōu)點(diǎn)的非剛體運(yùn)動(dòng)對應(yīng)點(diǎn)估計(jì)方法，而且已應(yīng)用在測量犬心臟的左心室體積膨脹系數(shù)上。有限元方法和UNGCK算法結(jié)合分析心臟表面應(yīng)變，其適用于心臟三維實(shí)體的有限元單元網(wǎng)格模型。

在曲面造型中，曲面在任一點(diǎn)附近的形狀與在該點(diǎn)曲面的主曲率的乘積即高斯曲率有關(guān)。該點(diǎn)與附近點(diǎn)的高斯曲率比較可以反映出該點(diǎn)附近的形狀變化，故可以用高斯曲率表達(dá)該點(diǎn)的形狀信息，對該點(diǎn)附近的形狀質(zhì)量進(jìn)行評判。高斯曲率計(jì)算原理圖圖10-3所示。三維空間的3點(diǎn)確定1個(gè)平面，在這一平面內(nèi)的這3點(diǎn)確定1個(gè)球，即確定該球的半徑和沿著球的切矢。設(shè)表示曲面S在一點(diǎn)P上的單位法矢，切S且經(jīng)過的平面與曲面相交成一條曲線。同樣，不經(jīng)過但經(jīng)過P點(diǎn)的平面與曲面同樣也可以相交成一條曲線。讓每一個(gè)法平面與一個(gè)方向及單位切矢對應(yīng)，即在曲線P點(diǎn)，一個(gè)法曲面曲率k_n，對應(yīng)一個(gè)位置。這個(gè)法曲面曲率隨著切平面繞的旋轉(zhuǎn)而變化。k_n存在最大和最小值，即為P點(diǎn)的主曲率。

圖10-3　高斯曲率計(jì)算原理圖

UNGCK算法具有較高的準(zhǔn)確率和較好的數(shù)值穩(wěn)定性，但要直接用在心臟有限元模型上還應(yīng)做相應(yīng)的改進(jìn)。我國學(xué)者劉復(fù)昌、夏德深在相關(guān)的論文中提出了一種基于左心室有限元模型的α剪枝平均UNGCK算法。論文的α剪枝平均UNGCK算法考慮到左心室有限元模型表面由許多三角形拼成，拓?fù)潢P(guān)系不規(guī)則，因此取在表面上與節(jié)點(diǎn)鄰近的一些點(diǎn)代替節(jié)點(diǎn)鄰域擬合曲面。在網(wǎng)格節(jié)點(diǎn)所在的擬合曲面上，以該網(wǎng)格節(jié)點(diǎn)為中心的鄰域內(nèi)重新找運(yùn)動(dòng)誤差最小的點(diǎn)作為對應(yīng)點(diǎn)，這樣找到的對應(yīng)點(diǎn)不再局限于網(wǎng)格節(jié)點(diǎn)。采用經(jīng)典的剛體配準(zhǔn)算法——最近點(diǎn)迭代（iterative closest point，ICP）算法來估計(jì)相鄰時(shí)刻對應(yīng)點(diǎn)的候選點(diǎn)，在一對多的對應(yīng)關(guān)系下，它選取距離最近的且該點(diǎn)法線方向與待匹配點(diǎn)法線方向夾角最小的點(diǎn)作為對應(yīng)點(diǎn)。α剪枝平均的好處就是可以去除數(shù)據(jù)集中的一些不正常點(diǎn)（即過大或過小點(diǎn)），而對正常的點(diǎn)改變很小，這樣可避免把不正常點(diǎn)以相同的權(quán)值和正常點(diǎn)作平均，比直接對數(shù)據(jù)點(diǎn)作平均要合理得多。他們提出基于有限元和α剪枝平均UNGCK算法的左心室應(yīng)變分析，利用α剪枝平均對UNGCK算法存在的局部運(yùn)動(dòng)函數(shù)與物體整體運(yùn)動(dòng)不一致的情況作修正，能直觀、有效地反映左心室在整個(gè)心臟收縮期的形變情況及應(yīng)變分布。

二、有限元模型的心臟運(yùn)動(dòng)建模

如何求心臟表面各單元節(jié)點(diǎn)的位移是求解心臟表面應(yīng)變分布的關(guān)鍵。顯然四面體元有12個(gè)自由度，每個(gè)節(jié)點(diǎn)有3個(gè)自由度，因此，對一個(gè)有n個(gè)節(jié)點(diǎn)的結(jié)構(gòu)而言，其整體剛度矩陣K是3n×3n的，可求解δ的結(jié)構(gòu)節(jié)點(diǎn)位移矢量，F(xiàn)的結(jié)構(gòu)節(jié)點(diǎn)載荷矢量。一旦求得未知的位移與節(jié)點(diǎn)力矢量，就可以用公式（1019）求得每個(gè)單元的應(yīng)力矢量。

公式（1019）中，σ是6×1的單元應(yīng)力矢量，δ是12×1的單元節(jié)點(diǎn)位移矢量，每個(gè)單元的σ形式為：

其中前三項(xiàng)為正應(yīng)力，后三項(xiàng)為切應(yīng)力，在實(shí)際的程序設(shè)計(jì)實(shí)驗(yàn)中，這些參數(shù)的選擇將根據(jù)要解決的問題進(jìn)行調(diào)整和取值。

在建立基于有限元方法的心臟實(shí)體模型中，初始時(shí)刻的心臟表面有限元模型全部由四面體單元組成，該模型由619個(gè)四面體單元組成，共771個(gè)節(jié)點(diǎn)。在實(shí)際進(jìn)行的程序檢測中，我們發(fā)現(xiàn)771個(gè)節(jié)點(diǎn)中的12個(gè)點(diǎn)是在心臟表面有限元模型的上邊緣，形成鋸齒，它們即所謂的孤立點(diǎn)，有的形成的四面體體積V為零的，有的沒有成為在我們所生成的四面體中的單元節(jié)點(diǎn)。因此必須把這些點(diǎn)先處理好，使其不得參與后面的有限單元的計(jì)算，以免影響求位移與應(yīng)力的正確性。通過應(yīng)用程序我們檢測到如果有四面體單元的體積為零的四面體，必須排除出去，否則會出現(xiàn)分母為零時(shí)得到無窮大的結(jié)果，顯然這是不合理的。經(jīng)調(diào)試程序時(shí)發(fā)現(xiàn)生成的四面體單元中共有4個(gè)四面體單元的體積為零。我們必須要把這四個(gè)四面體體積為零的四面體作上標(biāo)記，使其不參加后面的有限元方法的計(jì)算，以免影響整個(gè)程序計(jì)算的精度與速度。綜合上述考慮，那么對于心肌這種被認(rèn)為是近似不可壓縮的物質(zhì)，我們參考人體血流動(dòng)力學(xué)參數(shù)取其泊松比v都取0.4，彈性模量即楊氏模量E取50kPa，質(zhì)量密度ρ取1.07g/cm³。我們在有限元模型的心臟表面位移計(jì)算算法研究中，實(shí)際程序中采用的點(diǎn)是771個(gè)點(diǎn)減去這四個(gè)四面體體積為零的12個(gè)點(diǎn)后的759個(gè)節(jié)點(diǎn)，619個(gè)常應(yīng)變四面體單元來計(jì)算的。

求節(jié)點(diǎn)位移矢量的算法思路設(shè)計(jì)

在進(jìn)行有限元法求表面節(jié)點(diǎn)的位移矢量時(shí)，其中我們就要利用虛功原理與變分原理以及至關(guān)重要的邊界條件來進(jìn)行各初始條件的導(dǎo)入，進(jìn)而求得在載荷的邊界條件下心臟表面節(jié)點(diǎn)的位移矢量。在相關(guān)文獻(xiàn)中的論述中計(jì)算單位法線及高斯曲率時(shí)，采用等間隔分布的網(wǎng)格數(shù)據(jù)點(diǎn)（如與坐標(biāo)軸對應(yīng)的點(diǎn)）作為節(jié)點(diǎn)鄰域擬合曲面，也就是先將每層短軸方向的輪廓點(diǎn)擬合，得到短軸平面上形式統(tǒng)一的輪廓點(diǎn)數(shù)據(jù)；然后將長軸輪廓點(diǎn)數(shù)據(jù)轉(zhuǎn)換成短軸平面數(shù)據(jù)形式，得到形式完全一致的心臟內(nèi)輪廓表面線模型，考慮了肌纖維方向及其橫截面方向的差異，運(yùn)用了運(yùn)動(dòng)心臟模擬器驗(yàn)證本書模型的結(jié)果，簡言之就是等參元、等間隔分層長、短軸法。由于先前采用的四面體算法不是按分層等間隔來進(jìn)行考慮的，因此上述方法并不太適用于我們所建的心臟表面有限元模型，那么我們采用的是近似橢球體的各參數(shù)來進(jìn)行模擬。其思想是：α表示沿模型的基部到頂部方向的取值為-π/2≤α≤π/2，β表示沿模型圓周方向的取值，-π≤β≤π，其邊界載荷力F是基于系數(shù)函數(shù)cosβ，cosα，sinα的形式，向外、向內(nèi)的F分別表示舒張、收縮的載荷力均由橢球體的參數(shù)函數(shù)來確定?？紤]到進(jìn)行的是心臟表面的四面體模型，并由于心臟本身的功能特點(diǎn)及人體血流動(dòng)力學(xué)參數(shù)，心臟的收縮與舒張主要發(fā)生在下半部，即左、右心室，故應(yīng)著重使得分析中的中心坐標(biāo)點(diǎn)下移到心臟下半部的左、右心室的中心，使其滿足并接近真實(shí)的心臟收縮與舒張形態(tài)。