任 亮 邱天爽 郭 勇

(大連理工大學(xué)生物醫(yī)學(xué)工程系，大連 116024)

引言

超聲影像在醫(yī)學(xué)診斷和治療中的應(yīng)用被得到廣泛認(rèn)可［1］。與 X－ray、CT、MRI相比，超聲影像具有實(shí)時(shí)、多角度、無(wú)輻射損傷、廉價(jià)等優(yōu)勢(shì)，特別在穿刺活檢和射頻消融手術(shù)中得到廣泛應(yīng)用［2］。普遍使用的二維超聲影像存在部分容積效應(yīng)，使得病灶靶點(diǎn)在三維空間中的定位不準(zhǔn)確，影響診斷和治療效果;而三維超聲影像可以提供更精確和豐富的三維空間信息，有利于病灶靶點(diǎn)的定位和后續(xù)手術(shù)操作。目前，重建三維超聲影像的方法主要有3種［3］:二維超聲探頭機(jī)械掃描［4］、三維超聲探頭實(shí)時(shí)成像和手持二維超聲探頭掃描［5］。前兩種方法改裝后的探頭體積大，限制其在手術(shù)導(dǎo)航中的應(yīng)用，而手持掃描的方式僅需要進(jìn)行超聲探頭標(biāo)定。求解超聲圖像相對(duì)超聲探頭的空間變換矩陣的過(guò)程，稱(chēng)為超聲探頭標(biāo)定［6］。

目前，國(guó)內(nèi)外普遍采用基于平面N形模板的超聲探頭標(biāo)定方法［7－8］。該方法在電磁定位條件下，用探頭掃描N形標(biāo)定模板，得到超聲聲束與模板交點(diǎn)的超聲圖像坐標(biāo)和模板坐標(biāo);通過(guò)最小二乘法，分別標(biāo)定出電磁發(fā)射器相對(duì)標(biāo)定模板、超聲圖像相對(duì)電磁接收器的兩個(gè)未知空間變換矩陣。該方法能夠滿(mǎn)足手術(shù)導(dǎo)航的精度要求［9］，但是在電磁定位條件下，標(biāo)定算法包含兩個(gè)未知空間變換矩陣，共14個(gè)未知變量，平面模板空間信息少，需要采集大量的標(biāo)定圖像，計(jì)算量較大。另外，電磁定位設(shè)備易受手術(shù)導(dǎo)航環(huán)境中復(fù)雜的電磁場(chǎng)影響［10］，且電磁接收器需要連接電源線(xiàn)，這會(huì)限制探頭的活動(dòng)空間。

為此，筆者采用一種區(qū)別于文獻(xiàn)［7－8］的基于空間立體標(biāo)定模板的超聲探頭標(biāo)定方法，提出如下改進(jìn):一是在光學(xué)定位條件下，僅需要計(jì)算超聲圖像相對(duì)超聲探頭的一個(gè)未知空間變換矩陣，共8個(gè)未知變量，減少了未知量個(gè)數(shù);二是新設(shè)計(jì)一種由兩個(gè)互相垂直平面組成的空間立體標(biāo)定模板，增加聲束平面與模板交點(diǎn)的空間信息，且不需要已知模板的幾何尺寸，制作簡(jiǎn)單。實(shí)驗(yàn)表明，本方法僅需要兩幅標(biāo)定圖像即可完成標(biāo)定，計(jì)算量小。當(dāng)采用至少20幅標(biāo)定圖像時(shí)，標(biāo)定精度約為2.5 mm;當(dāng)采用至少40幅標(biāo)定圖像時(shí)，平均誤差約為1.0 mm，標(biāo)準(zhǔn)差約為2.4 mm。

1 方法

1.1 數(shù)學(xué)模型

超聲探頭標(biāo)定如圖1所示。筆者設(shè)定4個(gè)坐標(biāo)系，即標(biāo)定模板坐標(biāo)系M、攝像機(jī)坐標(biāo)系C、超聲探頭坐標(biāo)系D、超聲圖像坐標(biāo)系I，3個(gè)空間變換矩陣TCM、TDC和TID。其中TCM是坐標(biāo)系C到坐標(biāo)系M的變換矩陣，TDC是坐標(biāo)系D到坐標(biāo)系C的變換矩陣，TID是坐標(biāo)系I到坐標(biāo)系D的變換矩陣。

筆者提出標(biāo)定模板，由兩個(gè)互相垂直的平面OMABC和平面OMCDE組成。在采集超聲圖像時(shí)，探頭聲束平面與標(biāo)定模板的三坐標(biāo)軸分別交于點(diǎn)F、G、H。設(shè)定標(biāo)定模板坐標(biāo)系M、攝像機(jī)坐標(biāo)系C、超聲探頭坐標(biāo)系D、超聲圖像坐標(biāo)系I的坐標(biāo)軸，如圖1所示。光學(xué)定位設(shè)備由兩個(gè)型號(hào)相同的攝像機(jī)組成，利用文獻(xiàn)［11］的方法標(biāo)定攝像機(jī)。設(shè)標(biāo)定模板與探頭聲束交于點(diǎn)P，則標(biāo)定公式

式中，PM= ［XM，YM，ZM，1］T是點(diǎn) P 在坐標(biāo)系 M 下的空間坐標(biāo)，PI=［sx·u，sy·v，0，1］T是點(diǎn) P 在坐標(biāo)系I下的圖像坐標(biāo)，sx和sy是需要標(biāo)定的未知比例因子，變換矩陣TCM和TDC可以由光學(xué)定位設(shè)備獲得，TID為待求的變換矩陣，包含 tx、ty、tz、α、β、γ共6個(gè)未知空間變量。

在標(biāo)定過(guò)程中，坐標(biāo)系M和坐標(biāo)系C固定不動(dòng)，手持超聲探頭在攝像機(jī)視野內(nèi)移動(dòng)并掃描標(biāo)定模板。在超聲探頭采集超聲圖像的同時(shí)，采集探頭在攝像機(jī)中的圖像，求解超聲圖像坐標(biāo)系I與超聲探頭坐標(biāo)系D的變換矩陣TID的過(guò)程就是超聲探頭標(biāo)定過(guò)程。

1.2 標(biāo)定算法

如圖1所示，對(duì)于點(diǎn)F，根據(jù)式(1)有

按式(3)中左側(cè)的兩個(gè)零元素展開(kāi)，可以得到關(guān)于未知量 TID(tx，ty，tz，α，β，γ)、sx和 sy的兩個(gè)標(biāo)定方程。同理，對(duì)于Y軸和Z軸上的點(diǎn)G和點(diǎn)H，可以分別得到兩個(gè)標(biāo)定方程，所以每幅標(biāo)定圖像可以得到6個(gè)標(biāo)定方程。有8個(gè)未知量，當(dāng)標(biāo)定圖像數(shù)量N滿(mǎn)足6N≥8時(shí)，即采用至少2幅可以非線(xiàn)性?xún)?yōu)化求解。

圖1 超聲探頭標(biāo)定Fig.1 Illustration of ultrasound probe calibration

根據(jù)式(3)，對(duì)于超聲聲束與坐標(biāo)系M的X軸的交點(diǎn)，設(shè)

令 fX，2(i)=0 和fX，3(i)=0，分別對(duì)應(yīng)式(4)左側(cè)兩個(gè)零元素展開(kāi)得到標(biāo)定方程。同理，有fY，1(i)=0、fY，3(i)=0 和 fZ，1(i)=0、fZ，2(i)=0。構(gòu)造目標(biāo)函數(shù)

令 ψ = ［tx，ty，tz，α，β，γ，sx，sy］，當(dāng) N ≥2 時(shí)，堆疊標(biāo)定方程組 f(ψ)= ［f1(ψ);f2(ψ);…;fN(ψ)］，則

根據(jù) Levenberg－Marquardt算法，有

式中，η是設(shè)定的迭代收斂因子，I是單位矩陣。每次迭代后，更新J和Δf，直到Δψ趨近于0時(shí)，ψ即為筆者所求的標(biāo)定參數(shù)。

2 實(shí)驗(yàn)

2.1 仿真數(shù)據(jù)實(shí)驗(yàn)

按照式(2)，利用 Matlab隨機(jī)生成1組 TCM，100組TDC(i)、100組包含超聲聲束與標(biāo)定模板3個(gè)交點(diǎn)的超聲圖像坐標(biāo)PI(i)。根據(jù)式(6)，令每組實(shí)驗(yàn)3個(gè)交點(diǎn)F、G、H模板坐標(biāo)的零元素估計(jì)值為)=［fF，2(i);fF，3(i);fG，1(i);fG，3(i);fH，1(i);fH，2(i)］，理論值為fi(ψ)=[0，0，0，0，0，0]，標(biāo)定誤差為每組有3個(gè)交點(diǎn)，統(tǒng)計(jì)N組實(shí)驗(yàn)的標(biāo)定平均誤差為其中，2≤N≤100，評(píng)價(jià)式(8)的標(biāo)定平均誤差Emean和標(biāo)準(zhǔn)差隨標(biāo)定圖像數(shù)量N增加的變化趨勢(shì)，每組實(shí)驗(yàn)重復(fù)100次。

2.2 真實(shí)數(shù)據(jù)實(shí)驗(yàn)

2.2.1 實(shí)驗(yàn)材料

實(shí)驗(yàn)平臺(tái)如圖2所示，標(biāo)定模板由兩塊互相垂直的平整光滑的玻璃板A和B組成，長(zhǎng)寬高分別為250 mm×150 mm×200 mm，板厚為3 mm。在實(shí)驗(yàn)中，一臺(tái)PC機(jī)通過(guò)圖像采集卡與光學(xué)定位設(shè)備和超聲儀器連接，利用圖像采集軟件，保證在采集攝像機(jī)圖像的同時(shí)觸發(fā)采集超聲圖像，光學(xué)定位設(shè)備由兩個(gè)型號(hào)相同的攝像機(jī)組成，圖像分辨率均為1024像素×768像素，利用文獻(xiàn)［11］標(biāo)定攝像機(jī)后空間點(diǎn)重建精度可以達(dá)到1 mm。在模板表面，粘貼黑白棋盤(pán)格特征標(biāo)識(shí)物，識(shí)別棋盤(pán)格角點(diǎn)并重建角點(diǎn)的三維坐標(biāo)，根據(jù)解析幾何原理，計(jì)算攝像機(jī)坐標(biāo)系C到標(biāo)定模板坐標(biāo)系M的變換矩陣TCM;為保證超聲在水槽介質(zhì)中的傳播速度與在人體中的傳播速度一致，向水槽注入適量的50℃熱水［7］。采用Aloka Prosund α7超聲儀的腹部探頭，在探頭上粘貼黑白棋盤(pán)格特征標(biāo)識(shí)物。

2.2.2 實(shí)驗(yàn)過(guò)程

實(shí)驗(yàn)主要包括4個(gè)步驟。

步驟1:將標(biāo)定模板放入空水槽中固定，利用光學(xué)定位設(shè)備計(jì)算攝像機(jī)坐標(biāo)系C到標(biāo)定模板坐標(biāo)系M的變換矩陣TCM;

步驟2:向水槽中注入適量的50℃熱水，手持超聲探頭在標(biāo)定模板上方以不同的姿態(tài)掃描標(biāo)定模板，利用光學(xué)定位設(shè)備識(shí)別探頭上的黑白棋盤(pán)格角點(diǎn)，并重建角點(diǎn)的三維坐標(biāo)，繼而計(jì)算超聲探頭坐標(biāo)系D到攝像機(jī)坐標(biāo)系C的變換矩陣TDC(i)，通過(guò)軟件觸發(fā)方式在同一時(shí)刻采集超聲圖像，共采集100幅標(biāo)定圖像。

步驟3:手動(dòng)選取識(shí)別超聲聲束與標(biāo)定模板的3個(gè)交點(diǎn)F、G、H的圖像坐標(biāo)，如圖3所示。

步驟 4:利用本方法計(jì)算 TID(tx，ty，tz，α，β，γ)、sx和sy，并評(píng)價(jià)式(8)的標(biāo)定平均誤差Emean和標(biāo)準(zhǔn)差隨標(biāo)定圖像數(shù)量N增加的變化趨勢(shì)。

3 結(jié)果

3.1 仿真實(shí)驗(yàn)結(jié)果

圖2 實(shí)驗(yàn)平臺(tái)搭建。(a)標(biāo)定實(shí)物圖;(b)標(biāo)定模板Fig.2 Calibration setup for real experiments.(a)Calibration equipment;(b)Calibration pattern

圖3 交點(diǎn)F、G、H的超聲圖像Fig.3 F、G、H in ultrasound image

標(biāo)定圖像數(shù)量由2～100幅，每增加5幅圖像進(jìn)行1組實(shí)驗(yàn)，評(píng)價(jià)式(8)的標(biāo)定平均誤差Emean和標(biāo)準(zhǔn)差隨標(biāo)定圖像數(shù)量N增加的變化趨勢(shì)，每組實(shí)驗(yàn)重復(fù)100次，標(biāo)定結(jié)果如圖4所示。隨著標(biāo)定圖像數(shù)量的增加，標(biāo)定平均誤差和標(biāo)準(zhǔn)差均逐漸減小，且在僅有兩幅標(biāo)定圖像的情況下，標(biāo)定平均誤差小于5 mm，當(dāng)圖像數(shù)量多于40幅時(shí)，標(biāo)定平均誤差小于1 mm，標(biāo)準(zhǔn)差約為2 mm，說(shuō)明本算法具有較好的穩(wěn)定性。

圖4 標(biāo)定算法的仿真實(shí)驗(yàn)結(jié)果Fig.4 Calibration results with simulation data

3.2 真實(shí)實(shí)驗(yàn)結(jié)果

標(biāo)定圖像數(shù)量由2～100幅，每增加5幅圖像進(jìn)行1組實(shí)驗(yàn)，評(píng)價(jià)式(8)的標(biāo)定平均誤差Emean隨著標(biāo)定圖像數(shù)量N增加的變化趨勢(shì)，每組實(shí)驗(yàn)重復(fù)100次，標(biāo)定結(jié)果如圖5所示。采用至少20幅標(biāo)定圖像時(shí)的標(biāo)定平均誤差約為2.5 mm，采用至少40幅標(biāo)定圖像時(shí)的平均誤差約為1.0 mm，標(biāo)準(zhǔn)差約為2.4 mm。

本標(biāo)定方法與文獻(xiàn)［7－8］中各種方法的比較如表1所示。與基于電磁定位的方法相比，本方法在達(dá)到與Cambridge方法標(biāo)定精度一致的情況下，所用的標(biāo)定圖像數(shù)量為40幅，且在標(biāo)定精度方面優(yōu)于文獻(xiàn)［7］的方法和 Cross－wire、Single－wall、Threewire等3種經(jīng)典方法，說(shuō)明本算法具有較好的有效性。

圖5 標(biāo)定算法的真實(shí)實(shí)驗(yàn)結(jié)果Fig.5 Calibration results with real data

由此可見(jiàn)，本算法在光學(xué)定位條件下僅有8個(gè)未知量，所以只需要較少的標(biāo)定圖像即可完成標(biāo)定;且本設(shè)計(jì)的立體標(biāo)定模板增加了空間信息，有利于提高標(biāo)定精度。本標(biāo)定誤差大于文獻(xiàn)［1］的標(biāo)定誤差約0.1 mm，原因在于本標(biāo)定模板的制作過(guò)程簡(jiǎn)單，為保證兩平面的垂直和剛性，平面需要有一定厚度，這在超聲圖像中會(huì)造成一定的噪聲，使手動(dòng)識(shí)別超聲圖像中的特征點(diǎn)時(shí)產(chǎn)生誤差，影響標(biāo)定精度。

另外，隨著標(biāo)定圖像數(shù)量的增加，本真實(shí)實(shí)驗(yàn)誤差的減小速度快于仿真實(shí)驗(yàn)誤差的減小速度，原因在于真實(shí)實(shí)驗(yàn)中的模板制作、環(huán)境光線(xiàn)等對(duì)圖像識(shí)別引入不可避免的隨機(jī)誤差，圖像數(shù)量少時(shí)隨機(jī)誤差影響較大，隨著標(biāo)定圖像數(shù)量的增加，每組實(shí)驗(yàn)通過(guò)100次平均可以快速地減小隨機(jī)誤差，而仿真實(shí)驗(yàn)中不存在隨機(jī)噪聲，誤差的減小趨勢(shì)較為平緩，因此真實(shí)實(shí)驗(yàn)誤差的減小速度更快。

表1 真實(shí)圖像的標(biāo)定結(jié)果Tab.1 Calibration results with real data

4 討論

與文獻(xiàn)［7－8］相比，筆者利用光學(xué)定位技術(shù)獲取超聲探頭在空間中的姿態(tài)，需要求解的未知量個(gè)數(shù)減少為8個(gè)，進(jìn)而降低了計(jì)算量。同時(shí)，在文獻(xiàn)［8］的 Cross－wire、Single－wall和 Three－wire 等 3 種方法的基礎(chǔ)上，將標(biāo)定模板設(shè)計(jì)為由兩個(gè)互相垂直的空間平面組成，增加了標(biāo)定信息，僅需采集至少兩幅圖像即可完成標(biāo)定。這些都說(shuō)明，在光學(xué)定位下，通過(guò)增加標(biāo)定模板立體信息，可以減少標(biāo)定所需的圖像數(shù)量和計(jì)算量。

由于本研究只是識(shí)別超聲聲束與標(biāo)定模板在圖像中的交點(diǎn)，在人眼的視野中無(wú)法識(shí)別出該交點(diǎn)的真實(shí)位置，因此無(wú)法評(píng)價(jià)該交點(diǎn)測(cè)量的精度，本研究通過(guò)對(duì)目標(biāo)函數(shù)理論值的逼近程度來(lái)評(píng)價(jià)標(biāo)定精度。這種評(píng)價(jià)方式與文獻(xiàn)［1，7－8］中的Crosswire、Single－wall和 Three－wire 等方法一致。由表 1看出，本研究的仿真和真實(shí)實(shí)驗(yàn)結(jié)果在圖像數(shù)量多于40幅時(shí)，平均標(biāo)定誤差約為1 mm，優(yōu)于文獻(xiàn)［7－8］的精度，滿(mǎn)足超聲影響導(dǎo)航手術(shù)的需求，并且在光學(xué)定位下不需要考慮手術(shù)室內(nèi)電源線(xiàn)和電磁兼容的問(wèn)題，使可標(biāo)定的范圍更廣。

另外，可以考慮將本標(biāo)定模板退化為空間單獨(dú)的一個(gè)點(diǎn)，這樣可以利用光學(xué)定位識(shí)別該點(diǎn)空間坐標(biāo)作為真值，繼而評(píng)價(jià)標(biāo)定精度，但是需要采集較多數(shù)量的標(biāo)定圖像，這些是今后繼續(xù)研究的內(nèi)容。

5 結(jié)論

針對(duì)超聲探頭標(biāo)定問(wèn)題，筆者提出一種在光學(xué)定位條件下基于立體標(biāo)定模板的超聲探頭標(biāo)定方法，既減少了未知量個(gè)數(shù)，又增加了標(biāo)定空間信息。在標(biāo)定過(guò)程中，利用光學(xué)定位設(shè)備獲取模板坐標(biāo)系和探頭坐標(biāo)系的空間姿態(tài)，軟件觸發(fā)的同時(shí)采集超聲圖像和探頭在攝像機(jī)中的圖像，兩正交平面組成的立體標(biāo)定模板利用3組正交關(guān)系堆疊標(biāo)定方程。本方法僅需要2幅以上標(biāo)定圖像即可優(yōu)化求解出超聲圖像到探頭的空間變換矩陣，計(jì)算量小，操作簡(jiǎn)便，且所求未知量?jī)H有8個(gè)，少于同類(lèi)文獻(xiàn)的未知量數(shù)量。真實(shí)實(shí)驗(yàn)結(jié)果表明，采用至少40幅圖像時(shí)的標(biāo)定平均誤差約為1.01 mm，標(biāo)準(zhǔn)差約為2.4 mm，在標(biāo)定圖像數(shù)量較少的情況下具有較高的精度，為進(jìn)一步三維超聲重建奠定了基礎(chǔ)。

［1］ChaouiJ，Dardenne G，Hamitouche C，etal. Virtual movements－based calibration method of ultrasound probe for computer assisted surgery［C］//IEEE International Symposium on Biomedical Imaging:From Nano to Macro.Bostoa:IEEE，2009:1207－1210.

［2］徐靜，楊向東，朱森強(qiáng)，等.用于肝癌介入治療的術(shù)中三維超聲導(dǎo)航系統(tǒng)［J］.中國(guó)生物醫(yī)學(xué)工程學(xué)報(bào)，2007，26(5):719－723.

［3］白培瑞，白凈.基于二維陣列換能器的實(shí)時(shí)三維超聲成像技術(shù)綜述［J］.中國(guó)生物醫(yī)學(xué)工程學(xué)報(bào)，2011，30(6):938 －945.

［4］Nakamoto M，Nakada K，Sato Y，et al.Intraoperative magnetic tracker calibration using a magneto－optic hybrid tracker for 3 － D ultrasound－based navigation in laparoscopic surgery［J］.IEEE Transactions on Medical Imaging，2008，27(2):255 －270.

［5］Marie－Ange J，Louis－GillesD，Marie－HélèneRC，etal.Performance evaluation ofa medicalrobotic 3D－ultrasound imaging system［J］.Medical Image Analysis，2008，12(3):275－290.

［6］Boctor E，Viswanathan A，Choti M，et al.A novel closed form solution for ultrasound calibration［C］//IEEE International Symposium on Biomedical Imaging:Nano to Macro.Arlington:IEEE，2004:527－ 530.

［7］羅楊宇，徐靜，魯通，等.基于磁定位器的手動(dòng)三維超聲圖像標(biāo)定［J］.中國(guó)生物醫(yī)學(xué)工程學(xué)報(bào)，2008，27(2):250－254.

［8］Prager RW，Rohling RN，Gee AH，et al.Rapid calibration for 3－D freehand ultrasound［J］.Ultrasound in Medicine ＆ Biology，1998，24(6):855 －869.

［9］Birkfellner W，Watzinger F，Wanschitz F，et al.Systematic distortions in magnetic position digitizers［J］.Medical Physics，1998，25(11):2242 －2248.

［10］Mercier L，Lang? T，Lindseth F，et al.A review of calibration techniques for freehand 3 － D ultrasound systems［J］.Ultrasound in Medicine＆ Biology，2005，31(4):449－471.

［11］Zhang Zhengyou.A flexible new technique for camera calibration［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，2000，22(11):1330 －1334.