WO2013180269A1

WO2013180269A1 - 超音波撮像装置

Info

Publication number: WO2013180269A1
Application number: PCT/JP2013/065196
Authority: WO
Inventors: 鱒沢　裕; 慎太高野; 貞一郎池田; 栗原　浩
Original assignee: 日立アロカメディカル株式会社
Priority date: 2012-05-31
Filing date: 2013-05-31
Publication date: 2013-12-05
Also published as: JP2015163087A

Abstract

　医用超音波診断装置において、周波数に応じて変化する被検体媒質中の不均質特性に即した適応ビームフォーマの制御技術を提供する。　被検体からの超音波エコーを受信し、所定の受信焦点の位置に応じて受信信号を遅延させて第一の時間信号群を形成する。複数の部分帯域適応処理部２２－１～２２－Pは、それぞれ、第一の時間信号群から所定の周波数帯域の第二の時間信号群を抽出し、第二信号群に適応処理を行って適応重みを求め、適応重みにより第二の時間信号群を重み付けした後加算する。これにより、広い帯域の受信信号であっても、周波数帯域ごとに適切に適応処理することができる。

Description

超音波撮像装置

　本発明は、超音波を用いて被検体内の画像を撮像する超音波撮像装置、特に医用超音波撮像技術に関する。

　人体の各部位をはじめとする被検体の内部を、非侵襲的に画像化する医用超音波撮像装置が広く知られている。この装置では、圧電素子などの電気音響変換手段を複数内蔵する超音波探触子を被検体に接触させ、超音波パルス波を送受信することにより被検体の内部の情報を得ている。超音波撮像装置内の送波ビームフォーマは、超音波探触子内の複数の電気音響変換手段に対して、適切な駆動電気信号を受け渡し、生体内へ超音波パルスの送信を行わせ、所定の音響的伝搬効果を生じさせる。また、超音波パルスが被検体内で反射されることで生じたエコーは、再び超音波探触子で受信される。超音波探触子は、内蔵する複数の電気音響変換手段で超音波エコーを電気信号に変換する。超音波撮像装置内の受信ビームフォーマは、超音波探触子で得られた複数の受信信号を、所定の音響的伝搬過程の想定のもとに遅延処理や振幅重み付け処理をして、被検体の内部情報を像再生に資する信号に変換する。

　従来の受信ビームフォーマは、予め想定した音響的伝搬効果をもとに複数の受信信号の各々に与える遅延時間情報や振幅重みを決定している。例えば、複数の受信信号に対して受信信号値の如何に拘らず、伝搬音速を一様と仮定し、超音波探触子の電気音響手段と被検体内部との幾何学的な位置関係を元に距離を計算する。こうした受信ビームフォーマによる処理は、遅延加算処理と呼ばれる。遅延加算処理では、被検体内部の想定した幾何学的位置に、強く音波を後方散乱するものが存在する場合、遅延後の互いの干渉によって位相が揃うため、加算後の信号振幅値が増強される。逆に、想定と異なれば、遅延後の信号群の間の互いの位相差の違いが、加算されることにより信号振幅値を減弱させる。これらの信号振幅値が、像再生に用いられる。

　この遅延加算処理は、超音波のパルス波面の伝搬飛行時間の差が位相差を生み、干渉により信号振幅値に強弱が付くことを利用している。しかし、同じ受信信号に複数位置からの反射エコー信号情報が重畳している場合には、重畳の度合いを推定して考慮に入れることはできない。適応ビームフォーマ処理は、受信信号群を確率変数とみなし、統計的推定から受信信号群各々の遅延あるいは位相、重み付けを逐次、変更したのちに加算する処理を行うことで目的とする信号を好条件で得られる。適応ビームフォーマは、ソナーやレーダー、移動体通信などのリモートセンシングの分野で広く利用されている。ソナーの送受波器、レーダーならびにアンテナの配列（アレイ）においては、異なる方向への遅延加算処理の反復が空間フーリエ変換に関連付けられ、空間フィルタ技術、空間のスペクトル解析技術として用いられている（非特許文献１,非特許文献２,非特許文献３）。非特許文献１、非特許文献２によれば、Ｃａｐｏｎ法あるいはＭＶＤＲ（Ｍｉｎｉｍｕｍ　Ｖａｒｉａｎｃｅ　Ｄｉｓｔｏｒｔｉｏｎｌｅｓｓ　Ｒｅｓｐｏｎｓｅ）法、ＡＰＥＳ（Ａｍｐｌｉｔｕｄｅ　ａｎｄ　Ｐｈａｓｅ　ＥＳｔｉｍａｔｉｏｎ）法などがアルゴリズムとして列挙できる。これらの手法は，受信データのベクトルの相関行列あるいは共分散行列の推定を経て実現される。

　超音波撮像装置では、人体腹部の深部(凡そ１０～２０ｃｍ)まで有効視野が広がる肝臓や、表在性の甲状腺などを撮像対象として選ぶ場合がある。使用される超音波パルスの周波数は、２～１８ＭＨｚであり、約３オクターブの周波数の開きがある。圧電素子などで構成される超音波探触子の電気音響変換手段の共振周波数を代表的な中心周波数とすると、上記周波数の範囲の中心周波数に対する比帯域は５０～１００％と大きく、広帯域信号として扱う必要がある。例えば、体表から１０ｃｍ程度の深さまでの被検体視野を撮像する場合には，深さに応じて受信ビームフォーマの出力周波数帯域の上限を浅部から深部に向かうにつれて漸次低周波側に移動させることで画像の描出能を改善できることが知られている。また、距離方向の分解能を向上させるために、アレイ（超音波探触子の配列された電気音響変換手段）の個々の受信信号に対し、時間軸方向の周波数帯域を変化させる受信ダイナミックフィルタ技術や、受信信号を複数（おおよそ３～５程度）の帯域通過フィルタで周波数帯域の部分空間に分解し、受信信号の検波前あるいは検波後の信号の合成を制御する周波数コンパウンド技術が知られている。

Ｂ．Ａｌｌｅｎ，　Ｍ．Ｇｈａｖａｍｉ，"Ａｄａｐｔｉｖｅ　Ａｒｒａｙ　Ｓｙｓｔｅｍｓ"，Ｊｏｈｎ　Ｗｉｌｅｙ　＆　Ｓｏｎｓ　Ｌｔｄ（２００５）Ｐ．Ｓｔｏｉｃａ，Ｒ．Ｍｏｓｅｓ，"Ｓｐｅｃｔｒａｌ　Ａｎａｌｙｓｉｓ　ｏｆ　Ｓｉｇｎａｌｓ"，Ｐｒｅｎｔｉｃｅ　Ｈａｌｌ（２００５）Ｊ．Ｇ．ＭｃＷｈｉｒｔｅｒ，Ｔ．Ｊ．Ｓｈｅｐｈｅｒｄ，"Ｓｙｓｔｏｌｉｃ　ａｒｒａｙ　Ｐｒｏｃｅｓｓｏｒ　ｆｏｒ　ＭＶＤＲ　ｂｅａｍｆｏｒｍｉｎｇ",ＩＥＥＥ　Ｐｒｏｃｅｅｄｉｎｇｓ，　ｖｏｌ．１３６，Ｐｔ．Ｆ，Ｎｏ．２，ＡＰＲＩＬ（１９８９）

　上記リモートセンシング分野の適応ビームフォーマ技術は、多くの場合、受信信号源（後方散乱エコー源）の位置が、アレイに対して無限遠に仮定されている。医用超音波撮像装置は、超音波探触子内の超音波送受信素子のアレイの物理的な代表寸法に対する、受信信号源までの距離の比（Ｆ値）があまり大きくないため、Ｆ値が無限大を仮定する適応ビームフォーマ技術を適用するには工夫が必要である。近距離で点を仮定する受信信号源から超音波探触子に到達する際のエコーの波面が凹面になるため、超音波探触子の受信信号を受信焦点を中心とする凹面に合わせて遅延させ、合焦した遅延受信信号を得る。医用超音波診断装置では、時間軸方向のパルス幅が、画像の距離方向分解能に関係し、Ｆ値と波長が、方位分解能を支配するが、さらに合焦した遅延受信信号の時間軸方向とアレイの空間軸方向の関連付けを考慮することが望ましい場合もある。

　また、従来の受信ビームフォーマは、撮像対象に対する事前情報が与えられない限りは、音速一定でかつ超音波の波長に対して十分小さな組織が均質に連続する媒質を想定して信号を処理するように設計されている。しかし、実際の被検体の内部組織は、脂肪層や筋肉などの組織ごとに音速や密度の差異があり、超音波信号が減衰する。また、被検体内に超音波波長と同程度の寸法の組織構造が存在し、超音波を散乱する場合もある。このような超音波の減衰や散乱に起因して超音波波面の到達乱れを生じ、撮像能力が低下することがある。しかも、減衰や散乱による波面の到達乱れは、超音波の周波数に依存して変化する。

　本発明の目的は、医用超音波診断装置において、周波数に応じて変化する被検体媒質中の不均質特性に即した適応ビームフォーマの制御技術を提供することにある。

　本発明では、受信信号を遅延処理した第一の時間信号群を重み付け処理する適応処理部を備え、この適応処理部が、複数の部分帯域適応処理部を有する。複数の部分帯域適応処理部はそれぞれ、第一の時間信号群から所定の周波数帯域の第二の時間信号群を抽出し、第二の時間信号群に適応処理を行って適応重みを求め、適応重みにより第二の時間信号群を重み付けした後加算する。

　本発明では、受信信号を遅延処理後の第一の時間信号群を、周波数帯域ごとの第二の時間信号群に分け、それぞれ部分帯域適応処理部が適応ビームフォーマ処理するため、周波数に応じて変化する被検体媒質中の不均質特性に即した適応ビームフォームを行うことができる。

本実施形態の超音波診断装置の全体を示すブロック図。本実施形態の受信ビームフォーマ１７の構成を示すブロック図。図２の受信ビームフォーマの部分帯域適応処理部２２－１の構成を示すブロック図。図３の部分帯域適応処理部３２（ＱＲＤ－ＲＬＳアルゴリズム構成）の構成例を示すブロック図。図４のシストリックアレイ４５の演算部構成を示すブロック図。図５のシストリックアレイのセルの演算を示す図。図３の部分帯域適応処理部３２（ＳＭＩ法等によるＭＶＤＲ）の構成例を示すブロック図。本発明の適応ビームフォーマの時空間再標本化を説明する図。部分帯域適応処理部２２－１の時空間スナップショットとアンサンブル平均を説明する図。部分帯域適応処理部２２－Ｐの時空間スナップショットとアンサンブル平均を説明する図。可干渉性損失行列の対角成分と対角位置変数ｑを説明する図。可干渉性損失関数Ｃ_ｉおよび有効干渉長λ_ｉの部分帯域による違いを説明するグラフ。有効干渉長λ_ｉの部分帯域周波数による違いを説明するグラフ。可干渉性損失関数Ｃ_ｉおよび有効干渉長λ_ｉの受信開始からの経過時刻による違いを説明するグラフ。有効干渉長λ_ｉの経過時刻による変化を説明する図。可干渉性損失行列Ｂ_ｉへの口径位置変数ｑ’の設定を説明する図。口径制限関数Ｄ_ｉおよび有効口径幅α_ｉを説明する図。有効口径幅α_ｉの部分帯域周波数による違いを説明するグラフ。口径制限関数Ｄ_ｉおよび有効口径幅α_ｉの経過時刻による変化を説明する図。部分帯域ごとの有効口径幅α_ｉの経過時刻変化を説明する図。不均一散乱シミュレーションに用いた生体モデルを説明する図。シミュレーションの音源の送信波形を説明する図。シミュレーションの音源の周波数スペクトルを説明する図。不均一散乱の効果を説明する図。ＡＰＥＳ適応ビームフォーマの効果を説明する図。部分帯域ＡＰＥＳ適応ビームフォーマの効果を説明する図。部分帯域ＡＰＥＳ適応ビームフォーマのコヒーレント合成効果を説明する図。部分帯域ＡＰＥＳ適応ビームフォーマのコヒーレント合成効果を説明する図。部分帯域ＡＰＥＳ適応ビームフォーマのコヒーレント合成効果を説明する図。部分帯域ＡＰＥＳ適応ビームフォーマのコヒーレント合成効果を説明する図。部分帯域ＡＰＥＳ適応ビームフォーマのコヒーレント合成効果を説明する図。

　本発明では、被検体からの超音波エコーを受信する複数の受信素子と、所定の受信焦点の位置に応じて複数の素子の受信信号をそれぞれ遅延させて第一の時間信号群を形成する遅延部と、第一の時間信号群を重み付け処理する適応処理部とを有する。適応処理部は、複数の部分帯域適応処理部を有する。複数の部分帯域適応処理部はそれぞれ、第一の時間信号群から所定の周波数帯域の第二の時間信号群を抽出し、第二の時間信号群に適応処理を行って適応重みを求める。そして、この適応重みにより第二の時間信号群を重み付けした後、加算する。

　これにより、周波数帯域ごとに適応処理を行うことができるため、周波数に応じて変化する被検体媒質中の不均質特性に対応して、適応ビームフォームを行うことができる。

　なお、超音波波面の到達乱れがある場合には、適応アルゴリズムで想定する信号源のエコー信号の空間分布が、多重経路を経由して伝搬していることも想定に含めることが望ましい。すなわち、単一点の反射源からの理想的球面波のようなエコー波面が、空間軸方向と時間軸方向に分散した複数の小波面群の重畳として到達する場合をも考慮する。さらに、これらの到達乱れが超音波のパルスの波長すなわち時間周波数により異なった変化をすることを考慮することが望ましい。

　上述の複数の部分帯域適応処理部のうち１以上は、抽出する周波数帯域の一部または全部が、他の前記部分帯域適応処理部が抽出する周波数帯域と重なっていてもよい。

　複数の部分帯域適応処理部は、相互に独立して適応処理を行う構成とすることができる。これにより、周波数帯域ごとに独立して適応処理を行うことができるため、生体の非線形性による調波の発生に即した適応処理もできる。また、超音波パルスの送信を複数帯域に分けて指向性や焦点距離、方向を異ならせる周波数多重送信を周波数ごとに適切に適応処理することが可能になる。

　複数の部分帯域適応処理部は、相互に独立して、第二の時間信号群をそれぞれ時間軸方向および空間軸方向の少なくとも一方について再標本化した後、前記適応処理を行う構成としてもよい。

　それぞれの周波数帯域に応じて適切な標本化間隔の第二の時間信号群を生成できる。超音波探触子のアレイを構成する受信素子の指向性は、音波の放射面の寸法と波長の関係で決まる。波長は帯域を代表する周波数と、音速に依存する。撮像すべき被検体の超音波探触子近傍では受信ビームフォーミングに同時に関与すべき受信素子幅（受信口径）は、指向性の主極の角度範囲に受信焦点が含まれるよう周波数帯域に応じて幅を変化させることが望ましい。同じ焦点に対して低周波成分は大口径、高周波成分は小口径となる。数オクターブにわたる広帯域受信を行う探触子では、グレーティングローブを抑圧するために高周波数成分を考慮しておよそ０．５～１．５波長の適切な素子間隔に細分化されている。高周波成分を基準とすると、低周波数成分に対しては過度に細分化された状態となる。実際の素子のアレイより大きな間隔の仮想的なアレイを作り出すよう空間軸方向の再標本化を行って低周波成分に対して実質的に受信信号数数を低減できる場合がある。また、単一の広帯域処理では高周波数成分に対応して時間エイリアシングを起こさない時間標本化周期が設定されるが、時間軸方向の再標本化をより狭い部分周波数帯域ごとに行うことにより、周期を長くでき単位時間当たりの標本数を低減できる場合がある。このように、部分周波数帯域ごとに受信口径の制限、空間軸方向や時間軸方向の再標本化を行うことで装置演算資源を効率よく割り当てられる場合がある。

　複数の部分帯域適応処理部は、第二の時間信号群をもとに相関行列あるいは共分散行列を推定する。行列の推定にはアンサンブル平均を取るが、次数を第二の時間信号群の空間軸方向の標本点数より少なくした場合には、空間軸方向にもアンサンブル平均を取ることができる。時間軸方向および空間軸方向の少なくとも一方についてアンサンブル平均をとる適応処理を部分周波数帯域ごとに独立して構成してよい。また、アンサンブル平均の条件を第一の時間信号群の受信開始からの経過時間とともに変化させてよい。

　第二の時間信号群から相関行列等を推定する場合に、十分に安定な（確率収束した）期待値を得るのに適切な時間軸方向および／または空間軸方向のアンサンブル平均数や行列次数は、周波数帯域ごとに異なる場合がある。本発明では、第二の時間信号群を形成する上記構成によりこれを実現できる。また、アンサンブル平均の条件を第一の時間信号群の受信開始からの経過時間とともに変化させることにより、超音波探触子に対して被検体の浅い位置と深い位置とで、十分に安定な適応処理結果を得る条件を調整することができる。

　複数の部分帯域適応処理部は、相互に独立して、第二の時間信号群から相関行列あるいは共分散行列を求め、行列要素を重み付けした後、適応処理を行うことも可能である。行列要素に重み付けする処理は、第二の時間信号群に空間軸方向の標本点に対応して所定の重みづけベクトルの要素を空間方向の標本点ごとに乗じてから相関行列あるいは共分散行列を推定する方法と、相関行列あるいは共分散行列を先に推定してから、所定の可干渉性損失行列（ｃｏｈｅｒｅｎｃｅ　ｌｏｓｓ　ｍａｔｒｉｘ）との行列要素位置同士の積(Ｓｃｈｕｒ－Ｈａｄａｍａｒｄ積)を行う二通りが考えられる。前者では、重みづけベクトルの要素値を非負とし、ベクトルの中央要素位置の値が最も大きく両端の要素に向かって減少させ、所定値以上離れた要素ではゼロになるようにする。後者では可干渉性損失行列の行列要素の値を、非負とし、主対角成分の値が最も大きく、主対角成分から離れる副対角成分ほど小さくなり、主対角成分から所定値以上離れた副対角成分の行列要素はゼロであるように定める。また、複数の部分帯域適応処理部は、それぞれ独立して、第一の時間信号群の受信開始からの経過時間とともに、可干渉性損失行列の行列要素をゼロにする所定値を変化させることも可能である。

　このように相関行列あるいは共分散行列の行列要素を重み付けすることにより、周波数帯域の高低や、被検体の受信焦点の深さに応じて信号に重畳しやすいノイズや生体不均一の影響を抑制して、適応重みを求めることが可能になる。また、撮像すべき被検体の深度に応じて受信ビームフォーミングに関与すべき受信口径幅を第一の時間信号群の受信開始からの経過時間に応じて段階的に変化させる受信可変口径技術を用いる場合がある。受信口径幅に含まれる受信素子の数が少なく、口径幅の素子幅単位の変化に伴い、適応重みが経過時間に対して急激に変化する場合がある。相関行列あるいは共分散行列の行列要素経過時間に応じた連続的で滑らかな重み関数を元に重み付けすることにより、適応重みも滑らかに変化させることができる。

　複数の前記部分帯域適応処理部の出力にそれぞれ個別に、時間に依存する異なる重みを乗じて加算できる。加算は位相情報を保存するコヒーレント加算あるいは位相を保存せず包絡線振幅情報だけを加算するインコヒーレント加算部を有する構成としてもよい。

　複数の部分帯域適応処理部は、それぞれ独立した条件で、第一の時間信号群の受信開始からの経過時間とともに、相関行列あるいは共分散行列の次数（第二の時間信号群の空間軸方向の再標本化出力信号数）を低減し、第二の時間信号群について、時間軸方向および空間軸方向の少なくとも一方について求めるアンサンブル平均を行う標本数を増加させてもよい。

　これにより、受信開始からの経過時間とともに第二の時間信号群の信号から選ぶ行列次数を低減することにより、反射エコーへの生体不均一伝搬の影響が浅部から深部にわたり周波数帯域に応じて異なる変化をすると想定する場合に、対応することができる。また、相関行列あるいは共分散行列の次数を低減した場合にも、アンサンブル平均を行う標本数を増加させることで、安定な（確率収束した）適応処理を行うことができる。

　複数の部分帯域適応処理部は、それぞれ独立した条件で、第一の時間信号群の受信開始からの経過時間とともに、抽出する第二の時間信号群の信号数を低減することにより、相関行列あるいは共分散行列の次数を低減し、かつ、次数の低減に伴って前記重みづけベクトルや可干渉性損失行列の行列要素の値を変化させることも可能である。

　これにより、周波数に応じて反射エコー信号への生体不均一の影響が浅部から深部にわたり周波数に応じて異なる変化をする場合に、相関行列あるいは共分散行列の次数を低減することにより、演算資源（たとえば処理演算装置や回路の最大単位時間当たり演算能力）を有効に利用することができる。また、可干渉性損失行列を行列要素の値を変化させることにより、次数変化に伴う急激な適応重み変化を抑制できる。

　以下、具体的に本発明の実施形態について説明する。

　本発明を実施する超音波撮像装置の全体構成を図１を用いて説明する。超音波撮像装置は、操作者からの入力を受け付ける入力部１１と、制御部１２と、送信ビームフォーマ１３と、送信回路１４と、送受分離回路１５と、超音波探触子１０と、受信回路１６と、受信ビームフォーマ１７と、バックエンド処理部１８と、表示部１９とを備えて構成される。

　入力部１１は、スイッチ群やキーボード等からなり、操作者からの入力を受け付ける。制御部１２は、所定の時間間隔で超音波送信開始指令を送信ビームフォーマ１３に通知する。送信ビームフォーマ１３は、送信回路１４に複数の送信制御信号群を出力する。送信制御信号群は、送信開始時刻、送信パルス波形、送信振幅情報などを含む。送信回路１４は、必要に応じて図示しないディジタル－アナログ変換器などを用いて、送信制御信号群の情報に対応するアナログ連続電圧波形を生成し、内蔵する電力増幅回路によりこれを増幅して高電圧パルス信号群を生成し、送受分離回路１５に出力する。送受分離回路１５は、送信回路１４から受け取った高電圧パルス電圧信号群を超音波探触子１０に受け渡す。

　超音波探触子１０は、電気音響変換器アレイ１０ａを備えている。電気音響変換器アレイ１０ａは、電気信号（電圧波形）を機械応力信号（音圧波形）に変換する圧電体などで構成されたトランスデューサを複数配列した構成である。電気音響変換器アレイ１０ａの複数の受信素子は、送受分離回路１５から受け取った高電圧パルス信号群によりそれぞれ駆動され、送信パルス１０ｂを被検体内部に照射する。

　送信パルス１０ｂが、被検体内の反射体１０ｃにより後方散乱されることにより、パルスエコー１０ｄが発生する。パルスエコー１０ｄが超音波探触子１０に到達すると、電気音響変換器アレイ１０ａの複数のトランスデューサの各々により、パルスエコー１０ｄの音圧波形が電圧波形に変換されて受信電圧信号群となる。以下、受信時においては、トラスデューサを受信素子と呼ぶ。

　受信回路１６は、受信電圧信号群に対して、増幅、時間利得補償、アナログーディジタル変換等の処理を行う。これにより、受信電圧信号群は、離散時間で標本化されたディジタル受信信号群に変換される。本実施形態では以下、特に説明のない限り標本化されたディジタル受信信号群のことを単に「受信信号群」と呼称するものとする。

　受信信号群は、受信回路１６から受信ビームフォーマ１７に出力される。本実施形態では受信信号群はＮチャンネルであるとする。受信ビームフォーマ１７は、受信信号群から求めた相関行列あるいは共分散行列を元に推定した振幅重み（適応重み）を受信信号群に与える適応ビームフォーマ処理能力を有する。すなわち、受信ビームフォーマ１７は、受信信号群を構成するＮチャンネルの信号の各々に対して、時間遅延および／または搬送波の位相回転処理を行った後、適応重みにより重み付けを行った後、総和する。受信ビームフォーマ１７は、総和した信号をバックエンド処理部１８に出力する。また適応処理のための情報は制御部１２とのデータバス１２－１を介して行われる。

　バックエンド処理部１８は、受信ビームフォーマ１７の出力の直交検波、振幅対数圧縮する処理を行う。また、複数回の送受信における受信ビームフォーマ１７の出力を、検波前あるいは振幅対数圧縮後の状態で内蔵する記憶部に格納する。そして、格納した複数回の送受信における出力に対して所定の重み付け処理をした後、総和を求める操作を行う。また、受信ビームフォーマ１７の出力から画像（映像）を生成する処理等を行う。バックエンド処理部１８は、表示部１９へ画像（映像）信号を出力し、実時間で断層像や立体走査結果等を表示させる。

　本発明は、受信ビームフォーマ１７の構成と動作に特徴がある。本実施の形態の受信ビームフォーマ１７の構成について図２を用いてより詳しく説明する。

　受信ビームフォーマ１７は、受信波形記憶部２０と、遅延整相部２１と、適応処理部２２、加算器２３とを備えて構成される。適応処理部２２は、複数（ここではＰ個）の部分帯域適応処理部（ＡＤＦ_１～ＡＤＦ_Ｐ）２２－１，２２－２・・・２２－Ｐを有する。複数の部分帯域適応処理部２２－１～２２－Ｐは、周波数帯域で分割した受信信号についてそれぞれ適応重みを求め、受信信号の振幅を適応重みで重み付けする処理を行う。

　受信波形記憶部２０は、受信回路１６が出力するＮチャンネルの受信信号群を、所定の時間長に渡って記憶しながら、逐次、時間的に古い受信信号群を新しい受信信号群に更新する。遅延整相部２１は、受信波形記憶部２０から読み出した受信信号群の波形をチャンネル毎に遅延させる。チャンネル毎の遅延量は、電気音響変換器アレイ１０ａに到達する際の受信信号群の波面（凹面）を平面に補正するように定められている。具体的には、電気音響変換器アレイ１０ａを構成する個々の受信素子と像生成したい画素点（受信焦点）との間の、送信から受信までの超音波パルスの飛行経過時間を幾何的に求めた情報を元に定められる。この波形遅延処理は、受信波形記憶部２０に記憶された受信信号群の時間標本化周期より細かい時間精度で行う。

　遅延整相部２１は、遅延処理した受信信号群（第一の時間信号群）を適応処理部２２の複数の部分帯域適応処理部２２－１～２２－Ｐにそれぞれ出力する。部分帯域適応処理部２２－１～２２－Ｐはそれぞれ、受信信号群（第一の時間信号群）から所定の周波数帯域（部分帯域）の第二の時間信号群を抽出し、抽出した第二の時間信号群について、適応処理を行うことで適応重みを求める。この適応重みにより、受信信号群（第二の受信信号群）の振幅を重み付け（適応重み付与）処理する。部分帯域適応処理部２２－１～２２－Ｐの出力は、加算器２３にて総和される。

　部分帯域適応処理部２２－１～２２－Ｐが適応処理する部分帯域は、互いに一部または全部が重なるように設定してもよい。

　部分帯域適応処理部２２－１～２２－Ｐがそれぞれ処理する部分帯域は、制御部１２が設定する。制御部１２が設定する部分帯域は、予め定められた帯域でもよいし、操作者が入力部１１を介して設定した帯域でもよい。

　部分帯域適応処理部２２－１～２２－Ｐが適応処理する受信信号群のチャンネル数は、それぞれ全チャンネル（Ｎチャンネル）でもよいし、１チャンネル以上の所定の一部のチャンネルでもよい。図２では遅延整相部２１が全チャンネルの受信信号群（第一の時間信号群）を部分帯域適応処理部２２－１～２２－Ｐに出力し、部分帯域適応処理部２２－１～２２－Ｐが設定されているチャンネルをそれぞれ選択して第二の時間信号群を得て、適応処理する例を示しているが、遅延整相部２１がチャンネルを選択してそれぞれの部分帯域適応処理部２２－１～２２－Ｐに出力する構成にしてもよい。

　また、部分帯域適応処理部２２－１～２２－Ｐの適応処理は、同一の条件で行ってもよいし、別々の条件で行ってもよい。例えば、周波数帯域により異なる指向性を持たせた送信パルスを同時送受信（周波数多重送受信）し、得られた受信信号群を部分帯域適応処理部２２－１～２２－Ｐで周波数帯域ごとに処理する場合は、周波数帯域ごとに別の条件で適応処理することが効率が良い。また、受信信号群から生体非線形性や造影剤による調波成分の周波数帯域を他の周波数帯域から分離する場合も、それぞれに別の条件で適応処理を行う自由度を得ることができる。

　さらに、部分帯域適応処理部２２－１～２２－Ｐの適応処理は、互いに独立させて行わせることもできるし、相互に関連させて処理することもできる。例えば、Ｐ個の部分帯域の適応処理に必要なアンサンブル平均の標本数や相関行列の次数などの条件パラメタを選択した一部の帯域だけで直接指定し、その他の部分帯域のパラメタ値を、帯域の中心周波数や時間の単調変化する関数として直線や曲線で内挿する方法などがある。

　次に、部分帯域適応処理部２２－１～２２－Ｐの構成を図３を用いて説明する。図３は、1つの部分帯域適応処理部２２－１の構成を示している。他の部分帯域適応処理部２２－２～２２－Ｐの構成も、図３の部分帯域適応処理部２２－１と同じである。部分帯域適応処理部２２－１は、解析フィルタ３１、部分帯域適応フィルタ３２、合成フィルタ３３とを備えている。

　遅延整相部２１の出力するＮチャンネルの受信信号群（第一の時間信号群）は、解析フィルタ３１に入力される。解析フィルタ３１は、制御部１２の制御下で、所定の部分帯域を取り出すフィルタリング処理（帯域制限処理）を行う。これにより、図２の遅延整相部２１から受け取った受信信号群（第一の時間信号群）から、部分帯域適応処理部２２－１において適応処理すべき部分帯域の第二の時間信号群（部分帯域信号３４）を生成して、部分帯域適応フィルタ３２に出力する。

　このとき、解析フィルタ３１は、フィルタ処理の積和演算に同時に関係すべき複数チャンネル（チャンネル数ｎ、ｎ≦Ｎ）を全チャネルから選択する処理や、受け取った受信信号群の標本化周期とは異なる標本化周期で第二の時間信号群を再標本化することも可能である。電気音響変換器アレイ１０ａの複数の素子は、空間軸方向に物理的に分割されている。各素子の中央位置で代表される空間位置の系列が定められ、空間軸方向に標本化されたチャネル群を形成している。実際に素子が存在するＮチャネルの空間位置の系列より、空間周期（間隔）の大きな位置系列を作り出し、Ｍ（Ｍ≦Ｎ）チャネルにチャンネル数の低減を行う処理は、空間軸方向の再標本化である。また、アナログ－ディジタル変換時の標本化周期の時刻系列より長い周期の時刻系列で出力する処理は、時間軸方向の再標本化である。チャンネル数を低減する空間軸方向再標本化や、単位時間当たりの標本数を低減する時間軸方向再標本化を行うことにより、部分帯域適応フィルタ３２の相関行列あるいは共分散行列の次数を低減することができる。また、時間軸方向の再標本化を行うことにより時間方向のアンサンブル平均標本点の間の独立性が高まる。逆行列の直接計算やＱＲ分解を元にいわゆる逆行列の定理を計算過程に用いる適応フィルタでは、行列の次数の増加に対して処理量が３乗あるいは２乗の規模で増すため、予め空間軸方向および／または時間軸方向再標本化により次数を低減しておくことにより、解析フィルタ３１の再標本化による処理量増加を見込んでも、なお全体として処理量を低減することができる場合がある。

　部分帯域適応フィルタ３２は、第二の時間信号群（部分帯域信号３４）を適応フィルタ処理することにより適応重みを求め、重み付け処理する。重み付け処理後の部分帯域信号（部分帯域適応信号３５）を出力する。部分帯域適応フィルタ３２の具体的な構成については後述する。合成フィルタ３３は、部分帯域適応信号３５を、制御部１２の図示しない指令に基づき、他の部分帯域適応処理部２２－２～２２－Ｐの出力する部分帯域適応信号３５と共通の時刻系列において加算可能な信号に変換する処理を行う。

　解析フィルタ３１は、帯域通過フィルタ処理、直交復調処理と低域通過フィルタ処理との組み合わせ、および、デシメートやダウンコンバート等のデータレート変換処理、のうちの１以上を制御部１２の制御下で行なう。合成フィルタ３３は、解析フィルタ３１が直交復調処理を行った場合には、直交変調処理を行い、解析フィルタ３１がデシメートやダウンコンバート等のデータレート変換処理を行った場合には、オーバーサンプルやアップコンバート等のデータレート復元を行う。

　また、合成フィルタ３３は、部分帯域適応フィルタ３２の出力する部分帯域適応信号３５に対して、図示しない制御部１２の指令に基づき所定の重み付けを必要に応じて行なうことも可能である。

　非常に簡易な構成例としては、解析フィルタ３１を受信信号群から所定の部分帯域のみを通過させる帯域通過フィルタとし、合成フィルタ３３を、インパルス応答が時間軸方向１点の重み付け乗算器とする構成にできる。また、別の例としては、ベースバンド処理を実現するために、解析フィルタ３１は、参照信号を部分帯域の中心周波数とする直交復調器、デシメーションフィルタおよびダウンサンプラで構成することができる。この場合、合成フィルタ３３は、解析フィルタ３１のダウンサンプルに対応したオーバーサンプルを行うオーバーサンプラ、平滑フィルタ、および、第一の時間信号群の全帯域に従って定める基準周波数と部分周波数帯域の基準周波数との差の周波数を参照周波数とする直交変調器等で構成し、アップコンバートを適切に行なうようにする。また、部分周波数帯域のビーム出力をインコヒーレント合成する場合には、合成フィルタ３３は直交検波による包絡線検出を行っても良い。

　合成フィルタ３３の出力が検波されていなければコヒーレント加算出力として図２の加算器２３で適応ビーム出力が得られる。合成フィルタ３３の構成により検波されていれば加算器２３で周波数に関してインコヒーレント加算された適応ビーム出力が得られる。合成フィルタ３３が検波する場合はバックエンド処理部１８での検波処理は省かれる。

　このように、本実施形態によれば、受信信号群の周波数成分が数オクターブにわたる広帯域信号である場合にも、所定の部分帯域ごとに分割して、別々の部分帯域適応処理部２２－１～２２－Ｐで処理することができる。これにより、全ての周波数成分を包含して一つの適応フィルタで処理する場合よりも、生体伝搬不均一の周波数依存性等の影響を低減して、適切に重み付けした受信信号群を得ることができる。

　次に、部分帯域適応フィルタ３２の詳しい構成について図４を用いて説明する。部分帯域適応フィルタ３２としては、公知の適応フィルタを利用した構成にすることができる。本発明の部分帯域適応フィルタ３２は、相関行列あるいは共分散行列の推定に対して、その次数や空間軸方向および／または時間軸方向のアンサンブル平均処理の標本数、それらの更新条件や更新頻度を、直接あるいは間接に操作することができる構成であればどのような適応フィルタであってもよい。例えば、シスリックアレイ等のアルゴリズムと専用並列計算ハードウエアで部分帯域適応フィルタ３２を構成することも可能であるし、信号処理プロセッサなどを利用したソフトウェアで部分帯域適応フィルタ３２の全体を構成することもできる。本実施形態では、部分帯域適応フィルタ３２が、ＱＲ分解（ＱＲ　Ｄｅｃｏｍｐｏｓｉｔｉｏｎ：ＱＲＤ）アルゴリズムを用いるＭＶＤＲ適応フィルタ（ＱＲＤ－ＭＶＤＲ）である例と、ＱＲ分解を経ずに逆行列計算を行うＳａｍｐｌｅ　Ｍａｔｒｉｘ　Ｉｎｖｅｒｓｉｏｎ（ＳＭＩ）法等でのＭＶＤＲ適応フィルタ、ＡＰＥＳ適応フィルタである例とを以下説明する。

　まず、部分帯域適応処理部３２が、ＱＲＤ－ＭＶＤＲアルゴリズムを用いる場合について、図４および図５を用いて説明する。ＱＲＤ－ＭＶＤＲは、相関行列の逆行列演算処理を直接行わずに、目標値（たとえば制約条件の無いＭＶＤＲではゼロ）に対する出力残差のノルムが相似変換やユニタリ変換に対して不変であることを利用する処理である。

　図４のように、ＱＲＤ－ＭＶＤＲアルゴリズムの部分帯域適応フィルタ３２は、適応制御部４１と、部分帯域波形記憶部４２と、シストリックアレイ多重化器４３と、適応フィルタ処理演算を行うシストリックアレイ４５と、係数テーブル記憶部４４と、部分帯域適応出力記憶部４６とを備えて構成されている。シストリックアレイ４５は「セル」と呼ばれる同期式の処理手段から構成される。

　部分帯域波形記憶部４２は、図３の解析フィルタ３１の出力するチャンネル数Ｍの部分帯域信号３４を一定の時間の長さにわたり入力順に保持する。適応制御部４１は、部分帯域波形読み出しアドレス信号４１－１を部分帯域波形記憶部４２に出力する。これにより、部分帯域波形記憶部４２からチャンネル数Ｍの部分帯域信号３４の波形データが読み出され、シストリックアレイ多重化器４３に受け渡される。シストリックアレイ多重化器４３は、図５に示すように、部分帯域信号３４のデータをチャンネル毎に必要時間だけ遅延させながら分配して、シストリックアレイ３５を同期動作させるのに適切な入力データシーケンス５５に並べ替える。シストリックアレイ多重化器３３は、部分帯域信号３４の他に参照信号ベクトルも含めて入力データシーケンス５５を生成する。シストリックアレイ多重化器４３のデータパスの選択制御は、適応制御部４１の出力する多重化アドレス信号４１－２により指定される。また、シストリックアレイ４５に特有の動作モード指令（Ｍ_ＳＹＳ）４１－６も適応制御部４１により出力される。

　並行して適応制御部４１は、係数テーブル記憶部４４は係数多重化アドレス信号４１－３を出力する。係数テーブル記憶部４４は、適応処理前の重みづけベクトルの要素をチャネルごと乗算するための重みづけベクトル係数ｗ_ｉｎｂおよびｗ_ｉｎｒを、時間的な並べ替えを行いながらシストリックアレイ３５の一部の「セル」（図５のアレイ入力側セル群５６）に対して出力する。シストリックアレイ４５では、部分帯域信号３４のチャネルに対し重みづけベクトルによる重みづけが、係数テーブル記憶部４４から受け取った係数との積で行われ、引き続くＱＲ分解と最適重みの賦与が行われる。これにより、適応処理前にチャネル重みをつけた部分帯域の適応処理が実行される。

　シストリックアレイ４５の適応出力は、部分帯域適応出力記憶部４６に格納される。部分帯域適応出力記憶部４６は、シストリックアレイ４５の適応出力を、後段の処理に同期して出力するため所定の一定時間にわたって保持する。適応制御部４１は、部分帯域適応出力記憶部４６に適応出力整序アドレス信号４１－４）により時間的な整序を行って部分帯域適応信号３５を出力する。

　適応制御部４１の出力する部分帯域波形読み出しアドレス信号４１－１、係数多重化アドレス信号４１－３、適応出力整序アドレス信号４１－４は、適応制御部４１に内蔵されているシーケンステーブル記憶部４１－５のテーブルに従って演算され発生される。適応制御部４１のシーケンステーブル記憶部４１－５および係数テーブル記憶部３４に記憶されるテーブルの内容は、図１の制御部１２からデータバス１２－１を経由して図４の適応制御部４１および係数テーブル記憶部４４へと転送される。

　ここで図４のシストリックアレイ４５の構成についてＭ＝４の場合の例に図５を用いて説明する。シスリックアレイ４５の構成は、非特許文献３を含め多くの構成法があるので、本発明に関係する部分だけを簡単に説明する。シストリックアレイ４５は、境界セル５１、遅延セル５２、内部セル５３、最終セル５４の少なくとも４種の「セル」を含んで構成される。「セル」は所定の並列信号処理の繰り返し単位の纏まりを指し、個々の「セル」は物理的なハードウェア演算手段の分割単位を必ずしも意味しない。図５のシストリックアレイ４５には、小円の図形で表した最終セル５４は１個、比較的大きな円の図形で表した境界セル５１と同じものは４個、傾斜した長方形の図形で表した遅延セル５２と同じものは４個、正方形の図形で表した内部セル５３と同じものは１０個含まれている。

　これらの「セル」が演算回路の単位である必要はなく、たとえば境界セル５１、遅延セル５２、内部セル５３、最終セル５４の中から個数の異なる複数種の「セル」の部分集合をまとめて特定の物理的な演算手段の処理とし、それらを並列に配置して構成しても良い。

　図４のシストリックアレイ多重化器４３の出力する入力データシーケンス５５は、図５のように、部分帯域信号３４の同一時刻である入力データベクタ群｛ｘ_１１、ｘ_１２、ｘ_１３、ｘ_１４｝、｛ｘ_２１、ｘ_２２、ｘ_２３、ｘ_２４｝、｛ｘ_３１、ｘ_３２、ｘ_３３、ｘ_３４｝…、および参照信号ベクタ｛ｙ_１、ｙ_２、ｙ_３、…｝をそれぞれ所定の遅延量で整序したものである。入力データシーケンス５５は、シストリックアレイ４５の一辺に並ぶ境界セル５１や内部セル５３に対してそれぞれ入力される。図５において、下向き矢印の連続は、シストリックアレイ４５の全ての「セル」群が演算結果を隣接したセルに同期転送を行う時間を単位として入力遅延させることを表す。たとえば、シストリックアレイ４５へのｘ_１４の入力はｘ_１１に比べて３回の同期周期分遅延させて入力されることを表している。シストリックアレイ４５の三角形の上辺に沿う方向（チャネル番号方向）の入力ベクタ要素の分配、すなわち境界セル５１とその右に接続される内部セル５３を含む４つの内部セルのどれに対して入力データベクタの要素を転送するかは、図４のシストリックアレイ多重化器４３が適応制御部４１からの多重化アドレス信号４１－２により指定されることにより決定される。また、図５の入力データシーケンス５５に対応して入力ベクタ要素を配置するかは、図４の部分帯域波形記憶部４２が適応制御部４１から部分帯域波形読み出しアドレス信号４１－１により指定される。

　アレイ入力側セル群５６に含まれるセルは適応処理を行う前の重み付けベクトルの要素と入力データベクタ群の要素との要素ごとの積を行う。これらは入力データシーケンス５５に合わせて図４の適応制御部４１が指定する係数多重化アドレス信号４１－３により係数テーブル記憶部４４が図５に図示しないバスを経由してアレイ入力側セル群５６に含まれるそれぞれのセルの係数ｗ_ｉｎｂ，ｗ_ｉｎｒとなる。

　図５の境界セル５１、内部セル５３および最終セル５４の演算内容を図６に示す。境界セル５１、内部セル５３はアレイ入力側セル群５６に含まれない場合は係数ｗ_ｉｎｂおよびｗ_ｉｎｒは１となる。遅延セル５２は、一回の同期分だけ結果を遅延させて転送するセルである。

　図６に示すように、境界セル５１、内部セル５３および最終セル５４の演算内容は、図４の適応制御部４１の出力する動作モード指令（Ｍ_ＳＹＳ）４１－６によって交番する。図６において、ｘ_ｉｎは_、セルへの要素入力，ｘ_ｏｕｔは、セルからの要素回転出力，ｃは、Ｇｉｖｅｎｓ回転余弦、ｓは、Ｇｉｖｅｎｓ回転正弦、ｒは、上三角行列要素，γ_ｉｎは、変換係数入力，γ_ｏｕｔは、変換係数出力，βは、忘却係数（指数重み係数），μは、拘束条件値、ｗ_ｉｎｂおよびｗ_ｉｎｒは、重みづけベクトル係数である。重みづけベクトル係数ｗ_ｉｎｂ，ｗ_ｉｎｒ以外は、ＱＲＤ－ＭＶＤＲアルゴリズムの係数として広く知られており、具体的には、非特許文献３等に広く開示されているのでここでは説明を省略する。

　重みづけベクトル係数ｗ_ｉｎｂ，ｗ_ｉｎｒは、セルへの要素入力ｘ_ｉｎに乗じられる係数である。例えば、重みづけベクトル係数ｗ_ｉｎｂ，ｗ_ｉｎｒの値を、０または１とすれば、入力データベクトルの次数を増減したり、受信口径位置を設定したのと同様の効果が得られる。また、受信開始からの経過時刻に対して滑らかに増減すれば、受信口径位置を滑らかに移動させたのと同様の効果を得ることができる。

　次に、部分帯域適応フィルタ３２が、逆行列計算を行うＳＭＩ法によるＭＶＤＲ適応フィルタである例について説明する。この部分帯域適応フィルタ３２は、図７に示すように、図４の構成においてシストリックアレイ多重化器４３およびシストリックアレイ４５の代わりに、演算部７０を配置する。演算部７０は、単数あるいは複数の数値演算プロセッサで構成しても良く、以下に説明する数式を用いて演算を行う。部分帯域波形記憶部４２、係数テーブル記憶部４４、部分帯域適応出力記憶部４６は数値演算プロセッサに付随する主記憶装置で構成しても良い。係数テーブル記憶部４４に記憶されるデータは適応前にデータに重みづけするベクトル係数ｗ_ｉｎｂ，ｗ_ｉｎｒではなく、後述する可干渉性損失行列の要素となる。それ以外の構成は、図４の部分帯域適応処理部３２と同様であるので説明を省略する。

　図７の部分帯域波形記憶部４２には、Ｍチャンネル（Ｍ≦Ｎ）の部分帯域波形が保持されている。部分帯域波形記憶部４２は、適応制御部４１からの部分帯域波形読み出し読み出しアドレス信号４１－１により指定された、所定の標本化時刻ｔにおける所定の複数チャンネルのデータを演算部に受け渡す。全Ｍチャンネルの時刻ｔのデータは、式（１）のようにベクトルｘ（ｔ）として表すことができる。
　　　ｘ（ｔ）＝［ｘ_ｔ１，ｘ_ｔ２，…，ｘ_ｔＭ］^Ｔ　　　・・・・（１）

　式（１）において右上付き記号「^Ｔ」は転置を表し、ｘ（ｔ）は列ベクトルである。以下、本実施例においてベクトルは特に指定しない限り全て列ベクトルとする。ｘ（ｔ）は、図３の解析フィルタ３１の構成により、実数ベクトルあるいは、直交変調あるいは復調を受けた複素数ベクトルであってもよい。このベクトルはある瞬間時刻ｔでの電気音響変換器アレイ１０ａの出力データ標本値であり、説明の便を図るため以下「スナップショット」と呼称する。演算部３８では、このｘ（ｔ）を重み付けする適応重みを、スナップショットｘ（ｔ）の時間軸方向近傍のデータから推定する。
　　　ｗ（ｔ）＝［ｗ_ｔ１，ｗ_ｔ２，…，ｗ_ｔＭ］^Ｔ　　　・・・・（２）

　部分帯域適応フィルタ３２（以下、ＭＶＤＲ適応フィルタとも呼ぶ）の演算部７０の出力ｂｍ（ｔ）は、式（３）で表される。
　　　ｂｍ（ｔ）＝ｗ^Ｈ（ｔ）ｘ（ｔ）　　　・・・・（３）
ここで右上付き記号「^Ｈ」は行列やベクトルの共役転置（エルミート共役）を表わし、ｘが実数ベクトルである場合も虚数成分が０の複素数ベクトルとして解釈する。

　ＭＶＤＲアルゴリズムに必要な相関行列（空間共分散行列）Ｒ（ｔ）は「スナップショット」ベクトルｘ（ｔ）のダイアド（ｄｙａｄ）積ｘｘ^Ｈが形成する行列の期待値で決定される。期待値演算をＥ｛　｝で表すとすれば、相関行列（空間共分散行列）Ｒ（ｔ）は、式（４）で表わされる。
　　　Ｒ（ｔ）＝　Ｅ｛ｘｘ^Ｈ｝　　　・・・・（４）

　受信焦点の信号源とアレイの間の位置関係がもたらす位相差を補償する、ＭＶＤＲにおけるモードベクトル（あるいはステアリングベクトル）ａ_ｍは、適切な複素数Θ_１，Θ_２，…，Θ_Ｎと指数関数ｅｘｐ（）を用いて式（５）で表わされる。
　ａ_ｍ＝［ｅｘｐ(Θ_１)，ｅｘｐ(Θ_２)，…，ｅｘｐ(Θ_Ｍ)］^Ｔ・・・（５）
本実施例では、ｘ（ｔ）は、到来方向に合わせた遅延後の信号であるため、推定到来方向に合わせて位相で偏向する必要が無いため、Θ_１，Θ_２，…，Θ_Ｍは全てゼロとする。

　最適重みに対する拘束条件を式（６）とし、ＭＶＤＲの定義からコスト関数ε＝ｗ^Ｈ（ｔ）ｘ（ｔ）を最小化するｗ（ｔ）をＬａｇｒａｎｇｅ未定係数法で求めると、ｗ（ｔ）の最適解は、Ｒ（ｔ）の逆行列をＲ^－１（ｔ）として式（７）で表わされる。

　　　ｗ^Ｔ（ｔ）ａ_ｍ＝１　　　・・・・（６）
　　　ｗ（ｔ）＝Ｒ^－１（ｔ）ａ_ｍ／｛ａ_ｍ ^ＨＲ^－１（ｔ）ａ_ｍ｝　・・・・（７）

　一方、有限の標本数の期待値演算からは、相関行列（空間共分散行列）Ｒ（ｔ）に対応する相関行列（空間共分散行列）推定量^∧Ｒ（ｔ）が得られる。式（４）の相関行列Ｒ（ｔ）の推定量^∧Ｒ（ｔ）を得るために、適応重みｗ（ｔ）を推定したい時刻ｔの近傍で、スナップショットベクトルｘ（ｔ）について時間軸方向の有限範囲でアンサンブル平均を取り、元のベクトルｘより次数が少ない部分ベクトルｘ_Ｓを選び出し、それらのダイアド積を空間軸方向にもアンサンブル平均を行ってＲ（ｔ）を近似する推定量^∧Ｒ_Ｓ（ｔ）を式（８）で求める。記号「^∧」は推定量を表す。

　式（８）において、ｘ_Ｓ(t,ａ)は、ベクトルｘ（ｔ）において、先頭位置をa（a＝１，２，…，Ａ）とする連続Ｋチャンネル（Ｋ≦Ｍ≦Ｎ）の部分ベクトルである。時間平均化の幅Ｌは奇数であるが、偶数として―Ｌ／２≦ｌ≦Ｌ／２－１とするなどしてよい。元のスナップショットから空間軸方向についてアンサンブル平均を計算するために、得られる適応重みは次数を低下させて近似する。例えば、Ｍチャンネルの中から小さな連続Ｋチャンネルの部分ベクトルを取り出すとＡ＝Ｍ―Ｋ＋１個の空間アンサンブル平均が得られる。十分に確率収束して安定した推定量の相関行列^∧Ｒ_Ｓ（ｔ）を得るためには必要十分な数の独立したアンサンブル標本集合が必要である。

　式（５）のモードベクトルaの次元をＭからＫに縮小したモードベクトルａ_Sは、式（９）で表わされる。最適重み推定量^∧ｗ_Ｓ（ｔ）は、式（７）と同様の式（１０）で表わされる。

　　ａ_Ｓ＝［ｅｘｐ(Θ_１)，ｅｘｐ(Θ_２)，…，ｅｘｐ(Θ_Ｋ)］^Ｔ・・・（９）
　^∧ｗ_Ｓ(ｔ)＝^∧Ｒ_Ｓ ^－１(ｔ)^∧ａ_Ｓ／｛^∧ａ_Ｓ ^Ｈ　∧Ｒ_Ｓ ^－１(ｔ)^∧ａ_Ｓ｝　・・・（１０）

　ｘ（ｔ）をｘ_Ｓ(t,ａ)から近似する部分ベクトル推定量^∧ｘ_Ｓ（ｔ）は、式（１１）で表わされる。

　部分帯域適応フィルタ３２の演算部７０は、式（１２）により、適応重み推定量^∧ｗ_Ｓ（ｔ）で振幅の重み付け後の出力^∧ｂｍ_Ｓ（ｔ）を求めて出力する。
　^∧ｂｍ_Ｓ（ｔ）＝　^∧ｗ_Ｓ ^Ｈ（ｔ）^∧ｘ_Ｓ（ｔ）　　　　・・・・（１２）

また、ＭＶＤＲ法を改良したＡＰＥＳ法を用いて演算部７０が演算することも可能である。ＡＰＥＳ法では、式（５）のモードベクトルのa番目要素ｅｘｐ(Θa)を用いて、式（１３）の推定ベクトル^∧ｇ（ｔ）を計算する。

　求めた^∧ｇ（ｔ）を用いて、式（８）の推定量の相関行列^∧Ｒ_Ｓ（ｔ）を修正した相関行列^∧Ｒ_ＡＰＥＳ（ｔ）を求める。

　^∧Ｒ_ＡＰＥＳ（ｔ）＝　^∧Ｒ_Ｓ（ｔ）－^∧ｇ（ｔ）^∧ｇ^Ｈ（ｔ）　・・・・（１４）

　そして、式（１０）の^∧Ｒ_Ｓ（ｔ）を^∧Ｒ_ＡＰＥＳ（ｔ）に置き換えて、最適重み^∧ｗ_Ｓ(ｔ)を求める。他は、上記ＭＶＤＲ法の各式と同様にして、演算部７０により、適応重み^∧ｗ（ｔ）で振幅の重み付け後の出力^∧ｂｍ（ｔ）を求め、出力する。

　部分帯域適応処理部２２－１～２２－Ｐにおいて、解析フィルタ３１が、受信信号群を時間軸方向および／または空間軸方向について再標本化処理を行うことで、部分帯域適応フィルタ３２での処理次数を低減できる。この再標本化処理について以下詳しく説明する。

　複数の部分帯域適応処理部２２－１～２２－Ｐにそれぞれ設定される部分帯域は、予め定めた帯域または操作者が設定した帯域に設定可能である。例えば、解析フィルタ３１による時間軸方向の再標本化は、図２の遅延整相部２１の出力する遅延整相された第一の時間信号群（受信波形信号）に対し、部分帯域適応処理部２２－１～２２－Ｐの各部分帯域に応じて時間軸方向に所定の再標本化周期の第二の時間信号群が得られるように行う。再標本化出力の計算のために、第一の時間信号群の時間軸方向にｍ点の幅が用いられる。ｍは部分帯域の周波数に応じて必要な値を設定する。一例としては、帯域を代表する中心周波数が高い部分帯域では、小さなｍで時間軸方向に短い周期で再標本化出力し、中心周波数が低い部分帯域では、大きなｍで長い周期で再標本化出力するよう設定する。

　一方、空間軸の再標本化は、以下のように行う。図１の超音波探触子１０の内蔵する電気音響変換器アレイ１０ａの個々の電気音響変換器は、通常、距離方向または角度方向に一定の間隔で配列されている。この間隔は、送信、受信を行う超音波周波数に応じて決定されている。例として超音波周波数が２～１８ＭＨｚまで扱う場合を考えると、上限と下限の周波数の間で９倍の違いがある。上限周波数に対応して電気音響変換器アレイ１０ａの間隔を設定すると、下限の周波数にとっては９倍の空間標本化周波数で必要以上に細かく標本化されていることになる。そこで、部分帯域適応処理部２２－１～２２－Ｐのそれぞれの解析フィルタ３１は、設定されている部分帯域の周波数に応じた空間軸方向（電気音響変換器アレイ１０ａの分割方向、受信チャンネルの番号方向）での再標本化を行う。遅延整相部２１の出力する遅延整相された受信信号群である第一の時間信号群に対し、たとえば配列のチャンネル番号方向に部分帯域ごとに必要な所定の長さｎチャンネルにわたる空間軸方向の幅を用いる必要がある。ｎは部分帯域の周波数に応じて必要な値を設定する。一例としては、帯域を代表する中心周波数が高い部分帯域では、小さなｎで空間軸方向に短い周期で再標本化出力し、中心周波数が低い部分帯域では、大きなｎで長い周期で再標本化出力するよう設定する。

　上記時間軸方向および空間軸方向の再標本化を同時に実現するために、例えば解析フィルタ３１の演算は、遅延整相部２１の出力する遅延整相された受信波形信号である第一の時間信号群に対し、時間軸方向にｍ点、空間軸方向にｎ点の行列Ｆの要素との積和を計算する構成（２次元標本化関数、あるいは２次元補間関数）とする。これにより、時間と空間の再標本化および部分帯域への制限を同時に行うことができる。
　Ｆ＝｛ｆ_ｕv｝，ｕ＝１，２，…，ｍ；ｖ＝１，２，…，ｎ　・・・（１５）
式（１５）において、ｆ_ｕvは、行列Ｆのｕ行ｖ列の要素を示す。

　行列Ｆの各要素は、遅延整相部２１の出力する遅延整相された受信波形信号の標本化時刻系列と、物理的な受信素子による空間の標本化位置系列とにより定まる２次元格子点間の位置についてダウンサンプルやデシメートが行われた再標本化出力が得られる。出力したい点の近傍２次元格子点群からの時間軸方向、空間軸方向ずれ量に応じて、異なる行列要素値ｆ_ｕvを用いる。例えば時間軸方向に遅延整相部２１の出力する遅延整相された受信波形信号の時間周期の３回につき２回の出力で、かつ、チャンネルによる空間標本化の４チャンネル幅につき３チャンネルの出力といった２次元の再標本化あるいはレート変換を実現することができる。

　解析フィルタ３１は、それぞれの部分帯域適応処理部２２－１～２２－Ｐに設定されている部分帯域の周波数に応じて再標本化して、第二の時間信号群を得る。部分帯域適応処理部２２は、独立した適応処理で重み付けした第二の時間信号群を加算出力する。部分帯域適応処理部２２の出力時刻系列は、部分帯域適応処理部２２－１～２２－Ｐごとに異なる。これらＰ個の独立時刻系列出力を加算器２３で時間軸上の同一時刻系列で加算するために、合成フィルタ３３は、部分帯域適応フィルタ３２の出力時刻系列群に対して時刻系列の変換（レート変換）を行う。これにより、加算器２３は、部分帯域ごとに分割して部分帯域適応処理部２２－１～２２－Ｐで適応処理された時間信号群（部分帯域適応信号２５）を同一の時刻系列で加算することができる。

　上記のように解析フィルタ３１と合成フィルタ３を備えることにより、部分帯域適応処理部２２－１～２２－Ｐに設定されている部分帯域のそれぞれで、相関行列Ｒ（ｔ）の推定において時間軸方向のアンサンブル平均や逆行列計算を行う単位時間当たりの頻度を最小にする時刻系列が部分帯域適応処理部２２－１～２２－Ｐごとに異なり、結果的に適応出力の時刻系列が異なる場合でも、その出力を加算器２３により加算することができる。

　解析フィルタ３１による時間軸方向および空間軸方向の再標本化が部分帯域適応フィルタ３２へ与える効果を説明する。所定の時刻ｔでのスナップショット（式（１）のベクトル）の近傍時間でアンサンブル平均を十分な標本数で行えば、相関行列の推定量^∧Ｒ（ｔ）が良好に確率収束することが期待できる。

　部分帯域の周波数の周期に対して時間軸方向に過度に短い周期で標本化されているとデータの独立性が悪く平均化演算の効果が期待できなくなる。スナップショット（式（１）のベクトル）の次数を全チャネル数Ｎより少なくすると、同じ時刻で空間軸方向にも平均化が行える上に、逆行列の演算量も低減する。

　解析フィルタ３１は図２の遅延整相部２１の出力する遅延整相された第一の時間信号群の時刻番号ｔでＮチャンネルのスナップショットｘ（ｔ）を空間軸方向の標本化位置系列の間隔を元の物理的受信素子の間隔Ｓ_ＳのＳ_ｓｕｂ（＞１）倍で再標本化する。新たな空間位置系列では最大でＭ_ｓｕｂ（≦Ｎ／Ｓ_ｓｕｂ）の標本数が同時刻で得られる。この中から連続したＫ_ｓｕｂ（＜Ｍ_ｓｕｂ）チャンネルを選び出す。最終的に物理的なＮチャンネルから再標本化されたＫ_ｓｕｂチャネルのスナップショットがＭ_ｓｕｂ－Ｋ_ｓｕｂ＋１通り得られる。

　図２の遅延整相部２１の出力する遅延整相された受信波形信号のチャンネル番号（電気音響変換器アレイ１０ａの受信素子の番号）をｊとすれば、空間軸方向にＳ_ｓｕｂ倍の周期で再標本化して得た空間位置系列の再標本化チャンネル番号ｊ’との関係は式（１６）により表すことができる。
　ｊ’＝（ｊ－１）／Ｓ_ｓｕｂ＋１　　　　　　　　　　　　・・・（１６）
ここで、ｊ’＝１，２，…,Ｍ_ｓｕｂである。ｊ’の示す空間位置は、元のｊの空間位置系列位置の何れの位置にも一致しない場合がある。

　同様に解析フィルタ３１は図２の遅延整相部２１の出力する遅延整相された標本化周期Ｔ_Ｓの受信波形信号（第一時間信号群）をＴ_ｓｕｂ倍の周期に広げて再標本化した新たな時刻系列値を作り出す。再標本化時刻番号をｔ’とすれば、ｔとｔ’との関係は式（１７）で表される。
　ｔ’＝（ｔ－１）／Ｔ_ｓｕｂ＋１　　　　　　　　　　　　・・・（１７）

　解析フィルタ３１が出力する再標本化時刻番号ｔ’における再標本化チャネル番号ｊ’の信号値をｘ_{ｔ’,ｊ’}とする。再標本化時刻ｔ’＝ｔ_０を原点０とする整数変数ｌ’（値域はＬ’を奇数として－（Ｌ’－１）／２≦ｌ’≦（Ｌ’－１）／２）および再標本化チャンネルｊ’＝１を原点とする整数変数ａ’（値域はＡを自然数として１≦ａ’≦Ａ’）を用いて格子点座標（ｌ’，ａ’）を設定する。（ｔ’，ｊ’）＝（ｔ_０，１）を原点とする相対的な座標（ｌ’,ａ’）で指定する部分アレイベクトルｘ_Ｓｓｕｂ（ｌ’,ａ’）は、第一要素のチャネル番号をｊ’＝ａ’とし、以下順にｊ’が増加する方向に連続Ｋ_ｓｕｂチャネルを要素とする部分アレイベクトルｘ_Ｓｓｕｂ（ｌ’,ａ’）は式（１８）のようになる。
ｘ_Ｓｓｕｂ（ｌ’,ａ’）
＝［ｘ_{ｔ０＋l’ａ’}，ｘ_{ｔ０＋l’（ａ’＋１）}，…，ｘ_{ｔ０＋l’（ａ’＋Ｋｓｕｂ－１）}］^Ｔ
　　　　　　　　　　　　　　　　　　　　　　　　　　・・・（１８）

　つぎに、図８を用いて、解析フィルタ３１が、遅延整相部２１の出力する遅延整相された第一の時間信号群から上記部分アレイベクトルｘ_Ｓｓｕｂ（ｌ’，ａ’）の要素を成す信号値をｘ_{ｔ’,ｊ’}を生成する処理をさらに詳しく説明する。図８において、小円の格子配列は解析フィルタ３１に入力されるチャンネル番号ｊ（ｊ＝１～Ｎ）の各々の時間軸方向の時刻ｔ_０近傍（・・・，ｔ_０－５，・・・，ｔ_０－１，ｔ_０，ｔ_０＋１，・・・，ｔ_０＋５，・・・）の入力値を模式的に示したものある。これらは時刻番号ｔおよびチャンネル番号ｊで指定される格子点で標本化されており、たとえば解析フィルタ入力点８１は時刻番号ｔ＝ｔ_０－５、チャンネル番号ｊ＝２の標本値である。

　同じく、図８において正方形の格子配列は、解析フィルタ３１から出力される再標本化チャンネル番号ｊ’（ｊ’＝１～Ｍ_ｓｕｂ）の時間軸方向の再標本化時刻番号ｔ’＝ｔ_０’近傍（・・・ｔ_０’－３，・・・，ｔ_０’－１，ｔ_０’，ｔ_０’＋１，・・・ｔ_０’＋３，・・・）の出力値を模式的に示したものである。これらは再標本化時刻番号ｔ’および再標本化チャンネル番号ｊ’で指定される格子点で標本化されており、解析フィルタ出力点８２はｔ’＝ｔ_０’－３，ｊ’＝４の標本値である。小円の格子配列に対し、正方形の格子配列は空間軸方向の間隔はＳｓ×Ｓ_ｓｕｂ倍、時間軸方向の間隔がＴｓ×Ｔ_ｓｕｂ倍になっている。

　解析フィルタ３１の出力点８３（ｘ_{ｔ’,ｊ’}）の生成について例示する。出力点８３を得るために、入力の格子点基準位置８４の座標（ｔ_０－２，３）からチャンネル番号ｊの増加方向にｎ点、時刻番号ｔの増加方向にｍ点の標本値群８５と式（１５）の行列Ｆの要素｛Ｆ_ｕｖ｝とを積和する。図８の出力ｘ_{ｔ’,ｊ’}の位置をチャンネル番号ｊ’の増加方向および時刻番号ｔ’の増加方向に移動させながら、同時に行列｛Ｆ_ｕｖ｝の値を適切に変更し、｛Ｆ_ｕｖ｝と標本値群８５との積和を繰り返す。この積和の反復により、時間軸方向および空間軸方向に再標本化された遅延波形標本値を得る。

　部分帯域適応処理部２２－１～２２－Ｐごとに独立した条件で行う解析フィルタ３１の再標本化の処理に加え、部分帯域適応フィルタ３２での部分アレイベクトルの構成も独立した条件で行う。以下具体例を用いて説明する。

　部分帯域適応処理部２２－１において解析フィルタ３１出力から部分帯域適応フィルタ３２がアンサンブル平均のために部分アレイベクトルｘ_{Ｓｓｕｂ１}を構成する方法を図９に示す。同様に、部分帯域適応処理部２２－Ｐにおいて解析フィルタ３１出力から部分帯域適応フィルタ３２がアンサンブル平均のために部分アレイベクトルｘ_{ＳｓｕｂＰ}を構成する方法を図１０に示す。

　図９の正方形格子の縦軸は、解析フィルタ３１が時間軸方向に標本化周期をＴ_ｓｕｂ＝Ｔ_ｓｕｂ１倍した時刻系列の再標本化時刻番号ｔ_１’を、横軸は解析フィルタ３１が空間軸方向に標本化周期をＳ_ｓｕｂ＝Ｓ_ｓｕｂ１倍した空間格子間隔上のて再標本化チャンネル番号ｊ_１'（ｊ_１'＝１，２，…，Ｍ_ｓｕｂ１）を示す。アンサンブル平均は時間軸方向に幅Ｌ_１点，空間軸方向に幅Ａ_１点行う。再標本化時刻ｔ_１’＝ｔ_０１を原点０とする整数変数ｌ_１（値域は－（Ｌ_１－１）／２≦ｌ_１≦（Ｌ_１－１）／２）および再標本化チャンネルｊ_１’＝１を原点とする整数変数ａ_１（値域は１≦ａ_１≦Ａ_１）を用いて格子点座標（ｌ_１，ａ_１）を設定する。図９の格子点９１は（ｌ_１，ａ_１）＝（－１，６）である。アンサンブル平均のための部分アレイベクトル９３（ｘ_{Ｓｓｕｂ１}）の第一要素位置９１を格子点座標（ｌ_１，ａ_１）にある解析フィルタ３１の出力とし、ｊ_１'の増加方向に沿ってＫ_１個の格子点の出力値によりベクトルの要素が割り当てられる。第一要素の格子点座標を用いて部分アレイベクトル９３をｘ_{Ｓｓｕｂ１}（ｌ_１，ａ_１）と表記する。第一要素位置９１の格子点座標（ｌ_１，ａ_１）はアンサンブル平均範囲９２の中で設定される。図９から明らかなように、式（１９）の不等式関係を満たす。

　　Ａ_１＋Ｋ_１－１≦Ｍ_ｓｕｂ１　　　　　　　　　　　　　　・・・（１９）
　アンサンブル平均範囲９２は、ｊ_１'＝６、再標本化時刻番号ｔ_１’＝ｔ_０１を中心とする平均化と考える。部分帯域適応処理部２２－１のｔ_１’＝ｔ_０１における相関行列推定量^∧Ｒ_１（ｔ_０１）は、部分アレイベクトルｘ_{Ｓｓｕｂ１}（ｌ_１，ａ_１）のダイアド積の変数ｌ_１、ａ_１による２重総和から以下の式（２０）により計算する。

　部分ベクトル推定量^∧ｘ_{Ｓｓｕｂ１}（ｔ_０１）は以下の式（２１）で計算する。

　図１０の正方形格子の縦軸は、解析フィルタ３１が時間軸方向に標本化周期をＴ_ｓｕｂ＝Ｔ_ｓｕｂＰ倍した時刻系列の再標本化時刻番号ｔ_Ｐ’を横軸は解析フィルタ３１が空間軸方向に標本化周期をＳ_ｓｕｂ＝Ｓ_ｓｕｂＰ倍した空間格子間隔上のて再標本化チャンネル番号ｊ_Ｐ'（ｊ_Ｐ'＝１，２，…，Ｍ_ｓｕｂＰ）を示す。アンサンブル平均は時間軸方向に幅Ｌ_Ｐ点，空間軸方向に幅Ａ_Ｐ点行う。再標本化時刻ｔ_Ｐ’＝ｔ_０Ｐを原点０とする整数変数ｌ（値域は－Ｌ_Ｐ／２≦ｌ_Ｐ≦Ｌ_Ｐ／２－１）および再標本化チャンネルｊ_Ｐ’＝２を原点とする整数変数ａ_Ｐ（値域は１≦ａ_Ｐ≦Ａ_Ｐ）を用いて格子点座標（ｌ_Ｐ，ａ_Ｐ）を設定する。図１０の格子点１０１は（ｌ_Ｐ，ａ_Ｐ）＝（－１，４）である。アンサンブル平均のための部分アレイベクトル１０３（ｘ_{ＳｓｕｂＰ}）の第一要素位置１０１を格子点座標（ｌ_Ｐ,ａ_Ｐ）にある解析フィルタ３１の出力とし、ｊ_Ｐ'の増加方向に沿ってＫ_Ｐ個の格子点の出力値によりベクトルの要素が割り当てられる。第一要素の格子点座標を用いて部分アレイベクトル１０３をｘ_{ＳｓｕｂＰ}（ｌ_Ｐ，ａ_Ｐ）と表記する。第一要素位置１０１の格子点座標（ｌ_Ｐ，ａ_Ｐ）はアンサンブル平均範囲１０２の中で設定される。図１０から明らかなように、式（２２）の不等式関係を満たす。

　　　　　　Ａ_Ｐ＋Ｋ_Ｐ－１≦Ｍ_ｓｕｂＰ　　　　　　　・・・（２２）
　アンサンブル平均範囲１０２は、ｊ_Ｐ'＝５、再標本化時刻番号ｔ_Ｐ’＝ｔ_０Ｐ－１とｔ_Ｐ’＝ｔ_０Ｐの間を中心とする平均化であると考える。部分帯域適応処理部２２－Ｐのｔ_Ｐ’＝ｔ_０Ｐにおける相関行列推定量^∧Ｒ_Ｐ（ｔ_０Ｐ）は、部分アレイベクトルｘ_{ＳｓｕｂＰ}（ｌ_Ｐ，ａ_Ｐ）のダイアド積の変数ｌ_Ｐ、ａ_Ｐによる２重総和から以下の式（２３）により計算する。

　部分ベクトル推定量^∧ｘ_{ＳｓｕｂＰ}（ｔ_０Ｐ）は以下の式（２４）で計算する。

　図９と図１０に示した部分帯域適応フィルタ３２の処理では、空間軸方向のアンサンブル平均幅Ａ_１とＡ_Ｐ，時間軸方向のアンサンブル平均長さＬ_１とＬ_Ｐ，部分アレイベクトルのスナップショットの長さＫ_１とＫ_Ｐが異なる。また、図９では再標本化チャンネルＭ_ｓｕｂ１の全てを使用したのに対し、図１０では再標本化チャンネルＭ_ｓｕｂＰの限定範囲からスナップショットを作り出している。

　このように、本実施形態では、部分帯域適応処理部２２－１～２２－Ｐごとに、解析フィルタ３１が独立して再標本化を行うことにより、部分帯域適応フィルタ３２が空間、時間のアンサンブル平均化を効率よく自由度高く行える。

　生体不均一のために受信信号の不確定性が増し、部分帯域ごとに必要十分な確率収束を期するアンサンブル標本数を増す必要がある場合、部分アレイベクトルのスナップショットの長さＫ_１やＫ_Ｐを短くして空間のアンサンブル平均幅Ａ_１とＡ_Ｐ，あるいは時間のアンサンブル平均長さＬ_１とＬ_Ｐを大きくすることができる。また、その逆も行える。仮に、従来技術に見られるように単一の広帯域に全ての周波数成分が含まれている場合、高周波側などの生体不均一の影響をより受けやすい成分が相関行列推定値の確率収束の良否を支配することが予想される。これに対し、本実施形態では部分帯域ごとに独立した適応処理を行い、部分帯域ごとに自由度の高い部分アレイベクトルのスナップショットの長さと空間、時間のアンサンブル平均化を行うことができる。生体不均一の影響を受けやすい高周波側ではアンサンブル平均数を大きくし、低周波側では小さくすることにより、広帯域の受信信号を有効に適応処理することができる。

　上述した空間軸方向および時間軸方向の値Ａ_１，Ａ_Ｐ，Ｌ_１，Ｌ_Ｐ，Ｋ_１，Ｋ_Ｐ、およびアンサンブル平均範囲９２、１０２の位置などは、再標本化時刻番号ｔ_１’、ｔ_Ｐ’や再標本化チャンネル番号ｊ_１’、ｊ_Ｐ’の単位で変化する。時刻番号やチャンネル番号は格子点で量子化されている。このため、これらＡ_１，Ａ_Ｐ，Ｌ_１，Ｌ_Ｐ，Ｋ_１，Ｋ_Ｐ、およびアンサンブル平均範囲９２、１０２の値が、比較的小さな自然数である時にこれらを増減（すなわち設定値を変更）すれば、増減変化が引き起こす信号値の変化、確率収束の度合いや最終の重み推定値の変化分は大きく不連続変化をもたらすことがある。再標本化時刻番号ｔ_１’、ｔ_Ｐ’が小さい受信開始からあまり時間経過していない受信信号群は、超音波探触子１０の近傍部分に位置する被検体内の反射点１０ｃからのパルスエコー１０ｄの信号であり、電気音響変換器アレイ１０ａの各素子の持つ有効な指向性角度の範囲が限定されるため、ビームフォームに関与できる電気音響変換器アレイ１０ａから再標本化された受信チャネルの範囲位置や幅は変化する。これに対応させて上記Ａ_１，Ａ_Ｐ，Ｌ_１，Ｌ_Ｐ，Ｋ_１，Ｋ_Ｐの設定値やアンサンブル平均範囲９２、１０２の変更により、適応処理後の値が不連続に大きく変化する場合がある。

　そこで、適応処理における急激な確率収束や最終の重み推定値の変化を生じないようにするためダイアド積ｘ_{Ｓｓｕｂ１}ｘ_{Ｓｓｕｂ１} ^Ｈやｘ_{ＳｓｕｂＰ}ｘ_{ＳｓｕｂＰ} ^Ｈあるいはそれらの平均後の推定量の相関行列^∧Ｒ_１や^∧Ｒ_Ｐの行列要素を可干渉性損失行列によって重み付けする。

　部分帯域適応処理部２２－１～２２－Ｐの個々に対応する帯域を部分帯域番号ｉ（ｉ＝１，２，…，Ｐ）で識別する。可干渉性損失行列はその要素が部分帯域ごとに時間変化するのでＢ_ｉ（ｔ_ｉ’）と記載する。可干渉性損失行列Ｂ_ｉ（ｔ_ｉ’）を用いて、部分帯域番号ｉの部分アレイベクトルのダイアド積ｘ_{Ｓｓｕｂｉ}ｘ_{Ｓｓｕｂｉ} ^Ｈや相関行列推定量^∧Ｒ_ｉを重みづけする。Ｂ_ｉ（ｔ_ｉ’）の具体例としては行列要素を主対角成分は大きく、主対角成分から離れた副対角成分であるほど小さくなるよう重み付けし、具体的には、式（２５）で表される^∧Ｒ_{ＬＯＳＳｉ}（ｔ_ｉ’）を部分帯域の相関行列推定値^∧Ｒ_ｉ（ｔ_ｉ’）に代わって用いる。以下、^∧Ｒ_{ＬＯＳＳｉ}（ｔ_ｉ’）を本発明では修正相関行列推定値と呼称する。

^∧Ｒ_{ＬＯＳＳｉ}（ｔ_ｉ’）＝　Ｂ_ｉ（ｔ_ｉ’）◎　^∧Ｒ_ｉ（ｔ_ｉ’）
　　　　＝　Ｂ_ｉ（ｔ_ｉ’）◎Ｅ｛ｘ_{Ｓｓｕｂｉ}ｘ_{Ｓｓｕｂｉ} ^Ｈ｝　・・・・（２５）
ここで記号「◎」はアダマール－シューア（Ｈａｄａｍａｒｄ－Ｓｃｈｕｒ）積（行列の同位置の要素の積で行う積）である。行列の行番号をＩ、列番号をＪで表し、^∧Ｒ_ｉの要素表記を^∧Ｒ_ｉ（ｔ_ｉ’）＝｛ｒ_ｉＩＪ（ｔ_ｉ’）｝，Ｂ_ｉの要素表記をＢ_ｉ（ｔ_ｉ’）＝｛ｂ_ＩＪ（ｔ_ｉ’）｝とする時、Ｂ_ｉ（ｔ_ｉ’）◎^∧Ｒ_ｉ（ｔ_ｉ’）＝｛ｂ_ＩＪ（ｔ_ｉ’）ｒ_ｉＩＪ（ｔ_ｉ’）｝である。

　可干渉性損失行列Ｂ_ｉの要素値の設定について図１１を用いて説明する。可干渉性損失行列Ｂ_ｉは、Ｋ_ｉ×Ｋ_ｉの実対称行列である。可干渉性損失行列Ｂ_ｉの行番号をＩ（Ｉ＝１，２，…Ｋ_ｉ）、列番号をＪ（Ｊ＝１，２，…Ｋ_ｉ）とし、各行列要素位置に式（２６）による媒介変数ｑを割り当てる。
　　　ｑ＝｜Ｒｅ｛Ｉ｝―Ｒｅ｛Ｊ｝｜＋１．０　・・・（２６）
　式（２６）において記号「｜　｜」は絶対値演算、「Ｒｅ｛　｝」は自然数を実数化する関数を表す。これにより、図１１に示すように、可干渉性損失行列Ｂ_ｉの主対角成分１１０１には、ｑ＝１．０、主対角成分の両外側の二つの第一副対角成分１１０２、１１０３にはｑ＝２．０、以下、第ｎ副対角成分にはｑ＝ｎ＋１．０が割り当てられる。以下、この変数を本実施例では対角位置変数ｑ（変域は１．０≦ｑ≦Ｒｅ｛Ｋ_ｉ｝）とする。

　次に可干渉性損失行列Ｂ_ｉの要素値を求めるために対角位置変数ｑを変数とし、変数ｑに対する定義域が１．０≦ｑ≦λ_ｉ（ｔ_i’）で、関数値が非負の値域を持つ可干渉性損失関数Ｃ_ｉ（ｔ_i’,ｑ）を考える。ここでλ_ｉ（ｔ_i’）は部分帯域ｉの再標本化時刻ｔ_i’における有効干渉長である。λ_ｉ（ｔ_i’）は近傍のチャネル間で受信信号が可干渉性を持つと考える平均的な範囲を制限するものである。対角位置変数ｑが割り当てられる可干渉性損失行列Ｂ_ｉの同じ主対角、副対角要素の値は同じＣ_ｉ（ｔ_i’,ｑ）の値で定める。例えば、λ_ｉ（ｔ_i’）の関数である実数係数ｃ_ｉ（ｔ_i’）＝ｃ_ｉ（λ_ｉ（ｔ_i’））、円周率π、余弦関数（ｃｏｓ）を用いてＣ_ｉ（ｔ_i’,ｑ）＝ｃ_ｉ（ｔ_i’）×（１．０＋ｃｏｓ（π（ｑ－１．０）／λ_ｉ（ｔ_i’）））でモデルした場合を考える。対角位置変数ｑの値が増加するほど関数値は減少し、対応する位置の可干渉性損失行列Ｂ_ｉの要素値は小さく設定される。有効干渉長λ_ｉ（ｔ_i’）が可干渉性損失行列Ｂ_ｉの次数より小さくλ_ｉ（ｔ_i’）＜Ｒｅ｛Ｋ_ｉ｝である場合、可干渉性損失関数Ｃ_ｉ（ｔ_i’,ｑ）の定義域に含まれない対角位置変数ｑに対応する可干渉性損失行列Ｂ_ｉの要素値は全てゼロとする。ｑ＜λ_ｉ（ｔ_i’）となるｑの最大の整数値をＱ_ｉとすると、Ｑ_ｉよりｑの値が大きい位置にある副対角成分の要素１１０４、１１０５の値は全てゼロとする。これにより主対角成分１１０１の行列要素の値が大きな実数で、主対角成分から離れるにつれ行列要素の値が小さくなる可干渉性損失行列Ｂ_ｉを生成することができる。

　このような対角位置変数ｑの増加に対して連続減少し、時間変化する可干渉性損失関数Ｃ_ｉ（ｔ_i’,ｑ）を元にした可干渉性損失行列Ｂ_ｉによってダイアド積ｘ_{Ｓｓｕｂｉ}ｘ_{Ｓｓｕｂｉ} ^Ｈやそれらのアンサンブル平均推定値^∧Ｒ_ｉの行列要素を重み付けする効果は二つある。

　第一は、部分アレイベクトルｘ_{Ｓｓｕｂｉ}の長さ（次数）Ｋ_ｉの増減に対して時間や空間に対して連続的な変化を実現する効果がある。受信開始からの時間の経過（超音波探触子１０からの距離が遠くなる）に従って有効相関長（ｔ_i’）が減少するように想定してモデルする場合、部分アレイベクトルの長さＫ_ｉの値を再標本化時間ｔ_i’の増加に合わせて次第に減少させる。特にＫ_ｉが比較的小さい場合、時間変化する滑らかな可干渉性損失関数Ｃ_ｉ（ｔ_i’,ｑ）を係数に持つ可干渉性損失行列Ｂ_ｉによる重みづけが無ければ、部分アレイベクトルの長さＫ_ｉを１だけ減少する場合でも、大きな確率収束の変化や部分帯域適応重みの時間不連続が引き起こされる場合がある。

　第二は、予め高い可干渉性を想定して撮像システムが用意した部分アレイベクトルのスナップショットの長さＫ_ｉの時間的な増減に対し、撮像を開始してからの有効干渉長を調整できることである。Ｋ_ｉの時間的な増減の変更は適応ビームフォーマの演算処理量の増減や分散した記憶部、たとえばメモリやレジスタの配置や積和演算の処理負荷配分の変更に直接結びつくため、それらの異なる負荷配分状態間の遷移には制約が生まれる。ある状態でのＫ_ｉを一つ減少させれば、時間あるいは空間のアンサンブル平均数を増加させた状態を実現することができるが、アンサンブル平均数の増加による処理負荷の増加と、Ｋ_ｉを減少することによる直接あるいはＱＲ分解や逆行列演算のための処理負荷の減少とは、並列信号処理プロセッサやＦＰＧＡの回路の構成や動作状態に依存して等しい処理時間の交換関係にならない。演算処理が並列信号処理プロセッサ回路で、許容された処理時間内に予定された並列演算負荷分布のスケジュールの元で実装する場合や、動作時に回路論理を再構成可能なＦｉｅｌｄ　－Ｐｒｏｇｒａｍｍａｂｌｅ　Ｇａｔｅ　Ａｒｒａｙ（ＦＰＧＡ）で、事前に動作が検証された、限られた論理の構成（コンフィグレーション）状態で実現する場合、予め設定したＫ_ｉの中で良否の判定を反映させて可干渉性損失Ｃ_ｉ（ｔ_i’，ｑ）の値を動的に変更するだけで同じ演算資源分配手順を保持することができる。

　可干渉性損失数Ｃ_ｉ（ｔ_i ’，ｑ）の値のｃ_ｉ（ｔ_i’）、λ_ｉ（ｔ_i’）を受信開始からの時間ｔ_i’で変化させる場合に、ｑ方向に広がる幅であるλ_ｉ（ｔ_i’）を１．０よりも十分細かい変化量で減少させれば、１．０単位のＫ_ｉの減少より細かい中間的な関数形状が設定できる。時間に関してｃ_ｉ（ｔ_i’）、λ_ｉ（ｔ_i’）を動的に変更して実効的にＫ_ｉを１．０より十分細かい変化量で連続減少させたと見做すことができる。

　次に可干渉性損失関数Ｃ_ｉ（ｔ_i’，ｑ）および有効干渉長λ_ｉ（ｔ_ｉ’）について、二つの部分帯域ｉ＝１、ｉ＝Ｐを想定して、同一グラフに重ねたプロットを用いて図１２で説明する。

　図１２のグラフは、横軸が、対角位置変数ｑである。また、縦軸は、可干渉性損失関数Ｃ_ｉ（ｔ_i’，ｑ）の値域であり、非負である。対角位置変数ｑは媒介変数であるので、同一のｑの値に対して複数の可干渉性損失行列Ｂ_ｉの要素値が設定される。遅延整相部２１の出力時刻系列ｔと、部分帯域ｉ＝１、Ｐの再標本化時刻ｔ_１’、ｔ_Ｐ’は異なる標本化時刻系列である。式（１７）の関係により遅延整相部２１の出力時刻番号ｔ＝τ_１に相当する時刻の値をそれぞれｔ_１’（τ_１）、ｔ_Ｐ’（τ_１）とする。ｔ_１’（τ_１）、ｔ_Ｐ’（τ_１）が解析フィルタ２１の出力する再標本化時刻番号の系列に含まれず、可干渉性損失関数Ｃ_１（ｔ_１’（τ_１）；ｑ）、Ｃ_Ｐ（ｔ_Ｐ’（τ_１）；ｑ）の直接の標本値が存在しない場合には、前後の複数時刻における再標本化時刻値から適切に補間して求める。ここで表記Ｃ_１（ｔ_１’（τ_１）；ｑ）はｔ_i’の特定値ｔ_１’（τ_１）におけるｑのみの関数であることを示す。

　図１２において、実線１２０１は、部分帯域適応処理部２２－１の解析フィルタ２１で再標本化された時刻番号ｔ_１’（τ_１）における可干渉性損失関数Ｃ_１（ｔ_１’（τ_１）；ｑ）である。小円群は可干渉性損失関数Ｃ_１（ｔ_１’（τ_１）；ｑ）をｑが自然数の位置で標本化して可干渉性損失行列Ｂ_ｉの要素値とすることを表す。破線１２０２は、部分帯域適応処理部２２－Ｐの解析フィルタ２１で再標本化された時刻番号ｔ_Ｐ’（τ_１）における可干渉性損失関数Ｃ_Ｐ（ｔ_Ｐ’（τ_１）；ｑ）である。正方形群は可干渉性損失関数Ｃ_Ｐ（ｔ_Ｐ’（τ_１）；ｑ）をｑが自然数の位置で標本化して可干渉性損失行列Ｂ_ｉの要素値とすることを表す。

　可干渉性損失関数Ｃ_１（ｔ_１’（τ_１）；ｑ）は、有効干渉長λ_１（ｔ_１’（τ_１））でゼロになり、λ_１（ｔ_１’（τ_１））は部分帯域１での部分アレイベクトルｘ_{Ｓｓｕｂ１}の長さＫ_１よりも短く設定してある。可干渉性損失行列Ｂ_１の要素は式（２６）でｑ＝Ｑ_１より大きなｑになる要素位置ではゼロである。同様に、可干渉性損失関数Ｃ_Ｐ（ｔ_Ｐ’（τ_１）；ｑ）は有効干渉長λ_Ｐ（ｔ_Ｐ’（τ_１））でゼロになり、λ_Ｐ（ｔ_Ｐ’（τ_１））は部分帯域Ｐでの部分アレイベクトルｘ_{ＳｓｕｂＰ}の長さＫ_Ｐよりも短い。可干渉性損失行列Ｂ_Ｐの要素は式（２６）のｑ＝Ｑ_Ｐより大きなｑになる要素位置ではゼロである。システム制御上は再標本化時刻番号ｔ_１’（τ_１），ｔ_Ｐ’（τ_１）については長さＫ_Ｐを想定しておきながら撮像動作時にλ_１（ｔ_１’（τ_１））、λ_Ｐ（ｔ_Ｐ’（τ_１））を適応的な調整パラメータとして動的に変更することが可能である。なお、可干渉性損失関数Ｃ_１（ｔ_１’，ｑ）、Ｃ_Ｐ（ｔ_Ｐ’，ｑ）は上述のＣ_ｉ（ｔ_i’,ｑ）＝ｃ_ｉ（ｔ_i’）×ｃｏｓ（π（ｑ－１）／λ_ｉ（ｔ_i’））以外でも任意に設定してよい。関数の変数の定義域が単一の連続区間で定義でき、ｃ_ｉ（ｔ_i’）、λ_ｉ（ｔ_i’）で例示されたようなパラメータ群で定義域の圧縮や伸長が行え、関数の定義域の境界の内部で受信時間とｑに関して連続微分可能な関数であるのが好ましい。

　次に部分帯域ｉの周波数を中心周波数などを用いて代表周波数ｆ_ｓｕｂｉで表すものとする。図１３に、部分帯域適応処理部２２－１～２２－Ｐにそれぞれ設定されている再標本化時刻番号ｔ_i ’を共通する遅延整相部２１の出力時刻系列ｔに揃える場合を考える。図１２と同じく、遅延整相部２１の出力時刻番号ｔと部分帯域番号ｉの解析フィルタ３１の出力する再標本化時刻番号ｔ_i’の関係は式（１７）で表される。ｔ＝τ_１における各部分帯域の再標本化時刻ｔ_i’（τ_１）＝（τ_１－１）／Ｔ_ｓｕｂi＋１を求め、有効干渉長λ_ｉ（ｔ_i’（τ_１））に直接の標本値が存在しない場合には、前後の複数時刻における再標本化時刻での値から適切に補間して求める。部分帯域ｉの代表周波数ｆ_ｓｕｂｉは、一定の周波数間隔で定められるとは限らないので、代表周波数ｆ_ｓｕｂｉのそれぞれを連続な変数である部分帯域周波数ｆ_ｓｕｂの離散的な標本値（塗りつぶされた小円）であると考えてｆ_ｓｕｂの関数としてλ_ｉ（ｔ_i’（τ_１））をプロットした例を示す。図１３において１３０１は有効干渉長λ_１（ｔ_１’（τ_１））をｆ_ｓｕｂ＝ｆ_ｓｕｂ１のプロット値、１３０２は有効干渉長λ_２（ｔ_２’（τ_１））をｆ_ｓｕｂ＝ｆ_ｓｕｂ２のプロット値、１３０３は有効干渉長λ_Ｐ（ｔ_Ｐ’（τ_１））をｆ_ｓｕｂ＝ｆ_ｓｕｂＰのプロット値である。これらのプロットを通過する曲線１３０４で示すようなｔ＝τ_１における関数Ｂ_λ（ｆ_ｓｕｂ）を考え、部分帯域周波数ｆ_ｓｕｂが低周波数ほど有効干渉長λ_ｉ（ｔ_ｉ’（τ_１））を長く、高周波数ほど有効干渉長λ_ｉ（ｔ_ｉ’（τ_１））を短くした例である。有効干渉長が短いほど、可干渉損失行列Ｂ_ｉは、主対角成分およびそれにごく近い副対角成分のみが非負数で、他はゼロになる。高周波帯域の受信信号の可干渉性が劣化し、有効干渉長が短くなり、雑音も増す場合に、修正相関行列推定値^∧Ｒ_{ＬＯＳＳｉ}（ｔ_ｉ’）を頑健に確率収束させることが可能になる。

　図１４に、部分帯域ｉ＝１の可干渉性損失関数Ｃ_ｉの有効干渉長λ_１（ｔ_１’）を、遅延整相部２１の出力時刻番号がｔ＝τ_１，ｔ＝τ_２，ｔ＝τ_３（τ_１＜τ_２＜τ_３）と大きくなるにつれ短縮するように定めた例を示す。それぞれの時刻での可干渉性損失関数は、実線１４０１がＣ_１（ｔ_１’（τ_１））で小円が対角位置変数ｑに対する標本値，長破線１４０２がＣ_１（ｔ_１’（τ_２））で小三角が対角位置変数ｑに対する標本値，短破線１４０３がＣ_１（ｔ_１’（τ_３））で菱形が対角位置変数ｑに対する標本値である。すなわち、時刻番号ｔ_１’が大きくなるにつれ（超音波探触子１０からの距離が深くなるにつれ）、有効干渉長λ_１（ｔ_１’）が小さくなるため、可干渉性損失行列Ｂ_１の副対角成分の値が減少するように定められている。

　部分帯域ｉのそれぞれにおいて、受信開始からの経過時間に応じて有効干渉長λ_ｉ（ｔ_ｉ’）を小さくすることにより、深い領域からの受信信号の可干渉性が劣化して有効干渉長が短くなり、雑音も増す場合に、修正相関行列推定値^∧Ｒ_{ＬＯＳＳｉ}（ｔ_ｉ’）を頑健に確率収束させることが可能になる。

　図１５は、有効干渉長λ_ｉ（ｔ_ｉ’）が部分帯域周波数ｆ_ｓｕｂと遅延整相部２１の受信開始からの経過時間ある時刻番号ｔに関してどのように変化するかを概念的に示す。理解の便をはかるため、図１４の遅延整相部２１の出力時刻番号ｔ＝τ_１，ｔ＝τ_２，ｔ＝τ_３（τ_１＜τ_２＜τ_３）をここでも用いる。遅延整相部２１の出力時刻番号ｔと解析フィルタ３１の出力する再標本化時刻番号ｔ_i’の関係は式（１７）によりｔ_i’（ｔ）＝（ｔ－１）／Ｔ_ｓｕｂi＋１である。この関係に基づき、Ｔ_ｓｕｂi倍の周期から関連付けられるτ_１、τ_２、τ_３に対応する解析フィルタ３１の出力時刻番号をｔ_i’（τ_１）、ｔ_i’（τ_２）、ｔ_i’（τ_３）とする。異なる再標本化時刻ｔ_ｉ’の上で定義される有効干渉長λ_ｉ（ｔ_ｉ’）の受信開始からの時間変化曲線を、遅延整相部２１の出力する各部分帯域に共通の時刻番号ｔに対応してどのような傾向で変化させるか例示する。図１５の１５０１は有効干渉長λ_１（ｔ_１’（ｔ））の受信開始からの時間変化曲線、１５０２は有効干渉長λ_２（ｔ_２’（ｔ））の時間変化曲線、１５０３は有効干渉長λ_３（ｔ_３’（ｔ））の時間変化曲線である。いずれも受信開始から時刻番号ｔに対応して次第に減少し、ある一定値に漸近させている。また、図１３のｔ＝τ_１における関数Ｂ_λ（ｆ_ｓｕｂ）の曲線１３０４は、図１５では曲線１３０１、１３０２、１３０３とそれぞれ１５０１、１５０２，１５０３を交点とする。曲線１３０４は部分帯域周波数ｆ_ｓｕｂが大きくなるほど減少するよう設定する。このような調整の傾向から、λ_ｉ（ｔ_ｉ’）を時間と部分帯域周波数に関連付けた２変数関数を想定して調整するのが望ましい。ｆ_ｓｕｂとｔの２変数関数を本発明では有効干渉長関数λ（ｆ_ｓｕｂ，ｔ）と呼ぶことにする。

　有効干渉長関数λ（ｆ_ｓｕｂ，ｔ）の設定方法の一例は、ｆ_ｓｕｂとｔについての区間群に区切り、バイキュービック（ｂｉ－ｃｕｂｉｃ）補間関数などで全ての代表周波数ｆ_ｓｕｂｉとｔ_ｉ’について補間計算できるようする。より具体的には、適切な間隔であらかじめ設定した座標（ｆ_ｓｕｂ、ｔ）の群について、それらの格子点での値を外部入力に従って入力部１１より定める。有効干渉長λ_ｉ（ｔ_ｉ’）の計算では、ｆ_ｓｕｂ＝ｆ_{ｓｕｂｉ，}（１７）式のｔとｔ_ｉ’時間軸方向の再標本化関係からｔ_ｉ’に対応するｔを求めてバイキュービック補間でλ_ｉ（ｔ_ｉ’）を計算することができる。既に一例として挙げた可干渉性損失関数Ｃ_ｉ（ｔ_i’,ｑ）＝ｃ_ｉ（ｔ_i’）×（１．０＋ｃｏｓ（π（ｑ－１．０）／λ_ｉ（ｔ_i’）））の場合では、ｃ_ｉ（ｔ_i’）をλ_ｉ（ｔ_i’）の関数としており、可干渉性損失行列Ｂ_ｉの要素の行番号と列番号から式（２６）により求まる対角位置変数ｑとλ_ｉ（ｔ_ｉ’）が決定されると可干渉性損失関数Ｃ_ｉ（ｔ_i’,ｑ）の値が計算でき、部分帯域番号ｉの再標本化時刻ｔ_i’における可干渉性損失行列Ｂ_ｉの要素が計算できる。

　（可干渉損失行列Ｂ_ｉによる口径制限効果）
　次に可干渉性損失行列Ｂ_ｉには、第二の受信信号群に対し、空間軸方向の口径制限効果をさらに付与することが可能である。これを図１６を用いて説明する。部分帯域ｉに対する可干渉性損失行列Ｂ_ｉは、図１１と同様に次数がＫ_ｉ×Ｋ_ｉの実対称行列である。図１６のように、行列の縦横の中心線の対１６０１、１６０２を考える。図１６では次数Ｋ_ｉは偶数を例に図示しているため、中心線１６０１の位置は、自然数で与える行列の行番号Ｉに対して中間の１／２ずれたものが基準となる。同様に、中心線１６０２の位置は、自然数で与える列番号Ｊに対して中間の１／２ずれたものが基準となる。中心線の対１６０１、１６０２の位置から行番号、列番号の両端である１またはＫ_ｉへと離れるにつれ増加する変数ｑ’を割り当てる。本実施例ではｑ’を便宜上「口径位置変数」と呼称する。図１６の例ではｑ’＝１に相当する２行（１６０３，１６０４）と２列（１６０５，１６０６）、ｑ’＝Ｑ_ｉ’＝２に相当する２行（１６０７，１６０８）と２列（１６０９，１６１０）が示してある。

　図１６のようにＫ_ｉが偶数の場合、中心線の対１６０１、１６０２の位置は、行番号Ｉ、列番号Ｊの自然数列の成す格子の中間位置にあると考え、実数の口径位置変数ｑ’はＨを行番号Ｉまたは列番号Ｊの何れかであるとして式（２７）で表される。
　ｑ’＝｜（Ｒｅ｛Ｋ_ｉ｝＋１．０）／２．０－Ｒｅ｛Ｈ｝｜＋０．５
　　　　　　　　　　　　　　　　　　　　　　　　　　・・・（２７）
Ｒｅ｛　｝は自然数を実数にする関数、｜　｜は絶対値である。この場合のｑ’の値域は１．０≦ｑ’≦Ｒｅ｛Ｋ_ｉ／２｝である。

　Ｋ_ｉが奇数の場合は、行番号Ｉ、列番号Ｊの自然数列の成す格子上の位置に、中心線対１６０１、１６０２があると考え、口径位置変数ｑ’は、式（２８）で表される。
　ｑ’＝｜（Ｒｅ｛Ｋ_ｉ｝＋１．０）／２．０－Ｒｅ｛Ｈ｝｜＋１．０
　　　　　　　　　　　　　　　　　　　　　　　　　　　・・・（２８）
この場合のｑ’の値域は１．０≦ｑ’≦（Ｒｅ｛Ｋ_ｉ｝＋１．０）／２．０である。

　可干渉性損失行列Ｂ_ｉに口径制限効果を付与する場合は、たとえば中心線の対中心線対１６０１、１６０２を基準に口径位置変数ｑ’の増加に対し単調に減少し、再標本化時刻番号ｔ_i’により変化する関数Ｄ_i（ｔ_i’，ｑ’）を定め、Ｉ行Ｊ列の要素値に、Ｄ_i（ｔ_i’，ｑ’（Ｉ））×Ｄ_i（ｔ_i’，ｑ’（Ｊ））の重みを与える。関数Ｄ_i（ｔ_i’，ｑ’）を本発明では口径制限関数と呼称する。ｑ’がＱ_ｉ’より大きくなる行列要素値は、全てゼロとする。口径制限関数の時間変化を制御するため、新たに有効口径長α_ｉ（ｔ_ｉ’）を考える。α_ｉは受信開始からの再標本開時刻ｔ_ｉ’で変化する。例えば、口径制限関数Ｄ_i（ｔ_i’，ｑ’）は変数ｑ’が区間１．０≦ｑ’＜α_ｉ（ｔ_ｉ’）／２＋１．０ではＤ_i（ｔ_i’，ｑ’）＝１．０、区間α_ｉ（ｔ_ｉ’）／２＋１．０≦ｑ’≦α_ｉ（ｔ_ｉ’）＋１．０ではＤ_i（ｔ_i’，ｑ’）＝０．５ｃｏｓ（２π（ｑ’－１－α_ｉ（ｔ_ｉ’）／２）／α_ｉ（ｔ_i’））＋０．５と定義出来る。このＤ_i（ｔ_i’，ｑ’）は有効口径長α_ｉ（ｔ_ｉ’）の半分が１．０、残りが滑らかなテーパーを持つようにしている。

　図１１から図１５で説明した可干渉性損失行列Ｂ_ｉに口径制限効果を与えるため、再標本化時刻ｔ_i’におけるＢ_ｉのＩ行Ｊ列要素ｂ_iＩＪ（ｔ_i’）を式（２９）で表わす。
ｂ_iＩＪ（ｔ_i’）＝Ｄ_i（ｔ_i’，ｑ’（Ｉ））×Ｄ_i（ｔ_i’，ｑ’（Ｊ））
×Ｃ_ｉ（ｔ_i’，ｑ）　　　　　　　　・・・・（２９）
　ｑは式（２６）、ｑ’は式（２７）または式（２８）により計算できる。

　次に、口径制限関数Ｄ_i（ｔ_i’，ｑ’（Ｈ））と、有効口径長α_ｉ（ｔ_ｉ’）について、代表的な二つの部分帯域ｉ＝１、ｉ＝Ｐを想定して図１７を用いて説明する。

　図１７のグラフは、横軸が口径位置変数ｑ’である。また、縦軸は口径制限関数Ｄ_ｉ（ｔ_i’，ｑ’（Ｈ））の値域であり、０から１．０の間の値をとる。図１７では、遅延整相部２１の出力する時刻番号ｔ＝τ_１に対応する、式（１７）で計算される再標本化時刻番号ｔ_１’（τ_１）における部分帯域１の口径制限関数Ｄ_１（ｔ_１’（τ_１），ｑ’）を実線１７０１で表す。二重丸はｑ’の各値における標本値である。再標本化時刻番号ｔ_Ｐ’（τ_１）における部分帯域Ｐの口径制限関数Ｄ_Ｐ（ｔ_Ｐ’（τ_１），ｑ’）を破線１７０２で表した。二重の正方形はｑ’の各値における標本値である。可干渉性損失関数Ｃ_ｉ（ｔ_i’,ｑ’）の場合と同じように、ｔ_１’（τ_１）やｔ_Ｐ’（τ_１）が再標本化時刻の系列で標本化されていなければ、Ｄ_１（ｔ_１’（τ_１），ｑ’）やＤ_Ｐ（ｔ_Ｐ’（τ_１），ｑ’）は時系列上の前後の複数の標本値から補間して求められる。実線１７０１で示すＤ_１（ｔ_１’（τ_１），ｑ’）は、有効口径長α_１（ｔ_１’（τ_１））でゼロになり、α_１（ｔ_１’（τ_１））は、式（２７）のｑ’の値域の上限Ｋ_１／２よりも大きい。同様に、Ｄ_Ｐ（ｔ_Ｐ’（τ_１），ｑ’）は、有効口径長α_Ｐ（ｔ_Ｐ’（τ_１））でゼロになり、α_Ｐ（ｔ_１’（τ_１））も、式（２７）のｑ’の値域の上限Ｋ_Ｐ／２よりも大きい。

　このように、可干渉性損失行列Ｂ_ｉに口径制限効果を与えることにより、システム制御上は時刻番号ｔ_ｉ’については長さＫ_ｉを想定しておきながら、動作時に遅延整相部２１の出力する時刻番号ｔについてα_ｉ（ｔ_ｉ’（ｔ））を調整パラメータとして連続的に有効部分を変化させることができる。

　図１８に遅延整相部２１の出力する時刻番号ｔ＝τ_１における部分帯域周波数ｆ_ｓｕｂに対する有効口径長α_ｉ（ｔ_ｉ’（τ_１））の設定例を示す。各部分帯域の再標本化時刻ｔ_i’（τ_１）＝（τ_１－１）／Ｔ_ｓｕｂi＋１を求め、有効口径長α_ｉ（ｔ_i’（τ_１））に直接の標本値が存在しない場合には、前後の複数時刻における再標本化時刻での値から適切に補間して求める。部分帯域ｉの代表周波数ｆ_ｓｕｂｉは、一定の周波数間隔で定められるとは限らないので、代表周波数ｆ_ｓｕｂｉのそれぞれを連続な変数である部分帯域周波数ｆ_ｓｕｂの離散的な標本値（塗りつぶされた小円）であると考えてｆ_ｓｕｂの関数としてα_ｉ（ｔ_i’（τ_１））をプロットした例を示す。図１８において１８０１は有効口径長α_１（ｔ_１’（τ_１））をｆ_ｓｕｂ＝ｆ_ｓｕｂ１のプロット値、１８０２は有効口径長α_２（ｔ_２’（τ_１））をｆ_ｓｕｂ＝ｆ_ｓｕｂ２のプロット値、１８０３は有効口径長α_Ｐ（ｔ_Ｐ’（τ_１））をｆ_ｓｕｂ＝ｆ_ｓｕｂＰのプロット値である。これらのプロットを通過する曲線１８０４で示すようなｔ＝τ_１における関数Ｂ_α（ｆ_ｓｕｂ）を考え、部分帯域周波数ｆ_ｓｕｂが低周波数ほど有効口径長α_ｉ（ｔ_ｉ’（τ_１））を長く、高周波数ほど有効干渉長α_ｉ（ｔ_ｉ’（τ_１））を短くした例である。有効口径長α_ｉは受信素子の周波数に依存した指向性に関連して定められる。電気音響変換器アレイ１０ａの受信素子（トランスデューサ）は全ての部分帯域受信信号に対し同じ物理的寸法を持つため、相対的に低周波では広い受信指向角度、高周波には狭い受信指向角度を持つ。受信焦点位置を受信指向角範囲より外側に捉える素子を受信信号群に含めると、無関係の到来方向からの信号を主に含むようになり望ましくない。解析フィルタ３１の空間軸方向の再標本化周波数を決めるＳ_ｓｕｂと受信素子の指向性に依存して有効干渉長α_ｉが決定され、不要な受信応答を部分帯域周波数に応じて制限することができる。

　図１９に、部分帯域ｉ＝１の口径制限関数Ｄ_ｉの有効口径長α_１（ｔ_１’）を遅延整相部２１の出力時刻番号がｔ＝τ_１，ｔ＝τ_２，ｔ＝τ_３（τ_１＜τ_２＜τ_３）と大きくなるにつれ増加するように定めた例を示す。実線１９０１は口径制限関数Ｄ_１（ｔ_１’（τ_１））で二重丸が口径位置変数ｑ’に対する標本値である。有効口径長α₁（ｔ₁’（τ_１））はＫ_１／２よりもわずかに大きい。破線１９０２はＤ_１（ｔ_１’（τ_２））で二重三角が口径位置変数ｑ’に対する標本値である。有効口径長α₁（ｔ₁’（τ_２））はＫ_１／２よりも１，２程度大きい。一点鎖線１９０３がＤ_１（ｔ_１’（τ_３））で二重菱型が口径位置変数ｑ’に対する標本値である。有効口径長α₁（ｔ₁’（τ_３））はＫ_１／２よりも非常に大きく、二重菱型の標本値は全て１．０になっている。時刻番号ｔ_１’が大きくなるにつれ（超音波探触子１０からの距離が大きくなるにつれ）、受信焦点を指向角度に収めることができる素子数が増えるため、有効口径長α₁（ｔ₁’）は時間と共に増加する（１９０１から１９０２への変化）。α₁（ｔ₁’）／２をＫ₁／２－１が超過すると口径制限関数はＤ_１（ｔ_１’（τ_３））のようにｑ’に依存せず一定値１．０になる（１９０３）。このように、経過時間に応じて有効口径長を周波数に応じて適切に調整することにより、浅い部分の受信信号群では高周波信号の不要応答を削減して取得することができる。

　図２０は、有効口径長α_ｉ（ｔ_ｉ’）が部分帯域周波数ｆ_ｓｕｂと遅延整相部２１の受信開始からの経過時間ある時刻番号ｔに関してどのように変化するかを概念的に示す。理解の便をはかるため、図１４、図１９の遅延整相部２１の出力時刻番号ｔ＝τ_１，ｔ＝τ_２，ｔ＝τ_３（τ_１＜τ_２＜τ_３）をここでも用いる。遅延整相部２１の出力時刻番号ｔと解析フィルタ３１の出力する再標本化時刻番号ｔ_i’の関係は式（１７）により関連付けられる。τ_１、τ_２、τ_３に対応する解析フィルタ３１の出力時刻番号をｔ_i’（τ_１）、ｔ_i’（τ_２）、ｔ_i’（τ_３）とする。異なる再標本化時刻ｔ_ｉ’の上で定義される有効口径長α_ｉ（ｔ_ｉ’）の受信開始からの時間変化曲線を、遅延整相部２１の出力する各部分帯域に共通の時刻番号ｔに対応してどのような傾向で変化させるか例示する。図２０の２００１は有効口径長α_１（ｔ_１’（ｔ））の受信開始からの時間変化曲線、２００２は有効口径長α_２（ｔ_２’（ｔ））の時間変化曲線、２００３は有効口径長α_３（ｔ_３’（ｔ））の時間変化曲線である。いずれも受信開始の所定値から時刻番号ｔに対応して比例的に増加させている。また、図１８のｔ＝τ_１における関数Ｂ_α（ｆ_ｓｕｂ）の曲線１８０４は、図２０では曲線２００１、２００２、２００３とそれぞれ点１８０１、１８０２，１８０３を交点として交差する。曲線１８０４は部分帯域周波数ｆ_ｓｕｂが大きくなるほど減少するよう設定する。このような調整の傾向から、α_ｉ（ｔ_ｉ’）を時間と部分帯域周波数に関連付けた２変数関数を想定して調整するのが望ましい。ｆ_ｓｕｂとｔの２変数関数を本発明では有効口径長関数α（ｆ_ｓｕｂ，ｔ）と呼ぶことにする。

　有効口径長関数α（ｆ_ｓｕｂ，ｔ）の設定方法の一例は、ｆ_ｓｕｂとｔについての区間群に区切り、バイキュービック（ｂｉ－ｃｕｂｉｃ）補間関数などで全ての代表周波数ｆ_ｓｕｂｉとｔ_ｉ’について補間計算できるようする。より具体的には、適切な間隔であらかじめ設定した座標（ｆ_ｓｕｂ、ｔ）の群について、それらの格子点での値を外部入力に従って入力部１１より定める。有効口径長α_ｉ（ｔ_ｉ’）の計算では、ｆ_ｓｕｂ＝ｆ_{ｓｕｂｉ，}（１７）式のｔとｔ_ｉ’時間軸方向の再標本化関係からｔ_ｉ’に対応するｔを求めてバイキュービック補間でα_ｉ（ｔ_ｉ’）を計算することができる。既に一例として挙げた口径制限関数Ｄ_i（ｔ_i’，ｑ’）では変数ｑ’が区間１．０≦ｑ’＜α_ｉ（ｔ_ｉ’）／２＋１．０ではＤ_i（ｔ_i’，ｑ’）＝１．０、区間α_ｉ（ｔ_ｉ’）／２＋１．０≦ｑ’≦α_ｉ（ｔ_ｉ’）＋１．０でＤ_i（ｔ_i’，ｑ’）＝０．５ｃｏｓ（２π（ｑ’－１－α_ｉ（ｔ_ｉ’）／２）／α_ｉ（ｔ_i’））＋０．５と定義したが、変数の定義域をα_ｉ（ｔ_i’）の関数としており、可干渉性損失行列Ｂ_ｉの要素の行番号と列番号から式（２７）または式（２８）により求まる口径位置変数ｑ’とα_ｉ（ｔ_ｉ’）が決定されると口径制限関数Ｄ_ｉ（ｔ_i’,ｑ’）の値が計算でき、部分帯域番号ｉの再標本化時刻ｔ_i’における可干渉性損失行列Ｂ_ｉの要素が式（２９）で計算できる。この可干渉性損失行列Ｂ_ｉにより、推定値の相関行列（あるいは共分散行列）^∧Ｒ_{ＬＯＳＳｉ}の有効次元を、時刻番号ｔの変化に応じて連続的に変化させることができる。

　受信開始からの時刻番号ｔの進行に合わせ、生体内の不均一な伝搬で超音波エコーのコヒーレンスが次第に失われる場合がある。空間軸方向の部分アレイベクトルの次数Ｋ_ｉを低減すると限られたアレイ数あるいは空間再サンプル数の範囲でアンサンブル平均数を空間軸方向に増加させることができる。さらに時間軸方向でもアンサンブル平均を行い、確率収束を高めて頑健な推定値^∧Ｒ_{ＬＯＳＳｉ}を得ることができる。これらの条件の時間的変化の状況は、一般に周波数ごとに異なると考えられるので、可能な範囲内でアンサンブル平均数を増減することもできる。図１６では、部分アレイの長さＫ_ｉの中から中央部分だけの範囲の選択を行う例を示したが、中心線の対１６０１、１６０２を移動することで、部分アレイの長さＫ_ｉの範囲内でＫ_ｉより小さな部分アレイ長さを選択することも可能である。元の部分アレイ長さＫ_ｉの選択を図９、図１０のように固定しておきながらも、図１１の可干渉性損失関数Ｃ_i（ｔ_i’，ｑ）と図１６のような口径制限関数Ｄ_i（ｔ_i’，ｑ’）での口径制限とを階層的に組み合わせることで、滑らかなアンサンブル平均を自由度高く実現できる。これらは、周波数帯域ごと、時間ごとに変更できる。

　ソフトウェアによるシストリックアレイ３５の実装や、動的なコンフィギュレーションが可能なＦＰＧＡ回路などでの実現は、部分アレイ次数の低減により演算能力や演算時間を時間軸、空間軸方向の標本数の増加などに再配分を行うことができると考えられる。しかし、一般に演算処理が並列演算負荷を分布させた有限のスケジュール状態の元で実装されるような場合や、ＦＰＧＡで事前に有限の限られたコンフィグレーション状態の回路での実現では、部分アレイ次数と時空間平均化数の組み合わせの全てについて並列演算負荷の分布やコンフィグレーションを用意するのは実際的でない。

　並列演算負荷の分布や回路コンフィグレーション状態は、部分アレイ次数と時間および／または空間平均化数の増減に合わせて連続的に処理中に遷移するには不連続が伴う。同じ並列演算負荷分布やコンフィグレーション状態の中で滑らかにそれらの実効演算が変化することが望ましい。連続微分可能で滑らかな有効干渉長関数λ（ｆ_ｓｕｂ，ｔ）や有効口径長関数α（ｆ_ｓｕｂ，ｔ）のモデルに基づき可干渉性損失行列Ｂ_ｉの係数選択と階層的に組み合わせることで連続的な遷移を実現できる。

　時刻番号ｔの進行により、超音波エコーのコヒーレンスが失われる場合には、相関行列あるいは共分散行列に対角負荷を付与することにより、確率的に生起し得る逆行列計算あるいはそれに準ずる処理の計算不能状態を回避することができる。相関行列推定量^∧Ｒ_ｉ（ｔ_ｉ’）に式（３０）により、部分帯域ごとの仮想的な空間白色雑音電力σ_i ^２（ｔ_ｉ’）を与える。
　^∧Ｒ_ｉＥ（ｔ_ｉ’）＝^∧Ｒ_ｉ（ｔ_ｉ’）＋σ_i ^２（ｔ_ｉ’）Ｉ’・・・・（３０）
ここでＩ’は次数がＲ_ｉ（ｔ_ｉ’）と同じで対角成分の要素が全て１でそれ以外は０の単位行列とする。σ_i（ｔ_ｉ’）は時間ｔ_ｉ’に応じて変化するスカラ量である。

　これにより、送信直後の受信可変口径処理で使用する受信チャンネル数がごく少数に限られ、上記空間および／または時間のアンサンブル平均点数が十分な確率収束に至るよう確保することが不能な状況においても、直接逆行列計算を行う適応処理を破綻させることなく実現することができる。また、再帰型適応アルゴリズムにおいても、初期開始状態を提供することができる。これらは、周波数に依存するので周波数ごとのσ_i（ｔ’_ｉ）の設定により、異なる要請を充足することができる。

　以上の例では、受信開始からの経過時間に応じて、部分帯域ごとの相関行列あるいは共分散行列の時間または空間の更新間隔を増加させることで、逆行列計算あるいはそれに準ずる処理の時間当たり更新頻度を低下させ、単位時間当たり演算負荷を低減することができる。特に、逆行列計算処理あるいはそれに準ずる処理は、アルゴリズムによっては部分アレイ長さに対して３乗の規模で演算量が増加するため、高い演算負荷を必要とする。相関行列あるいは共分散行列を、より大域的深さ範囲で緩慢な変化をするものと仮定する場合に、十分な確率収束を得た頑健な適応処理を、部分周波数帯域ごとに独立して行うことができる。

　図１の入力部１１にこれらの有効干渉長関数λ（ｆ_ｓｕｂ，ｔ）、有効口径長関数α（ｆ_ｓｕｂ，ｔ）、可干渉性損失関数Ｃ_ｉ（ｔ_i’,ｑ）、口径制限関数Ｄ_i（ｔ_i’，ｑ’）の関数式やそれらの計算式を上記格子点での入力値から補間計算できるよう装置を構成しても良い。入力値や補間関数を入力部１１より入力し、それらの妥当性を確認するための図１２、図１３、図１４、図１５のようなグラフや関数を制御部１２やバックエンド処理部１８で計算して表示部１９で表示することもできる。また、それらにより可干渉性損失行列Ｂ_ｉの要素や可干渉性損失行列Ｂ_ｉの要素から計算した情報を必要に応じて撮像断面などと同時表示され、周波数や受信開始からの時間に対していかなる条件の元で適応ビームフォーミング処理を用いて撮像されているかが装置使用者に明白になるように装置を構成しても良い。これらの条件の表示により、適応処理を受けた画像が生成される際の頑健性（Ｒｏｂｕｓｔｎｅｓｓ）条件などを明確に提示することができる。

　＜シミュレート実験による傍証＞
　本実施形態の部分帯域ごとに部分帯域適応処理部２２－１～２２－Ｐで適応処理する構成の有効性を、時間領域差分法、有限差分時間領域法（ＦＤＴＤ法）のモデル計算で得た音波伝搬シミュレートにより確認した。その結果を示す。適応フィルタとしては、相関行列を用いるＭＶＤＲ適応ビームフォーマの変形であるＡＰＥＳビームフォーマを用いた。

　図２１は、シミュレート計算に用いた生体モデル２１０１である。生体モデル２１０１は、２次元であり、縦４５ｍｍ、横３５ｍｍである。生体モデル２１０１の第一媒質領域２１０３に、複数の第二媒質領域２１０２を互いの重複を許して分布させた。第二媒質領域２１０２は、０ｍｍ～１ｍｍの間の乱数発生により定めた半径を持ち、第一媒質領域２１０３内の、乱数発生により決定した中心座標位置を持つ円形領域である。第一媒質、第二媒質は、ともに密度１ｇ／ｃｍ^３、減衰率は０．６ｄＢ／ＭＨｚ／ｃｍとした。第一媒質の音速は１５４０ｍ／ｓで第二媒質の音速は１４５０ｍ／ｓで減衰率は０．６ｄＢ／ＭＨｚ／ｃｍとした。二つの媒質の違いは音速だけである。また、必要に応じて第二媒質を配置しない比較例の均質な生体モデルも用いた。

　また、電気音響変換器アレイ１０ａの模擬として受信面配列モデル２１０５を置いた。受信面配列モデル２１０５は、直線アレイとし、横幅が０．３０８ｍｍの受音面を持つ９６個の受信素子の受信面が配置周期０．３２２ｍｍで並んでいるとした。音響伝搬計算には、受信素子の材質は計算要素とせず、第一媒質領域２１０３の計算要素のままであり、個々の受信素子表面にある計算要素格子に到達した時間波形応答を総和して、受信素子の個々の受信信号時間応答を模擬することにより、電気音響変換器アレイ１０ａの特性の影響が反映しないよう理想化した。

　実際の超音波診断装置のパルスエコー法では、超音波の送信及びエコー受信過程を経るが、本発明の効果を的確に示すために、ここではエコー受信過程だけを模擬した。エコーの発生を模擬するため、モデル２１０１中にはそれぞれ同時に円筒波パルスを発生する音源２１０４を配置した。受信ビームフォーマ１７の動作は計算により求めた。音源２１０４は、電気音響変換器アレイ１０ａの受信面配列モデル２１０５の表面を基準に、深さ（距離）２０ｍｍ、３０ｍｍ、４０ｍｍ、領域の横幅中央を基準に１０．３ｍｍ離した２×３＝６点の格子位置に設定されている。仮想的な９６個の受信素子で得られる受信信号群の中から、Ｎ＝３２チャンネルを選び出して適応ビームフォーミングを模擬した。９６個の受信素子の中から３２個を選ぶ位置を順次横に移動させながら像再生を行う。

　６個の音源２１０４は、図２２に示す時間波形２２０１および図２３に示す時間フーリエスペクトル（周波数スペクトル）２３０１を示す中心周波数５ＭＨｚの音圧パルスを発生するとして計算した。図２１の生体モデル２１０１は、第二媒質領域２１０２の最大直径の値域（凡そ０ｍｍ～１ｍｍ）に図２２の音圧パルス２２０１の持つ第一媒質中での波長０．３ｍｍが含まれ、図２３の周波数スペクトル２３０１の周波数の各部分帯域、および、音源２１０４と受信面配列モデル２１０５の表面との間の波の伝搬飛行経路の違い、に応じて不均一散乱の影響が、周波数および深さで蓄積変化するよう期してモデルを構成した。

　音場伝搬シミュレーションで得られた受信信号群を、Ｈａｎｎｉｎｇ包絡線の中心周波数が部分帯域の中心周波数がｆ_ｓｕｂ１＝３．０ＭＨｚ、ｆ_ｓｕｂ２＝３．８ＭＨｚ、ｆ_ｓｕｂ３＝４．６ＭＨｚ、ｆ_ｓｕｂ４＝５．４ＭＨｚ、ｆ_ｓｕｂ５＝６．２ＭＨｚの５つの受信信号群にバンドパスフィルタで分けた。部分帯域は互いに部分的に重畳している。５つの部分帯域ごとの、空間平均化のための部分空間ベクトルの長さＫ_ｉの設定値をｋｓｕｂとし、ｋｓｕｂ＝４，８，１６の条件で前方後方空間平均化と時間平均化を行うＡＰＥＳ適応ビームフォーマ処理を行った。前方後方空間平均化は、部分アレイによる空間軸方向のアンサンブル平均を部分アレイを反転させて逆向きにも行うものである。また、画像値の統計比較を行うため、全ての条件でビームフォーマ処理は同じ画素位置および縦横の画素数で計算した。

　図２４にｆ_ｓｕｂ１＝３．０ＭＨｚ、ｋｓｕｂ＝４の場合の像を等高線図で示す。第二媒質による不均一散乱の影響が無い場合で、適応処理を行わず、口径重み付与も行わない遅延加算法により得られた画像内の最大値を用いて、図２４の等高線の値は規格化されている。図２４の等高線は、この遅延加算だけの最大値に対して－２０ｄＢの水準を太線、－４０ｄＢの水準を細線で表示している。また、図２４は、ｄＢ値の文字表示の判読の便を図るため、横方向に引き伸ばしてある。図２１の点音源２１０４を含め、６つの点音源の近傍位置に－２０ｄＢの等高線領域を確認できるが、第二媒質による不均一散乱の影響により形状が崩れた点像となっている。深部の音源像になるに従い、－４０ｄＢの点像が、深さ方向（Ｒａｎｇｅ、図２１の深さに一致）およびアレイ方向（Ａｚｉｍｕｔｈ、図２１のアレイ位置）に領域が拡大していることを確認できる。また、図２１の６つの点音源設定位置と比べ、図２４の－２０ｄＢの等高線領域が、深い側（Ｒａｎｇｅの大きい側）へ移動している。これは、第一媒質の音速を想定した遅延整相処理の設定音速１５４０ｍ／ｓに比べ、音速の遅い第二媒質が空間を大域的に見て低い音速媒質のように変化させたためと理解できる。

　次に、空間軸方向の再サンプル（再標本化）の効果を確認数する為、同じｆ_ｓｕｂ１＝３．０ＭＨｚ、ｋｓｕｂ＝４の場合について、空間軸方向のみに再サンプルを行い、９６チャンネルをその４分の３に低減を行なってから、前方後方空間平均化と時間平均化を行うＡＰＥＳ適応ビームフォーマ処理を行った。その結果を図２５に示す。空間軸方向の再標本化には、カイザー窓を用いた。図２５に示される結果は、前方後方空間平均化数は減少したが、再サンプルを行わない場合の図２４の結果と比べて、ほぼ同等の結果を得られることが確認できた。空間軸方向の異なる再サンプル状態で得られる適応処理後の結果の妥当性については、適応処理後の目的に合わせて適切に選択される。

　次に、部分帯域の周波数による比較を行った結果を示す。部分帯域の中心周波数ｆ_ｓｕｂ１、ｆ_ｓｕｂ３、ｆ_ｓｕｂ５の３帯域を選び出し、部分アレイベクトルの長さｋｓｕｂは全ての部分帯域で１６点とした。空間の再サンプルは行っていない。それぞれの部分帯域の画像の画素値を最大０ｄＢから０．５ｄＢごとの標本区間群に帰属させて計数し、－６０ｄＢ～０ｄＢまでのヒストグラム分布を求めた。３つの部分帯域のヒストグラム曲線を重ねて描いたものを図２６に示す。縦軸は頻度(Incidence)、横軸は画素値すなわち対数規格化した強度(Normalized　Intensity)である。部分帯域の中心周波数がｆ_ｓｕｂ１の画像のヒストグラム曲線は一点鎖線、ｆ_ｓｕｂ３は鎖線、ｆ_ｓｕｂ５は実線で描いてある。

　図２６により、部分帯域の中心周波数が増加するに従って、ヒストグラム曲線が０ｄＢ側に寄って分布していることがわかる。このことから周波数が増すにつれ、画素全体で高輝度側に寄ったことが確認できる。図２６では判別が困難になるため重ねて図示しないが、ｆ_ｓｕｂ２、ｆ_ｓｕｂ４を重ねた場合でも周波数に対する傾向は同じであった。ヒストグラム最頻値の画素値と最大値の間が大きく離れることを画像の諧調ダイナミックレンジの改善指標とするならば、部分帯域の周波数が増すにつれ分布が大きな画素値に偏っており、分布の山の部分がほぼ画像の背景水準を形成するとすれば、点音源像の識別を劣化させていることになる。

　この結果から、部分帯域周波数について分割したＡＰＥＳ適応ビームフォーマ処理結果を加算器２３でコヒーレント合成（そのまま加算）あるいはインコヒーレント合成（絶対値をとって加算）する場合、ヒストグラム最頻値の部分の大半の画素が周波数に関して無相関な情報で形成されると考え、画素値の分布の周波数で異なる特徴量を考慮してから再合成するのが有利であると考えられる。図２６の場合、部分帯域ごとのＡＰＥＳ適応ビームフォーマ処理出力に対し、振幅の対数圧縮前に高い部分帯域周波数であるほど小さな重み係数を乗じてコヒーレントあるいはインコヒーレント加算して最頻値付近の分布を周波数間で揃えて加算合成したり、対数圧縮した後のものに重み係数を乗じて周波数ごとに対数軸を相対的に伸縮してからインコヒーレント合成することで画像の強調を行うことができる。

　次に、部分帯域の中心周波数ｆ_ｓｕｂ１、ｆ_ｓｕｂ３、ｆ_ｓｕｂ５の３帯域のそれぞれについて、不均一散乱の影響が異なる３つの点音源の深さ位置の近傍それぞれについて画像を深さの区間に分け、周波数ごと、それぞれの深さでの最大値で規格化したヒストグラム曲線を図２７～図２９に示す。３つの点音源の近傍の深さの区間とは、図２５における深さ方向（Ｒａｎｇｅ）の２～３ｃｍ、３～４ｃｍ、４～４．５ｃｍである。これらの区間の開始深さで代表して図２７～図２９の凡例では、順にＤｅｐｔｈ＝２［ｃｍ］、Ｄｅｐｔｈ＝３［ｃｍ］、Ｄｅｐｔｈ＝４［ｃｍ］と記している。

　図２７は部分帯域の中心周波数ｆ_ｓｕｂ１、図２８はｆ_{ｓｕｂ３、}図２９はｆ_ｓｕｂ５の結果を示す。図２７～図２９において、一点鎖線は最も浅い深さ２ｃｍの音源、鎖線は中間の深さ３ｃｍの音源、実線は最も深い４ｃｍの音源の区間の結果である。図２７～図２９の何れも、ヒストグラムは音源が深くなるほど分布が高輝度の０ｄＢ側に寄る傾向が得られた。さらに、同じ深さごとに比較すると部分帯域の中心周波数が高いと分布が高輝度の０ｄＢ側に寄る傾向が見いだせる。たとえば、最も浅い２ｃｍの音源の区間（Ｄｅｐｔｈ＝２［ｃｍ］）のヒストグラムの最頻値は、図２７の中心周波数ｆ_ｓｕｂ１で尖頭位置２７０１からおおよそ－４２ｄＢ近傍、図２９の中心周波数ｆ_ｓｕｂ５で尖頭位置２９０１からおおよそ－３８ｄＢ近と読み取ることができる。最も深い４ｃｍの音源の区間（Ｄｅｐｔｈ＝４［ｃｍ］）のヒストグラムの最頻値は図２７の中心周波数ｆ_ｓｕｂ１で尖頭位置２７０２からおおよそ－２８ｄＢ近傍、図２９の中心周波数ｆ_ｓｕｂ５で尖頭位置２９０２からおおよそ－２２ｄＢ近傍と読み取ることができる。こうした例から、ヒストグラム最頻値あるいは他の指標で画素値の分布範囲を考慮してから再合成することで、より望ましい合成を実現できると考えられる。

　図２６からは周波数ごとに分割されたＡＰＥＳ適応ビームフォーマ処理結果の合成には、周波数ごとの重み付けが有効であろうこと、図２７～図２９の結果からは、それらの周波数ごとの重み付けの配分は、深さごとに変化させるべきであることが示された。周波数ごとに重み付けを行わずに合成すれば最も最頻値付近の分布が高輝度側の周波数成分が画像の諧調を支配する可能性がある。また、深さごとに変化するよう対数圧縮した後の出力に重み係数を乗じて合成すれば、深部での諧調ダイナミックレンジの分布を低輝度側へ広げ、目視による判定や読影を良好にすることが期待できる。

　次に空間平均化のための部分空間ベクトルの長さＫ_ｉの設定値であるｋｓｕｂに対する依存について、図３０に比較した例を示す。部分帯域の中心周波数ｆ_ｓｕｂ１、最も深い４ｃｍの近傍で規格化したヒストグラムの比較を示す。ｋｓｕｂは、一点鎖線が４、破線が８、実線が１６である。ヒストグラムの区間ごとの観察では部分的な順序逆転はあるが、全体としてｋｓｕｂが小さいほど分布が高輝度の０ｄＢ側に寄る傾向が得られた。ｋｓｕｂが増すと適応ビームフォーマの計算負荷が大きく増す場合には、ｋｓｕｂの設定値を周波数ごと、深さ（すなわち受信時刻ごと）に効率よく設定することが重要である。ある限られた演算回路や演算時間を仮定すると、これらの周波数ごと、時刻ごとの応答の違いを考慮してたとえば高い周波数の部分帯域にｋｓｕｂを増すなどの配分を採ることができる。あるいは、生体非線形の音響伝搬効果や非線形の調波成分の発生を積極的に利用する造影剤の撮像を考慮して特定の帯域にｋｓｕｂの増減を与える、すなわち空間平均化のための部分空間ベクトルの長さＫ_ｉの調整状態を作り出すことができる。

　図２４、図２５では、部分帯域の中心周波数をｆ_ｓｕｂ１、ｋｓｕｂを４点にした結果について4分の３（０．７５）倍の空間再サンプルを行う前後の画像の等高図で示したが、これらの２画像の深さ全体について、ヒストグラムを重ねて描いた結果を図３１に示す。点線が図２４の空間再サンプルを行わない（No　resample)場合、実線が図２５の空間再サンプルを行なった(resample　0.75)場合である。空間再サンプルを行わない場合の最頻値３１０１は、－４１ｄＢ程度、行った場合の最頻値３１０２は、－４３ｄＢ程度と読み取ることができ、空間再サンプルを行なうことで分布が－２ｄＢ移動していることが示された。これらの結果から部分帯域ごとに空間再サンプルを行なうことで良好な合成結果が得られることが期待される。

　以上のシミュレーション結果から、生体を模擬した不均一伝搬の周波数依存性を示す例を模擬し、本発明の効果と有効性について傍証を示した。

１０　超音波探触子
１０ａ　電気音響変換器アレイ
１０ｂ　送信パルス
１０ｃ　被検体内の反射体
１０ｄ　パルスエコー
１１　入力部
１２　制御部
１２－１　データバス
１３　送信ビームフォーマ
１４　送信回路
１５　送受分離回路
１６　受信回路
１７　受信ビームフォーマ
１８　バックエンド処理部
１９　表示部
２０　受信波形記憶部
２１　遅延整相手段
２２　適応フィルタ部
２２－１～２２－Ｐ　部分帯域適応処理部
２３　加算器
３１　解析フィルタ
３２　部分帯域適応処理部
３３　合成フィルタ
３４　部分帯域信号（第二の時間信号群）
３５　部分帯域適応信号
４１　適応制御部
４１－１　部分帯域波形読み出しアドレス信号
４１－２　多重化アドレス信号
４１－３　係数多重化アドレス信号
４１－４　適応出力整序アドレス信号
４１－５　シーケンステーブル
４１－６　動作モード指令（Ｍ_ＳＹＳ　）
４２　部分帯域波形記憶部
４３　シストリックアレイ多重化器
４４　係数テーブル記憶部
４５　シストリックアレイ
４６　部分帯域適応出力記憶部
５１　境界セル
５２　遅延セル
５３　内部セル
５４　最終セル
５５　入力データシーケンス
５６　アレイ入力側セル群
ｘ_ｉｎ　セルへの要素入力
ｘ_ｏｕｔ　セルからの要素回転出力
ｃ　Ｇｉｖｅｎｓ回転余弦
ｓ　Ｇｉｖｅｎｓ回転正弦
ｒ　上三角行列要素
γ_ｉｎ　変換係数入力
γ_ｏｕｔ　変換係数出力
β　忘却係数（指数重み係数）
μ　拘束条件値
ｗ_ｉｎｂ，ｗ_ｉｎｒ，重みづけベクトル係数
７０　演算部
８１　解析フィルタ入力点
８２　解析フィルタ出力点
８３　解析フィルタ３１の出力点（ｘ_{ｔ’,ｊ’}）
８４　入力の格子点基準位置
８５　解析フィルタ入力の標本値群
９１　部分アレイベクトルの第一要素位置
９２　部分帯域１のアンサンブル平均範囲
９３　部分帯域１の部分アレイベクトル（ｘ_{Ｓｓｕｂ１}）
Ｋ_１、Ｋ_ｉ、Ｋ_Ｐ　部分アレイベクトルのスナップショットの長さ
Ｌ_１、Ｌ_Ｐ　時間アンサンブル平均数（標本数）
１０１　部分アレイベクトルの第一要素位置
１０２　部分帯域１のアンサンブル平均範囲
１０３　部分帯域１の部分アレイベクトル（ｘ_{ＳｓｕｂＰ}）
１１０１　主対角成分
１１０２、１１０３　第一副対角成分
１１０４、１１０５　Ｑ_ｉよりｑの値が大きい位置にある副対角成分の要素
１２０１　可干渉性損失関数Ｃ_１（ｔ_１’（τ_１）；ｑ）
１２０２　可干渉性損失関数Ｃ_Ｐ（ｔ_Ｐ’（τ_１）；ｑ）
１３０１、１３０２、１３０３　有効干渉長の周波数依存関数Ｂ_λの標本点
１４０１、１４０２、１４０３　Ｃ_１（ｔ_１’,ｑ）の時間変化
１５０１、１５０２、１５０３　Ｃ_ｉ（ｔ_ｉ’,ｑ）の周波数変化
１６０１、１６０２　ｑ’の計算基準となる中心線の対
１６０３～１６１０　ｑ’の値に対応する要素位置
１７０１　口径制限関数Ｄ_１（ｔ_１’，ｑ’）
１７０２　口径制限関数Ｄ_Ｐ（ｔ_Ｐ’，ｑ’）
１８０１、１８０２、１８０３　有効口径長の周波数依存関数Ｂ_αの標本点
１９０１、１９０２、１９０３　Ｄ_１（ｔ_１’,ｑ’）の時間変化
２００１、２００２、２００３　Ｄ_ｉ（ｔ_ｉ’,ｑ’）の周波数変化
Ａ_１、Ａ_Ｐ　空間アンサンブル平均数（標本数）
Ｂ_ｉ　可干渉性損失行列
Ｂ_λ　ｔにおける有効干渉長の周波数依存関数
Ｂ_α　ｔにおける有効口径長の周波数依存関数
Ｃ_１、Ｃ_ｉ、Ｃ_Ｐ　可干渉性損失関数
Ｄ_１、Ｄ_ｉ、Ｄ_Ｐ　口径制限関数
Ｈ　行列の行番号または列番号
Ｉ　行列の行番号
Ｉ’単位行列
Ｊ　行列の列番号
ｆ_ｓｕｂ１、ｆ_ｓｕｂ２、ｆ_ｓｕｂｉ、ｆ_ｓｕｂＰ　部分帯域代表周波数
ｑ　対角位置変数
ｑ’口径位置変数
α_１、α_ｉ、α_Ｐ　有効口径長
λ_１、λ_ｉ、λ_Ｐ　有効干渉長
２１０１　計算モデル
２１０２　第一媒質領域
２１０３　第二媒質領域
２１０４　点音源
２１０５　受信面配列モデル
２２０１　音圧波形
２３０１　周波数スペクトル
２７０１、２７０２、２９０１、２９０２、３１０１、３１０２　尖頭位置

Claims

　被検体からの超音波エコーを受信する複数の受信素子と、所定の受信焦点の位置に応じて前記複数の受信素子の受信信号をそれぞれ遅延させて第一の時間信号群を形成する遅延部と、前記第一の時間信号群を重み付け処理する適応処理部とを有し、
　前記適応処理部は、複数の部分帯域適応処理部を有し、
　複数の前記部分帯域適応処理部はそれぞれ、前記第一の時間信号群から所定の周波数帯域の第二の時間信号群を生成し、前記第二の時間信号群に適応処理を行って適応重みを推定し、前記推定した適応重みにより前記第二の時間信号群を重み付けして加算したものを前記部分帯域適応処理部の出力とし、複数の前記部分帯域適応処理部の出力を加算することを特徴とする超音波撮像装置。
　請求項１に記載の超音波撮像装置において、複数の前記部分帯域適応処理部のうち１以上は、生成する前記周波数帯域の一部または全部が、他の前記部分帯域適応処理部が生成する周波数帯域と重なっていることを特徴とする超音波撮像装置。
　請求項１または２に記載の超音波撮像装置において、前記複数の部分帯域適応処理部は、相互に独立して前記適応処理を行うことを特徴とする超音波撮像装置。
　請求項１ないし３のいずれか１項に記載の超音波撮像装置において、前記複数の部分帯域適応処理部は、相互に独立して、前記第二の時間信号群をそれぞれ時間軸方向および空間軸方向の少なくとも一方について再標本化した後、前記適応処理を行うことを特徴とする超音波撮像装置。
　請求項１ないし４のいずれか１項に記載の超音波撮像装置において、前記複数の部分帯域適応処理部は、相互に独立して、前記第二の時間信号群について、時間軸方向および空間軸方向の少なくとも一方についてアンサンブル平均をとった後、前記適応処理を行うことを特徴とする超音波撮像装置。
　請求項１ないし４のいずれか１項に記載の超音波撮像装置において、前記複数の部分帯域適応処理部は、相互に独立して、前記第二の時間信号群から相関行列あるいは共分散行列を求め、前記相関行列あるいは共分散行列の行列要素を重み付けした後、前記適応処理を行うことを特徴とする超音波撮像装置。
　請求項６に記載の超音波撮像装置において、前記相関行列あるいは共分散行列の行列要素の重み付けの処理を、前記相関行列あるいは共分散行列と、所定の可干渉性損失行列との行列要素ごとの積を求めることにより行うことを特徴とする超音波撮像装置。
　請求項７に記載の超音波撮像装置において、前記可干渉性損失行列の行列要素の値は、非負であり、主対角成分が最も大きく、前記主対角成分から離れる副対角成分であるほど小さくなり、前記主対角成分から所定値以上離れた副対角成分上の行列要素はゼロであることを特徴とする超音波撮像装置。
　請求項５に記載の超音波撮像装置において、前記複数の部分帯域適応処理部は、それぞれ独立して、前記アンサンブル平均の条件を前記第一の時間信号群の受信開始からの経過時間とともに変化させることを特徴とする超音波撮像装置。
　請求項８に記載の超音波撮像装置において、前記複数の部分帯域適応処理部は、それぞれ独立して、前記第一の時間信号群の受信開始からの経過時間、周波数、および、前記可干渉性損失行列の行列要素の位置に依存する関数、の計算手段を有し、前記計算手段は、前記受信開始からの時間、前記複数の部分帯域適応処理部を代表する周波数、行列要素の位置に応じた前記可干渉性損失行列を求めることを特徴とする超音波撮像装置。
　請求項１０に記載の超音波撮像装置において、前記関数の入力手段を備え、前記可干渉性損失行列に関する情報を表示する手段を有することを特徴とする超音波撮像装置。
[規則91に基づく訂正 16.07.2013]　
　請求項１に記載の超音波撮像装置において、複数の前記部分帯域適応処理部の出力に、それぞれ個別に、時間に依存する異なる重みを乗じて、コヒーレント加算あるいはインコヒーレント加算する加算部をさらに有することを特徴とする超音波撮像装置。
　請求項５に記載の超音波撮像装置において、複数の前記部分帯域適応処理部はそれぞれ独立した条件で、前記第一の時間信号群の受信開始からの経過時間とともに、前記第二の時間信号群から生成する前記相関行列あるいは共分散行列の次数を低減し、前記第二の時間信号群について、時間軸方向および空間軸方向の少なくとも一方について求めるアンサンブル平均数を増加させることを特徴とする超音波撮像装置。