本發(fā)明涉及虛擬手術軟組織力反饋模型建模，用于實現(xiàn)在虛擬手術中生物軟組織力反饋模型與操作者的精確性實時性交互。
背景技術：
：一個理想的生物軟組織模型能夠準確的將軟組織位移情況及受力大小即時反饋給操作者，這就要求生物軟組織模型具有高精確性及實時性。而在虛擬手術中，高效的軟組織力反饋物理模型建立對虛擬手術系統(tǒng)的成功起著至關重要的作用。目前，基于物理驅動的軟組織模型日漸成熟，大致可分為兩類：一種是基于網絡結構，主要有質點-彈簧模型和有限元模型。質點-彈簧模型(MS)的優(yōu)點是原理簡單、計算量比較小和實時性好，缺點是形變具有局部性質，并且對模型中彈簧剛度系數(shù)設置合理的材料參數(shù)是沒有理論依據(jù)，很容易造成軟組織在形變仿真的時候，出現(xiàn)失真問題。有限元模型(FEM)的優(yōu)點是計算精度高，但是由于求解過程中需要組裝剛度矩陣，計算量比較大，如果進行大形變或者切割的時候，網格模型的拓撲結構發(fā)生復雜的變化，使得形變實時性差；另一種是無網格結構，主要方法有粒子系統(tǒng)、有限球方法和無網格方法。HieberSE等基于線彈性力學模型，利用光滑粒子流(SPH)無網格法，不需要網格拓撲結構，實時性高，可以滿足軟組織大范圍的形變，但存在穩(wěn)定性差以及不具備粘彈性等問題，該方法常用于快速粒子流或氣體等情況的仿真。比如，流血、爆炸和煙霧等。SuzukiS的填充球模型，是一種基于幾何驅動的無網格法，不包含軟組織的生物特性，不能滿足虛擬手術仿真的高精確度要求。無網格法是機械工程領域中的數(shù)值方法，近年來發(fā)展迅速，出現(xiàn)了如無網格伽遼金法(EFG)、徑向基函數(shù)法(RBF)和多尺度重構核粒子方法(MRKP)等多種無網格法。它采用一組離散的場節(jié)點來構造問題域和邊界，通過插值來生成每個場節(jié)點的形函數(shù)，在不需要任何點與點之間拓撲關系的情況下近似逼近問題域內未知的場函數(shù)，最后求解離散系統(tǒng)方程，得到每個場節(jié)點的位移、應力和應變。無網格法常用于求解懸臂梁等機械固體受力后的形變，在軟組織形變方面并沒有形成一套成熟的理論系統(tǒng)。由于無網格法本身存在的特性使得其可以應用與軟組織形變仿真中。其優(yōu)點：第一，自適應性強，采用點云結構模型，不需要點與點之間復雜的拓撲結構，可適用于大形變或者切割。第二，建立模型只需要節(jié)點的信息。第三，計算結果光滑連續(xù)，適合力反饋設備。但是由于無網格存在計算效率低下的問題，所以在實時性方面并不是很理想。綜上所述可知，基于網絡結構的質點-彈簧模型(MS)結構簡單實時性好，但其在形變上具有局部性質，精確性不足。有限元計算精度高但計算量大形變實時性差。而無網格結構具有非常高的精確性及自適應性滿足虛擬手術中反饋精確要求，但是無網格法依然面臨計算效率低下，在實時性方面有所缺陷。技術實現(xiàn)要素：本發(fā)明的目的是提出一種基于麥夸特算法的徑向基無網格軟組織數(shù)據(jù)的力反饋模型建模方法。針對現(xiàn)有生物軟組織物理模型缺陷，考慮生物軟組織各向異性特性，由徑向基無網格法一次性獲得以軟組織受力點為中心，面上有限個等間距點受力位移數(shù)據(jù)，使用麥夸特法擬合軟組織節(jié)點力-位移模型，再以點至線，以線帶面，建立生物軟組織力反饋模型。在保留徑向基無網格法精確性的同時，又大大簡化計算量，保證了交互的高實時性，實現(xiàn)虛擬手術中生物軟組織力反饋模型與操作者的精確性實時性交互。本發(fā)明是通過以下技術方案實現(xiàn)的。生物軟組織第hαβ節(jié)點有如下系統(tǒng)：以受力點O為坐標原點，其中x、y、z分別為hαβ節(jié)點橫軸方向位移量、縱軸方向位移量及垂直水平面方向位移量，F(xiàn)為受力點O反饋力大小，為隨機誤差項。利用徑向基無網格法獲取以受力點O為中心，間距為R的H個矩形點陣節(jié)點的軟組織輸出力F(k)和輸入位移x(k),y(k),z(k)數(shù)據(jù)。(在現(xiàn)實手術中生物組織受力形變的組織直徑范圍在1～4cm左右)。α×β＝H考慮生物軟組織力與形變的特點可將系統(tǒng)寫成：其中x、y、z分別為hαβ節(jié)點橫軸方向位移量、縱軸方向位移量及垂直水平面方向位移量，p為系統(tǒng)待估參數(shù)。開始取α＝1，β＝1。步驟(1)由徑向基無網格法獲得各節(jié)點軟組織力與形變輸入輸出數(shù)據(jù)，采樣當前第hαβ節(jié)點輸入輸出數(shù)據(jù)。獲取L組軟組織數(shù)據(jù)。步驟(2)為避免過度擬合，利用交叉驗證法先確定hαβ節(jié)點系統(tǒng)自變量階數(shù)na、nb、nc，而確定了自變量階數(shù)的同時也確定了參數(shù)p的個數(shù)m。步驟(3)設置待估參數(shù)p初值設定阻尼因子d＝10ud(0)(d(0)>0)。步驟(4)設定迭代次數(shù)u＝0，進行第一次迭代，利用L組軟組織數(shù)據(jù)、參數(shù)初值p(0)及阻尼因子初值d(0)通過麥夸特算法參數(shù)估計解得p值，將求得參數(shù)p代入系統(tǒng)原函數(shù)式，計算殘差平方和J(0)。步驟(5)將步驟(4)計算獲得參數(shù)p賦予p(0)，即p(0)＝p，令u＝1即d＝10d(0)進行第二次迭代，再執(zhí)行步驟(4)解得新的參數(shù)p并代入系統(tǒng)原函數(shù)式，計算新的殘差平方和J(1)。比較J(1)與J(0)，增加迭代次數(shù)，直至J(1)<J(0)為止。步驟(6)繼續(xù)增加阻尼因子d的數(shù)量級即u＝u+1，反復執(zhí)行步驟(5)，直至新計算獲得參數(shù)p與上一次計算參數(shù)即p(0)之差絕對值|Φ|≤δ允許最小誤差為止。此時hαβ節(jié)點系統(tǒng)參數(shù)估計完成。步驟(7)令α＝α+1回到步驟(1)重新執(zhí)行，即對節(jié)點矩陣縱向建模，當時，步驟(7)結束。步驟(8)令α＝1，β＝β+1回到步驟(1)重新執(zhí)行，即對節(jié)點矩陣橫向建模，當時，步驟(8)結束。至此建立了H個矩陣節(jié)點輸入輸出系統(tǒng)，即建立以受力點O為中心面積為HR2生物軟組織力反饋矩陣模型。本發(fā)明的優(yōu)點：以麥夸特法快速擬合徑向基無網格法獲得的軟組織輸入輸出數(shù)據(jù)，建立生物軟組織力反饋模型，該方法結構簡潔、計算方便，在保留徑向基無網格法精確性的同時，又大大簡化計算量，保證了交互的高實時性，實現(xiàn)虛擬手術系統(tǒng)軟組織力反饋模型高效建模。附圖說明圖1為發(fā)明流程框圖。圖2為力反饋示意圖。具體實施方式本發(fā)明將通過以下實施例作進一步說明。生物軟組織第hαβ節(jié)點有如下系統(tǒng)：以受力點O為坐標原點，其中x、y、z分別為hαβ節(jié)點橫軸方向位移量、縱軸方向位移量及垂直水平面方向位移量，F(xiàn)為受力點O反饋力大小，為隨機誤差項。利用徑向基無網格法獲取以受力點O為中心，間距為R的H個矩形點陣節(jié)點的軟組織輸出力F(k)和輸入位移x(k),y(k),z(k)數(shù)據(jù)。(在現(xiàn)實手術中生物組織受力形變的組織直徑范圍在1～4cm左右)。α×β＝H考慮生物軟組織力與形變的特點可將系統(tǒng)寫成：令原系統(tǒng)函數(shù)式可以寫成其中x、y、z分別為hαβ節(jié)點橫軸方向位移量、縱軸方向位移量及垂直水平面方向位移量，p為系統(tǒng)待估參數(shù)。通過徑向基無網格算法提供第hαβ節(jié)點的M組軟組織力與形變數(shù)據(jù)：(xi1,xi2,...,xina;yi1,yi2,...,yinb;zi1,zi2,...,zinc;Fi),i=1,2,...,M]]>開始取α＝1，β＝1。步驟1采樣當前由徑向基無網格法獲得的生物軟組織第hαβ節(jié)點輸入輸出數(shù)據(jù)。獲取L組軟組織數(shù)據(jù)：(xi1,xi2,...,xina;yi1,yi2,...,yinb;zi1,zi2,...,zinc;Fi),i=1,2,...,L]]>步驟2為避免過度擬合，先確定hαβ節(jié)點系統(tǒng)自變量階數(shù)na、nb，而確定了自變量階數(shù)的同時也確定了參數(shù)p的個數(shù)。確定自變量階數(shù)具體過程描述如下：(a)利用交叉驗證法，有L個樣本，對于na階水平方向位移自變量x，設縱軸位移、垂直水平位移自變量y，z為常量C、D，對x做L次擬合，將第i(i＝1,2,…L)次觀測值剔除，計算殘差平方和,尋找到使其最小的階數(shù)A即其中b-i,a為刪掉第i個觀測之后的估計值。(b)設定橫向位移x、垂直水平位移z為常量C、D，同理可求得nb。(b)設定水平位移自變量x、y為常量C、D，同理可求得nc。(d)初步描述系統(tǒng)為：F=p1+p2x1+p3x2+...+pna+1xna+pna+2y1+pna+3y2+...+pna+nb+1ynb+pna+nb+nc+1znc---(1)]]>步驟3設置待估參數(shù)p初值設定阻尼因子d＝10ud(0)(d(0)>0)；設定迭代次數(shù)u＝0即d＝d(0)，對于阻尼因子d＝0即高斯-牛頓法，在高斯-牛頓法中參數(shù)初值選擇的合適與否,將決定迭代計算的工作量,甚至迭代成功與否。即使有較好的選擇初值的方法,也難免出現(xiàn)發(fā)散的情況，而麥夸爾特法由于阻尼因子的引入，在一般的初值選擇條件下都可收斂于所求結果。步驟4設定迭代次數(shù)u＝0，進行第一次迭代，利用L組軟組織數(shù)據(jù)、參數(shù)初值p(0)及阻尼因子初值d(0)通過麥夸特算法參數(shù)估計解得p值，將求得參數(shù)p代入系統(tǒng)原函數(shù)式，計算殘差平方和J(0)。麥夸特求解待估參數(shù)具體過程描述如下：(a)將f(xi,yi,zi,p)在p(0)處按泰勒級數(shù)展開并略去二次及二次以上的項得：f(xi,yi,zi,p)≈f(xi,yi,zi,p)+∂f(xi,yi,zi,p)∂p1|p=p(0)(p1-p1(0))+∂f(xi,yi,zi,p)∂p2|p=p(0)(p2-p2(0))+...∂f(xi,yi,zi,p)∂pm|p=p(0)(pm-pm(0))---(2)]]>(b)式(2)中除待估參數(shù)(p1,p2,…pm)外，都是已知。對此使用最小二乘殘差平方原理：J=ΣiL{Fi-[f(xi,yi,zi,p(0))+Σj=1m∂f(xi,yi,zi,p(0))∂pj|p=p(0)(pj-pj(0))]}2+dΣjm(pj-pj(0))2---(3)]]>其中d≥0為阻尼因子，當d＝0時就是牛頓-高斯法，他對迭代初值的選擇要比麥夸特法更加嚴格。(c)通過將殘差平方J分別對p1,p2,…pm一階偏導設為0，以達到將其最小化的目的：∂J∂pk=2Σi=1L[Fi-f(xi,yi,zi,p(0))+Σj=1m∂f(xi,yi,zi,p(0))∂pj|p=p(0)(pj-pj(0))]·∂f(xi,yi,zi,p(0))∂pk|p=p(0)+2d(pk-pk(0))=0,k=(1,2,...,m)---(4)]]>(d)令ajk=∂f(xi,yi,zi,p)∂pk|p=p(0)·Σi=1L∂f(xi,yi,zi,p)∂pj|p=p(0)ajy=∂f(xi,yi,zi,p)∂pj|p=p(0)·Σi=1L(f(xi,yi,zi,p)-Fi)j=1,2,...,m;k=1,2,...,m---(5)]]>可將式(4)寫成：(a11+d)(p1-p1(0))+a12(p2-p2(0))+...+a1m(pm-pm(0))=α1ya21(p1-p1(0))+(a22+d)(p2-p2(0))+...+a2m(pm-pm(0))=α2y...am1(p1-p1(0))+am2(p2-p2(0))+...+(amm+d)(pm-pm(0))=amy---(6)]]>(f)解得P=p1p2···pm=p1(0)p2(0)···pm(0)+a11+da12...a1ma21a22+d...a2m···am1am2...amm+d---(7)]]>步驟5將步驟4計算獲得參數(shù)p賦予p(0)，即p(0)＝p，令u＝1即d＝10d(0)進行第二次迭代，再執(zhí)行步驟4解得新的參數(shù)p并代入系統(tǒng)原函數(shù)式，計算新的殘差平方和J(1)。比較J(1)與J(0)，增加迭代次數(shù)，直至J(1)<J(0)為止。具體過程描述如下：比較：若J(1)<J(0)，步驟5結束；若J(1)≥J(0)，取u＝2即d＝100d(0)，再執(zhí)行步驟4重解p并計算新殘差平方J(1)，若新計算的J(1)<J(0)步驟5結束，反之則再取u＝3即d＝1000d(0)再執(zhí)行步驟4…，如此反復迭代直至J(1)<J(0)為止。在式(6)中，aij為定值，所以當d越大時Φ＝|p-p(0)|越小，殘差平方和J收斂，而d相對初值d(0)增加的越大即u越大迭代的次數(shù)越多。選擇的界限是看殘差平方和是否下降。步驟6繼續(xù)增加阻尼因子d的數(shù)量級即u＝u+1，反復執(zhí)行步驟5，直至新計算獲得參數(shù)p與上一次計算參數(shù)即p(0)之差絕對值|Φ|≤δ允許最小誤差為止，使得估計參數(shù)波動范圍忽略不計。此時hαβ節(jié)點系統(tǒng)參數(shù)估計完成。步驟7令α＝α+1回到步驟1重新執(zhí)行，即對節(jié)點矩陣縱向建模，當時，步驟7結束。步驟8令α＝1，β＝β+1回到步驟1重新執(zhí)行，即對節(jié)點矩陣橫向建模，當時，步驟8結束。至此建立了H個矩陣節(jié)點輸入輸出系統(tǒng)，即建立以受力點O為中心面積為HR2生物軟組織力反饋矩陣模型。當前第1頁1 2 3