2006年[ 技術開発研究助成 ] 成果報告 : 年報第20号

高速超音波3次元動態計測用演算システム

研究責任者

田村 安孝

所属:山形大学 工学部 情報科学科 教授

共同研究者

赤塚 孝雄

所属:山形大学 工学部 応用生命システム工学科 教授

共同研究者

湯浅 哲也

所属:山形大学 工学部 応用生命システム工学科  教授

共同研究者

柳田 裕隆

所属:山形大学 工学部 情報科学科 助手

概要

1. はじめに
動いている3次元の生体組織を高速で連続して観測できる超音波画像計測システムの画質向上と高速化の手法を確立する。現在、超音波の2次元センサアレイ、広い範囲の情報を高速に観測する符号化送信、演算処理により画像を形成する合成開口処理という3つの手法を組み合わせた高速撮像システムの研究を進めている。本研究では、対象の画質の向上と計算の高速化のための演算方式を確立しシステムを構築して実験的に検証する。

2. 研究の概要
2.1 研究の必要性
体内の臓器や胎児などは3次元の形状を持ち常に運動している。これらの対象を少ない負担で連続して観測できる3次元超音波診断装置が求められている。しかし、従来の方式では超音波送信と受信を多数回繰り返す必要があるため、3次元の動態撮像への適用は困難であった。申請者らは、超音波センサを2 次元に配置した2次元センサアレイ、超音波を符号により変調して送信する符号化超音波、数値的に画像を形成する合成開口技術の3つを組み合わせた高速3次元撮像システムを提案し、そのプロトタイプを試作した1)、2)。これを用いた実験では、生体組織を模擬した評価用の対象(超音波ファントム)や生体組織の3次元観測が可能であることを確認した3)。しかし、医用診断への適用には、画質の向上と像再生のための計算の高速化が不可欠であることが指摘されている。2次元センサアレイ、合成開口、符号化超音波の組み合わせは、次世代の超音波診断装置の手法として有力であり、海外では報告例も多い。画質の改善と計算の高速化の可能性を示すことは、日本におけるこの分野の競争力の強化と早期の実用化のために急務である。

2.2 解決しようとする課題
(1) 画質の向上
反射率の微妙な違いを描出するため、画像のコントラストの向上が課題となる。合成開口方式では超音波素子の配置間隔が疎らな「疎アレイ」を使用し、ビーム形成に関わる素子の個数も多くできない。これが画質の低下の主な原因である。また、高速なデータ収集のために採用した符号化超音波も信号の歪みなどがあるとコントラストの低下につながる。
(2) 計算の高速化
数値化されたデータから計算により像を形成する方式では、計算システムの性能とコストが常に問題となる。特に、動態の計測をする場合、データ量が増えるため高速計算システムは必須のものとなる。

2.3 提案する解決策
画質の向上には、時間軸方向の情報を利用した数値的なフィルタリングにより対象の動きを検出しコントラストを向上する像再生の手法を開発する。計算の高速化は2つのアプローチを取る。一つはデータを高速に処理するため、複数のパーソナルコンピュータ(PC)を高速ローカルエリアネットワーク(LAN)で相互接続して並列動作させる「並列PC クラスタ」の構築である。市販のユニットを組み合わせて16 ノードまでのLinux ベースのシステムを構築し、この上で信号処理ソフトウェアを開発する。もう一つの方法は、演算装置のハードウェア化である。本研究の期間では、ハードウェア化に向いた計算方式を開発する。

3. 三次元画像系列の画質向上
3.1 検証用システム
本研究では、小規模な2次元アレイセンサを用いる現有の検証用システム(図1)を用いた2)。このシステムは64 個までの送信素子を駆動し、32 個までの受信素子の受信データを収集できる。超音波周波数は2MHz である。アレイセンサは約1mm 角の素子を40×40(1600 個)のマトリクス上に配置し,その内の144 個に結線したものである。

3.2 像再生の原理
この撮像システムでは,2値の直交関数系であるWalsh 関数で変調した送信信号を用いる。位置xTj( j ? 0, ... ,NT ?1)に置かれたNT 個の送信素子を,互いに異なるWalsh 関数の1周期で変調された信号で駆動する。ここで, NT は2の冪乗( k NT ? 2 )とする。対象からのエコーを記録すると,Walsh 関数と送信素子の組み合わせを変更して次の送受信サイクルを実行する。これを必要な回数だけ繰り返す。今,送信サイクルに番号付けしてp ? 0,1,2,?とする。p 回目の送受信サイクルでj 番目の送信素子を駆動する信号を( ) ( ) u t pj とし,対象からのエコーを位置xRi(i ? 0, ... ,NR ?1)に置かれたNR個の受信素子で受信して得られる信号を( ) ( ) r t pi とする。使用するWalsh 関数は,から始めて,Sylvester の展開,をk ?1回繰り返して得られるS 型Hadamard 行列の各行の要素をクロック信号に同期して読み出すことで発生される。このWalsh 関数を,これと同期する超音波信号と乗算することで変調信号を得る。使用する超音波の中心周波数をf0とし,このクロック信号の周期?tは超音波の周期1/ f0 の整数倍に設定する。送信素子とWalsh 関数の対応関係は,送信サイクルの番号と送信素子の番号を2進数表現し,そのビット毎の排他的論理和により対応するHadamard 行列の行番号を決定する方法を採用している。このとき, ( ) ( ) u t pj は,で与えられる。ここで, ? は2つの数値をk ビットの2進数としてビット毎の排他的論理和を取る演算を示すものとし,g(t)はで表現される幅?tの正弦波パルスとする。p 番目の送受信サイクルでの,空間内の位置xの画素を求める操作は次式のように表すことが出来る。となる。ここで,は,i 番目の受信波形( ) ( ) r t pi とj 番目の送信波形( ) ( ) u t pj の相互相関関数であり, ( ) ( ) r t pi を( ) ( ) u t pj にマッチした相関フィルタで処理した出力と解釈できる。送信波形( ) ( ) u t pj が互いに直交していて,しかも各々の自己相関関数が鋭いピークを持つ広帯域信号であれば,式(5)は送信素子を独立に駆動して得られたNT × NR 個のエコーデータを用いる合成開口による焦点形成の演算と解釈できる。送信と受信のビーム形成はエコーデータに対して数値的に実行される。この方式では,1回もしくは少数回の送受信で3次元像再生に必要なデータを収集できるため高速の撮像が可能である。また,送信及び受信の両方で,全ての奥行に対して焦点を形成できる。

3.3 画質の向上 4)
このシステムでは1回の送受信毎に3次元像が再生できる。しかし,得られた画像系列のフレーム毎の画像は比較的アーチファクトレベルが高いものになる。これは,広いダイナミックレンジを必要とする応用で問題となる。対象が静止している場合は、式(3)により送信素子とWalsh 関数との対応関係を変更しながら繰り返し送受信を行い、得られた画像を重ね合わせることで、アーチファクトの相対レベルを低減できる。一方、対象が運動している場合,対象と送信素子・受信素子の間の相対的な運動により,画像系列の同じ空間的位置での位相が時間(送受信サイクル番号)とともに回転する。この現象は,パルスドプラ法で見られるものと同じである。このため,相対速度が大きくなると,単なる画像の加算だけでは反射体が存在する位置での画素が同位相で加算されなくなり、アーチファクト低減の効果が得られなくなる。そこで,対象の移動速度を推定し,その影響を補償した重ね合わせを行う必要がある。本研究で検討したのは,画素時系列をフーリエ変換し,得られたピークの値から動き補償された画像を得る方法である。この方法では,ある幅を持った解析用時間窓の中で複素の画素時系列をフーリエ変換する。与えられた周波数範囲の中でフーリエ変換の絶対値の最大値IMAXと,そのときの周波数fMAX を像再生の空間格子の各点で求める。得られたIMAX(x)が解析用窓の中で動き補償された画像になり,fMAX (x)が合成されたビームに沿った方向の速度の推定値となる。これは,パルスドプラ法を我々の高速3次元システムで実行したものと考えることができる。

3.4 実験とシミュレーション
水中や生体組織を対象とし,奥行150mm の場合,送受信の時間間隔は0.2ms まで短くできる。つまり,原理的には,150mm の距離までの動く対象を1秒間に5000 コマの最大速度で観測することが可能である。ただし,試作したシステムの現在の仕様では,送受信の時間間隔の下限は20ms となっている。送受信の間隔を20ms に設定し,運動する物体の撮像を行った。生体類似の対象として,運動するバルーンを用いた。直径10mm のバルーンを水,超音波プローブ用ゼリー,ポリマー粉末の混合物を満たしポンプにより出し入れすることで膨張収縮運動をさせている。図2-a は1 回の送受信により得られた画像であるが,アーチファクトの中に埋もれていてバルーンの形状は認識できない。一方,得られた3次元像の時系列(16 フレーム)からドップラーシフトを検出して得られた図2-bの画像では,バルーン表面が描出されている。検証用システムでは取得できるデータ量に上限があるため、評価用実験装置による実験を想定してコンピュータシミュレーションを行った。条件は実験と同様に設定し、対象としては、微小な点状の反射体が集まって球体となり、膨張収縮運動を繰り返しているものを想定し、受信波形を計算して3 次元画像系列を再生した。球体のサイズなどは実験で用いたものと同程度となるように設定した。また、送信回数は、運動の様子が一通り観測できるよう、バルーンの膨張・収縮サイクルの約1/2 周期となる188 回( 187 × 20[ms]=3.74[s])とした。以上のようにして再生した3 次元画像系列の一部を図4に示す。ただし、送信20 回毎にサンプルしている。図3からは、アーチファクトが多く見られるものの、バルーンが最小の状態から膨張し、最大の状態に近づいていく動きが確認できる。以上のような画像系列に対して、提案した手法を用いて動き補償を行った結果を示す。図4は、p=15 から送信20 回毎にサンプルした画像である。画質系列の画質が向上していることがわかる。

4. 計算の高速化
4.1 PC クラスタによる高速演算システム 5)~7)
PC クラスタとは、並列に接続された複数のPCによりプログラムを実行するシステムである。安価なコストでスーパーコンピュータ並みの計算能力を発揮できる事が特徴である。本研究では、超音波の3次元高速撮像システム用のPC並列クラスタを構築し性能を評価した。並列像再生プログラムは、通常の像再生プログラムにMPI(Message Passing Interface)を実装し、各PC において同一のプログラムにて実行される、SPMD(Single Program Multiple Data)タイプとして作成した。MPI では、各PC へとrank という値が定義される。rank の値が0 のPC をroot rank と言い、PC クラスタの中で親ノードと呼ばれ、1 以上のrank の値を持つPC は子ノードと呼ばれる。MPI関数はこのrank の値を用いることにより、通信制御を行っている。

4.2 演算システムの構成
新しく下記のスペックを持つPC8台を、すべて同じ構成で構築した。CPU:Athlon64 3000+Mem:PC3200 512MB*2OS:Fedora Core3 MPICH-1.2.6像再生プログラムは大きく見て相互相関・遅延加算の2つの計算プロセスが存在する。本研究の並列像再生プログラムでは以下の流れで処理を行っている。
1) 親ノードから子ノードへと初期データを転送
2) 各ノードにて相互相関を行う。分割方式は照射回数分割。
3) 各ノードの値においてバタフライ演算を行い、値を統合していく
4) 各ノードにて遅延加算を行う。分割方式は空間分割。奥行き毎、各ノードが計算。
5) 値をroot rank へと統合。画像を生成。
MPI に関するプログラムの流れを示すと、
1) 初期値をMPI_Bcast エコーデータをMPI_Scatter にて、各ノードへと転送
2) 各ノードにて相互相関をし、結果をMPI__Allreduce にてバタフライ演算
3) 各ノードにて遅延加算をし、結果をMPI_Gather にて、root ノードへと転送
4) 出力はroot ランクで行う。
となる。計算の配分は、総てのPC に均等に計算させるように設定。素子数が増加することにより発生する、メモリ枯渇、スワップを最小限に抑えるべく、メモリの確保をすべて最初に行うのではなく、必要なときに必要なだけ確保するようにした。

4.3 シミュレーション
PC の台数を1台から8台へと増やすことで、プログラムの並列化の効果を確認した後に、PCの台数を8台と固定し、シミュレーション用データを用いて素子を32ch から増やしていく。これにより、実用的な受信素子数でのプログラムとPC クラスタの限界を測る。演算時間は、各PC の処理時間の平均値をとり、それを3回実行した平均値で評価した。送信素子数を32 に固定し、受信素子数を32 から1024 まで変化させて計算時間を評価した結果を表1および図5 に示す。送信素子32 個・受信素子512 個までは、安定した動作時間を確認でき、受信素子の増加に対して演算時間は比例している。しかし、受信素子1024 個において、reduce の実行時間が急激に増加している。この原因は、各ノードで発生したデータのswap によるものだと考えられる。実際に、最大1GB のswap を確認できた。これが異常な遅延の原因ではないかと考えられる。

5. 計算のハードウェア化の検討 8)
従来のシステムでは、像再生アルゴリズムはC言語でコーディングされたソフトウェアとして実装され、CPU により像再生演算がなされている。しかしこのシステム構成では像再生の処理速度はCPU の性能に依存してしまい、高性能なCPU を用いればある程度の処理の高速化は期待できるが効果的な高速化は期待することができないうえに性能あたりのコストが増加してしまうので限界がある。そこで効果的な高速像再生処理を実現するために、像再生演算アルゴリズムのハードウェア実装を検討した。

5.1 像再生のアルゴリズム
ハードウェア化に適した像再生手法として、マッチドフィルタによる手法を検討した。この手法では、像再生の格子点はアレイの中心から扇形に広がる直線上に配置される。同一直線上の格子点での画素値は同一のフィルタ演算により求められるため、演算の高速化が可能である。図6に従来の手法による再生像とマッチドフィルタによる再生像を示す。対象はマーカとして直径1mm の針を挿入した20cm 角の豚肉ブロック内部の像である。マッチドフィルタの像では、針の像がボケている。これは、奥行の分割方式を変更することで改善可能と考えている。この像再生手法のハードウェア化を検討し、設計を進めている。図7 は設計したハードウェアの構成である。

6. まとめ
本研究の成果は以下の通りである。
(1)動きの情報をドプラシフトとして検出し、周波数領域でのピークの値から画素を求める方法により、画像系列のコントラストが改善できることを確認した。
(2)PC クラスタにより、受信素子数を512 まで増加しても対応できることを確認した。
(3)ハードウェア化に適した像再生手法を検討し、ハードウェアの設計に着手した。
運動する3次元物体のドプラシフトを視野内で同時に観測できるシステムは世界的にも例が無いと考えられる。ドプラシフト検出による画像系列の画質向上は簡易な手法ではあるが、ある程度の効果があったので、より洗練された動き検出・補償の方法に発展させたい。クラスタはシミュレーションや実験データの処理に活用し、最終的には演算装置のハードウェア化により、実用化に結び付けたいと考えている。