計測自動制御学会東北支部 40 周年記念学術講演会
2004 年 12 月 22 日，23 日 日本大学工学部

講演番号：A08

血流の超音波計測融合シミュレーションにおけるエイリアジングの除去
○船本健一

早瀬敏幸

西條芳文

山家智之

（東北大学）

Elimination of Aliasing in Ultrasonic-Measurement-Integrated Simulation of Blood Flows*K. Funamoto, T. Hayase, Y. Saijo and T. Yambe (Tohoku University)
Abstract－Ultrasonic-Measurement-Integrated (UMI) Simulation, in which feedback signals from the ultrasonicmeasurement make the simulation to converge to the real blood flows, is a key to provide detailed information ofreal blood flows for development of an accurate diagnosis or treatment for serious circulatory diseases. However,aliasing in the ultrasonic blood velocity measurement causes UMI simulation to converge to an erroneous result. In
this paper, we have investigated the detection and the correction of aliasing in UMI simulation. As a result, it isrevealed that the artificial force in the feedback of UMI simulation can be used as an index to detect and correct thealiasing in original ultrasonic measurement.
Key Words: Blood flow, Measurement-integrated simulation, Ultrasonic measurement, Aliasing, Aneurysm

1

緒言

バック則の影響を明らかにした[8,9]．また，実際の超音波計測データとの融合を2次元数値シミュレーションで行い，計算結果を計測データに近づけることが可能であることも示した[8]．しかし，超音波計測融合シミュレーションは，計測誤差の一つであるエイリアジングも再現してしまい，その領域においてフィードバック量は異常に大きな値を示した．そこで，本研究では，2次元の数値シミュレーションにおいてフィードバック量を指標としたエイリアジングの検知を行い，2通りの方法によりその補正を行った結果について述べる．

近年，循環器系疾患と血行力学の関連性が指摘されている[1-4]．そのため，動脈瘤のような重篤な循環器系疾患の正確な診療を確立するためには，血管内の圧力分布や壁せん断応力などの詳細な情報が不可欠になってくる．
現在，医療現場における循環器系疾患の診断では，超音波診断装置が最も広く用いられている．超音波を用いる方法では，カラードプラ法により，生理的状態における血管形状と血流動態を非侵襲かつリアルタイムに表示することができる[5]．しかし，これは，血液の流速の計測にドプラ効果を利用するため， 2 計算方法
3次元性の強い血流に対し，ドプラ速度（速度の超音波ビーム方向の成分）しか計測することができない． 2.1 超音波計測融合シミュレーション
一方，血管内の血流の詳細な情報を提供する方法
超音波計測融合シミュレーションはFig. 1に示す
として，MRI，CTおよび超音波診断装置により得らように，超音波計測，数値シミュレーションおよびれた血管の実形状を用いる数値シミュレーションが， フィードバック則の3つの構成要素から成る．本解析
近年，盛んに行われている[6,7]．しかし，数値シミ手法では，通常の数値シミュレーションと同様に，ュレーションは，モデリングの精度や，有限領域の推定した境界条件を用いるが，超音波計測により得計算においては各瞬間の境界条件の設定に本質的なた血流のドプラ速度Vsと，数値シミュレーションに問題を抱えている．そのため，一般的に，実際の血より得たドプラ速度Vcを比較し，その誤差に基づく流と良く似た血流構造が得られるものの，完全に一フィードバック信号を加えながら計算を行う[8]．超致するものは得られない．
音波計測はカラードプラ法を用いて行う．Figure 2に，これまでに著者らは，計測と計算の両者がそれぞ
本研究対象とした，下行大動脈に発症した大動脈瘤れ有する問題点を解決するために，超音波による血の，経食道超音波走査により得られたカラードプラ
流計測と，コンピュータによる数値シミュレーショ画像を示す．Figure 2では，超音波プローブの位置（右ンを一体化した，超音波計測融合シミュレーション図の扇形の頂点）に近づいてくる血流に暖色系，遠[Ultrasonic-Measurement-Integrated (UMI) simulation]
ざかる血流に寒色系を用い，血管形状を表すBモードを開発してきた[8]．本手法は，超音波計測と数値シ画像に重ねて血流を表現してある．このカラードプ
ミュレーションの間に生じる誤差を，計算過程においてフィードバックし，計算結果を計測データに近づける．これまでに，その有用性の検証を2次元数値シミュレーションにより行い，境界条件やフィード

Fig.2 Schematic diagram of transesophageal
echography and obtained color Doppler image aroundan aortic aneurysm (Center frequency: 4.4MHz, pulse
repetition frequency: 4kHz)

Fig. 1 Block diagram of UMI simulation73ラ画像の輝度情報からドプラ速度Vsを得る．
Ultrasonic beam direction
本研究の数値シミュレーションに用いる2次元血j+1
管モデルは，カラードプラ画像中のBモード画像を2us値化し，そのデータを65×40の一様格子に割り当てる
Vs
Control volume
ことにより作成した（Fig. 3参照）．計算に用いるパVc
ラメータは，超音波計測時に同時に測定したデータ
ucuj等を基に，Table 1のように設定した．血液は密度ρ =33-31.0×10 [kg/m ]，粘度µ = 4.0 ×10 [Pa·s]のニュートンv
流体であるとした．また，上流端に与える血流量の時間変化は，Olufsenらの論文中のMRによる計測デー
タを参考にモデル化したものを用いた[10]．基礎方程j-1
式には2次元非圧縮性Navier-Stokes方程式と連続式をi+1
i-1用いる．以下の計算では，上流端の管形状を円形とi仮定し，その直径Dを代表長さとし，上行大動脈におFig. 4 Computational grid and definition of calculation
ける最大レイノルズ数の文献値Remaxと血液の動粘度
of feedback signal
νを用いて諸量を無次元化した．
らをSIMLER法に類似の手法により解いた[11-12]．
本研究では，超音波計測からは流れ場の限られた超音波計測融合シミュレーションにおけるフィー情報しか得られないことを考慮し，予め行った数値
ドバック則では，Fig.
3中に白い四角で示す，フィーシミュレーションの結果を基準解と定義する．ここ
ドバック領域内部の全格子点をフィードバックの代では，汎用のアプリケーション（FLUENT 6.1.22,表点として設定し，ドプラ速度Vに関して計測データ
Fluent Inc., Lebanon, NH, USA）により上行大動脈から（基準解）との誤差を求める．そして，それに比例腹部大動脈の間（分岐は無視）の血流の定常計算をした超音波ビーム方向の外力に相当するフィードバ行い，それにより得られた下行大動脈内の速度分布をFig. 3の血管形状の上流境界に与えて計算を行い， ック信号fvを，Navier-Stokes方程式に加えながら計算を行う．また，fvの作用による圧力への影響を考慮し，その数値解を計測データに相当する基準解と定義す
フィードバック信号s
連続の式から導出される圧pを，る．一方，超音波計測融合シミュレーションでは，力方程式に生成項として加える．それぞれの値は以上流境界条件として一様平行流を与える．基準解の下の式で求められる．
情報を内部にフィードバックしながら行う超音波計(1)f v = − K v ρ (Vc − Vs ) u ′max ∆S
測融合シミュレーションの結果は，基準解に収束す(2)s p = − K p ρ (Vc − Vs ) ∆Sることが期待される．ここで，いずれの数値シミュ
レーションにおいても，下流端は自由流出，壁面で
ここで，Kv，Kpはゲイン(無次元)，ρは流体の密度，はすべりなしの条件を与え，時間刻みは0.01 sとした． V ，V はそれぞれ基準解および超音波計測融合シミュ
sc
基礎方程式は有限体積法を用いて離散化され，それレーションによるドプラ速度， ∆Sはコントロールボ
リューム界面の面積，u’maxは流量最大時における上流境界の流速の空間平均値を表す．
2.2 エイリアジングの検出および補正方法
エイリアジングは，超音波ドプラ計測におけるドプラ偏移周波数が超音波のパルス繰り返し周波数の2分の1倍を超えたとき（計測可能なドプラ速度の最大値と最小値Vmax，Vminを超えたとき）に起こり，実際の速度とは異なる方向や大きさをもつ速度であるとみなされる現象である．例えば，ある格子点におけるドプラ速度VがV > Vmaxである場合，その位置におけるドプラ速度はVmin+(V-Vmax)であるとみなされる．V < Vminの場合も同様の現象が起きる．本研究では，基準解の計算結果に対し，Vmax，Vminの値をそれFig. 3 Computational grid and monitoring region inwhich monitoring points for feedback are arranged
Table 1 Computational conditions
Heart rate
Cardiac cycle T

0.87 Hz(52 beat/min)
1.15 sCardiac output

9.17×10-5 m3/s
(5.5 l/min)

Entrance flow

6.42×10-5 m3/s
(3.85 l/min)
0.74 m/s
28.25×10-3 m
4.0×10-6 m2/s

Maximum mean velocity u’max
Entrance vessel diameter D
Kinematic viscosity ν(a) Nonaliased data(b) Aliased data
Fig. 5 Comparison of color Doppler images ofnonaliased and aliased standard solution at peak flow74

ぞれ±25cm/sであると定義し，エイリアジングを作為は，これまでの研究結果に基づき，(Kv, Kp) = (0.4, 0.4)的に発生させた．Figure 5に，エイリアジング導入すに決定した[8]．また，通常の数値シミュレーション
る前後の基準解を比較して示す．計測可能領域を超において，その解が周期的な解に落ち着くまでに 26
えるドプラ速度は，それとは反対方向の速度とみな周期を要したことから，計算結果の評価は，すべてされるため，カラードプラ画像では，不適切な色付
26 周期目の結果を用いて行った．
けがなされる．複数の時刻におけるフィードバック
エイリアジングを含む基準解を用い，Correction B
領域内のカラードプラ画像をFig. 6(a)に示し，エイリを用いて超音波計測融合シミュレーションを行ったアジングが導入された領域をFig. 6(b)にボリューム場合について，計算開始から 1 心周期目の，Fig. 6(a)レンダリングして示す．Figure 6(a)に示す断面のうち， と同じタイミングにおけるフィードバック量の大きt = 0.15 s（収縮期の加速時）および0.40 s（収縮時のさを Fig. 6(c)に示す．Figure 6(a)から，エイリアジン減速時）の画像にはエイリアジングが生じている．
グが生じていた領域において，フィードバック量がエイリアジングの検知の方法としては，フィード大きくなっていることがわかる（図の矢印参照）．
バック量fvを指標として用い，その値が，fmin ≤ fv ≤ fmax
フィードバック量の大きさを指標とし，前節に述べた方法を用いて検知されたエイリアジングの領域を，の範囲にない場合を，異常なデータであるとして検知することにする．ここで，式(1)を参考にし，fmin =Fig. 6(d)にボリュームレンダリングして示す．エイリ
KvρVminu’max∆S，fmax = KvρVmaxu’max∆Sである．検知しアジングを作為的に導入した領域[Fig. 6(b)]と，フィたドプラ速度に関しては，2通りの補正を行う．一つードバック量を指標としてエイリアジングを検知しは，異常であると検知した基準解のドプラ速度Vsを， た領域[Fig. 6(d)]は完全に一致している．すなわち，数値シミュレーションの結果Vcに置き換え，フィーフィードバック量を指標として用いることにより，エイリアジングを検知することができた．
ドバック量をゼロにする方式である（Correction A: Vs= Vc）．もう一つは，エイリアジングの影響を考慮し，検知したエイリアジングに対し，2 通りの補正を行適正な基準解のドプラ速度Veを推定し，フィードバった超音波計測融合シミュレーションの結果を評価するため，以下の平均誤差ノルムを定義する．
ックに用いる方式である（Correction B: Vs = Ve）．ここで，Veは以下の式で与えられる．
11
en (a ) =acn (t ) − asn (t ) dt
(4)⎧⎪Vmax − (Vmin − Vs ) ( f v < f min < 0)N n amaxT T
Ve = ⎨(3)⎪⎩Vmin − (Vmax − Vs ) (0 < f max < f v )ここで，a は速度成分 u, v, w, V，圧力 p および速度ベ∑3∫

本研究におけるフィードバックの最適ゲインの値

クトル u を表し，amax は各パラメータに対する規格化のための値，T は 1 心周期，n はフィードバックの代表点番号，N はその総数，| · |はスカラー値に対し(a) Color Doppler images at t = 0.15, 0.40 and 0.80s(b) Volume rendering of aliasing region

計算結果および考察(c) Feedback intensity in the first cardiac cycle of UMI simulation
with aliasing Color Doppler images at t = 0.15, 0.40 and 0.80s(d) Volume rendering of the region where aliasing is
detected by UMI simulation with correction B

Fig. 6 Distribution of variables in the monitoring region in one cardiac cycle

75Table 2 Comparison of average error normvu
Standard solution
Ordinary simulation
UMI simulation using aliased data
- with no correction- with correction A- with correction B
UMI simulation using nonaliased data

010101

2.758
0.447
0.369
0.369

1.625
0.212
0.162
0.162

1.373
0.180
0.136
0.136

29.086
1.787
1.031
1.031

Proc. R. Soc. Lond., 117, pp. 109-159, 1979.
[2]

D.N. Ku, D.P. Giddens, C.K. Zarins, and S. Glagov, Pulsatileflow

andatherosclerosis

inthe

humancarotid

bifurcation―Positive correlation between plaque locationand low and oscillating shear stress, Atherosclerosis, 5, pp.
293-302, 1985.
[3]

D.P. Giddens, C.K. Zarins, and S. Glagov, The role of fluid
mechanics in the localization and detection of atherosclerosis,
J. Biomech. Eng., 115, pp. 588-594, 1993.

[4]

H.S. Harpual, B. Talei-Faz, L. Ansdell, E.J. Chaloner, A.
Sweeny, A. Grass, and M. Adiseshiah, Aneurysm wall shear
stress and tendency to rupture are features of physical wall
properties: an experimental study, J. Endovasc. Ther., 9, pp.
665-675, 2002.

[5]

K. Ferrara and G. Deangelis, Color flow mapping,
Ultrasound Med. Biol., 23, pp. 321-345, 1997.

[6]

D.A. Steinman, Image-based computational fluid dynamics
modeling in realistic arterial geometries, Ann. Biomed. Eng.,30, pp. 483-497, 2002.

結言

本研究では，血流の超音波計測融合シミュレーシ
ョンの，超音波計測において生じたエイリアジングの検知と補正に関する検討を行った．Correction B を用いることにより，エイリアジングを含んだ基準解を用いた場合でも，エイリアジングを含まない基準解を用いた場合と同じ，最良の結果が得られた．
Correction B はより精度の良い結果を提供できるが，実際の血流に対して応用する場合は，エイリアジン
グ以外の計測誤差も影響するため，Correction A の方がよりロバストな方法であると考えられる．以上のことから，本手法により，元の超音波計測の精度を向上できる可能性があることが示された．

[7]

M.J. Starmans-Kool, A.V. Stanton, S. Zhao, X.Y. Xu, S.A.M.
Thom, and A.D. Hughes, Measurement of hemodynamics inhuman carotid artery using ultrasound and computational
fluid dynamics, J. Appl. Physiol., 92, pp. 957-961, 2002.

[8]

K. Funamoto, T. Hayase, A. Shirai, Y. Saijo, and T. Yambe,
Fundamental study of Ultrasonic-Measurement-Integrated
simulation of real blood flow in the aorta, Ann. Biomed. Eng.,(to be published).

[9]船本健一, 早瀬敏幸, 西條芳文, 山家智之, 超音波計測
融合血流シミュレーションのフィードバック則に関する検討, 日本流体力学会年会 2004 講演論文集, pp.
538-539, 2004.

謝辞[10] M.S. Olufsen, C.S. Peskin, W.Y. Kim, E.M. Pedersen, A.

本研究の一部は，日本学術振興会特別研究員奨励費により行われた．また，数値シミュレーションは，東北大学流体科学研究所未来流体情報創造センターのスーパーコンピュータ，SGI Origin 2000 を使用して行われた．

Nadim,andJ.

Larsen,Numericalsimulation

andexperimental validation of blood flow in arteries withstructured-tree outflow conditions, Ann. Biomed. Eng. 28, pp.
1281-1299, 2000.
[11] S.V. Patankar, Numerical Heat Transfer and Fluid Flow,

参考文献
[1]

p0
1てはその絶対値，速度ベクトルに対してはその 1 乗
ノルム |u|+|v| を表し，c, s は超音波計測融合シミュレーションの結果とエイリアジングを含まない基準解であることをそれぞれ示している．各平均誤差ノルムの値を，通常のシミュレーションの値で規格化して Table 2 に示す．補正を行わなければ，エイリアジングを含まない基準解に対する誤差は，通常のシミュレーションにおいて生じる誤差よりも大きい値となった．一方，超音波計測融合シミュレーションに Correction A および B のいずれの補正方法を導入した場合においても，その計算精度には改善が見られた．特に，Correction B を用いた場合には，エイリアジングを含まない基準解をフィードバックした超音波計測融合シミュレーションと全く同じ結果が得られた．この場合，基準解に対する誤差は，各速度成分 u, v および V に対しては，通常のシミュレーシ
ョンの 37%, 16%および 14%にそれぞれ軽減されており，最も良い改善が見られた．ここで，圧力 p に関しては，計算精度に改善が見られないが，これは通常のシミュレーションにおいて生じる誤差自体が非常に小さいため，大きな問題ではないといえる．

4VHemisphere Pub. Corp.,Washington DC/New York, 1980.

C.G. Caro, J.M. Fitzgerald, and R.C. Schroter, Atheroma and[12] T. Hayase, J.A.C. Humphrey, and R. Greif, Mini-manual forarterial wall shear: observation, correlation and proposal of a

ROTFLO2, Dept. Mech. Eng. Rep. FM-90-1, Univ. Calif.

shear-dependent mass transfer mechanism for atherogenesis,

Berkeley, 1990.

76