JPH0425017B2 - - Google Patents

Description

【発明の詳細な説明】 発明の背景 Ａ 技術分野 本発明は、超音波を物体に送信し、物体の内部
からの反射超音波を受信して、物体内部の音響特
性を測定する超音波測定方法およびその装置の改
良に係り、特に物体内部の超音波伝播に伴う減衰
に関する情報をより正確に提供する超音波測定方
法およびその装置に関する。 Ｂ 先行技術とその問題点 超音波測定技術は現在金属探傷、魚群探知、医
療診断分野等広範囲にわたつて利用されている。
中でも医療用の超音波断層装置の最近の発展には
目をみなるものがある。 超音波断層装置は原理的には、パルスエコー法
を用いており、被測定物体としての生体内へ送信
された超音波パルスが生体内部の音響インピーダ
ンスの異なる境界で反射する現象を利用し、この
反射波（エコー）を受信していわゆるＢモード法
による生体の断層像を表示するものである。従つ
てこのエコーには生体内部での超音波の減衰、音
響インピーダンス、音速等の様様な情報が含まれ
ている。しかし従来の装置では、これらの各情報
が明確に分離されず、単にエコーの振幅を表示し
ているにすがない。 具体的には、生体内の音速を一定と仮定し、さ
らに生体内の超音波伝播による減衰は、いわゆる
STC（Sensitivity Time Control）回路と呼ばれ
ている回転によつて任意的に補正をしたエコー振
幅値を輝度に変調し、これをブラウン管上に断層
像として表示しているにすぎない。従つて、得ら
れた断層像は、生体内部の音響インピーダンス境
界面の２次元的分布を定性的に画像化したものと
なり、必然的に生体組織の位置や形に関する形態
情報がその利用の中心となつている。つまり、生
体組織の特性である減衰度、音速等の測定はなさ
れていないのが現状である。 生体組織の減衰情報を得ようとする試みがいく
つか報告されている。後で詳しく述べるように、
エコー信号には生体組織伝播による減衰と音響イ
ンピーダンスの異なる境界での反射強度の２つの
情報が含まれており、両者はいずれも、未知であ
る。したがつて厳密にこの２つの影響を分離する
ことは、今のところ極めて困難であると言わざる
を得ない。 反射強度が超音波の周波数に依存しないと仮定
した場合には、被測定物体の同一部分について複
数の異なる周波数の超音波を送受信しエコーの各
周波数成分の音圧比を測定すれば、反射強度の影
響をなくして伝播による減衰係数を求めることが
可能となる。このような仮定は、超音波の波長に
比べて十分大きな広がりをもつ音響境界、例えば
平面反射板の場合に成立する。しかし実際の生体
組織では、使用する超音波の波長程度あるいはそ
れ以下の大きさの散乱体が存在し得るので、この
仮定は生体組織全体を考えたとき必ずしも成立す
ることは考えにくい。 このような状況を鑑みて本出願人は、生体の超
音波伝播による減衰情報を音響境界での反射強度
の影響を少なくして測定する方法および装置を特
願昭55−49571として出願している。 更に前記出願を改良し、複数の異なる周波数
で、被測定物体からのエコー強度を測定し、これ
によつて被測定物体の減衰係数および減衰係数の
周波数依存性を近似的に測定することが出来る方
法および装置を特願昭58−229853として出願して
いる。しかしながら後者の出願における発明で求
められた減衰係数およびその周波数依存性は、被
測定物体の反射強度の周波数依存性が測定する周
波数の範囲で一定であればその誤差を少なく出来
るが、一定でない場合には反射強度による誤差を
少なく抑えることが難しくなることが予想され
る。 発明の目的 本発明の目的は、上記先行技術の問題点を解決
し、被測定物体の減衰係数および減衰係数の周波
数依存性を、この物体の反射強度による影響をよ
り少なくすることで、より正確に測定出来る超音
波測定方法およびその装置を提供することにあ
る。 本発明によれば、複数の異なる周波数の超音波
パルスを被測定物体に送信し、物体内より反射さ
れた超音波パルスのエコーを検出し、この検出さ
れた複数の周波数のエコーを情報処理することに
よつて物体の超音波特性に関する定量的な情報を
得る超音波測定方法において、２つの異なる周波
数より測定された減衰係数の平均値、およびれら
２つの周波数を含む３つの異なる周波数より測定
された減衰係数に関連した値と、その周波数の関
数とから、物体の減衰係数に関する情報を算出す
る。 本発明の１つの態様によれば、減衰係数に関す
る情報は、被測定物体内における超音波の減衰係
数を含む。 本発明の他の態様によれば、減衰係数に関する
情報は、被測定物体内における超音波の減衰係数
の周波数依存性を含む。 本発明によればまた、複数の異なる周波数の超
音波パルスを被測定物体に送信する手段と、物体
内より反射された超音波パルスのエコーを検出す
る手段と、検出された複数の周波数のエコーを情
報処理し物体の超音波特性に関する定量的な情報
を得る手段と、得られた定量的な情報を表示する
手段とを有する超音波測定装置において、情報を
得る手段は、２つの異なる周波数の減衰係数の平
均値、およびこれら２つの周波数を含む３つの異
なる周波数の減衰係数に関連した値と、その周波
数の関数を求め、これらの平均値、関連した値お
よびその周波数の関数値より物体の減衰係数に関
する情報を算出する。 発明の具体的説明および作用 次に添付図面を参照して、本発明による超音波
測定方法およびその装置の実施例を詳細に説明す
る。 第２図において超音波探触子１から被測定物体
４に放射された超音波パルス１０は、音響特性不
連続面２および３によつてそれぞれ反射され、そ
れぞれの反射パルス２０および３０が探触子１に
よつて検出される。 いま、パルス放射から反射パルス受信までの時
間をｔとすると、超音波探触子１から音響特性不
連続面までの距離ｘは、音速をＶとして次式で与
えられる。 ｘ＝Ｖ・ｔ／２ ここで周波数で放射されたパルス１０の強度
をI₀（）、受信した反射パルスの強度をＩ（、
ｘ）とすると、近似的に次式が成立する。 Ｉ（，ｘ）＝I₀（）・^{a( ,x)}・ｇ（ｘ）・x^b(
x)・exp｛−4∫^x ₀α（，ｘ）dx｝……(1) (1)式の両辺の自然対数をとると、 lnI（，ｘ）＝lnI₀（）＋ａ（，ｘ）・ln＋l
n〔ｇ（ｘ）・x^b(x)〕−4∫^x ₀α（，ｘ）dx……(2) ここでα（、ｘ）は被測定物体４の超音波伝
播による減衰係数であり、^{a( ,x)}・ｇ（ｘ）は、周
波数依存性を考慮に入れた音響特性不連続面の反
射強度である。 ａ（、ｘ）の値は、波長（λ＝Ｖ／）より
充分大きい音響特性不連続面では、ａ（、ｘ）＝
０であり、波長より充分小さい音響特性不連続面
では、ａ（、ｘ）＝４である。従つてある特定の
周波数帯域内ではａ（、ｘ）は一定であり、生
体については０ａ（、ｘ）４であると考え
られる。またx^b(x)は、反射波が広がることに起因
して超音波探触子１の位置での反射強度が弱まる
結果を考慮したもので、充分に広い音響特性不連
続面では、ｂ（ｘ）＝０、小さい音響特性不連続面
では、ｂ（ｘ）＝−２である。従つて一般に−２
ｂ（ｘ）０である。但し、以上の考察では、放
射超音波１０は理想的に細いペンシル・ビームと
し、また近接した複数の音響特性不連続面からの
反射波相互間の干渉の効果は無視している。 エコー強度の分布から情報処理をして被検査物
体４の超音波特性に関する定量的情報を得る原理
は、(2)式から出発する。 数学的見通しを良くするため、ここで仮りにエ
コー強度Ｉ（、ｘ）が距離ｘおよび周波数に
ついて連続的に得られるとし、(2)式の両辺を、ま
ず・lnで一階偏微分すると、 ∂lnI（，ｘ）／∂（ln）＝∂lnI₀（）／∂
（ln）＋∂a（，ｘ）／∂（ln）ln＋ａ（
，ｘ）／（１＋n）−4∫^x ₀∂α（，ｘ）／∂（
ln）dx……(3) となり、さらに(3)式をｘで微分して整理すると、 ∂a（，ｘ）／∂（ln）＝−１／４∂²lnI（
，ｘ）／∂x∂（ln）＋J₂……(4) J₂＝１／４∂²a（，ｘ）／∂x∂（ln
）ln＋１／４１／（１＋ln）∂a（，ｘ）／∂x
が得られる。 さらに(4)式を・lnでさらに微分すると、 ∂²α（，ｘ）／∂（ln）²＝−１／４∂³lnI（
，ｘ）／∂x∂（ln）²＋J₃ ……(5) J₃＝∂J₂／∂（ln）＝１／４∂³a（，ｘ）／∂
x∂（ln）²ln＋１／２１／（１＋ln）∂²a（
，ｘ）／∂x∂（ln） ＋１／４ ２＋ln／²（１＋ln）³∂a（，ｘ
）／∂x 被検査物体の超音波減衰係数が周波数の幅β
（ｘ）乗に比例すると α（、ｘ）＝α₀（ｘ）・〓^(x) ……(6) となり、従つて(4)式の左辺は次のようになる。 ∂α₀（ｘ）〓^(x)／∂（ln）＝∂α₀（ｘ）
〓^(x)／∂∂／∂（ln）＝β（ｘ）α（，ｘ
）／／１＋ln＝β（ｘ）α（，ｘ）／（１＋ln
）……(7) (5)式の左辺は次のようになる。 ∂²α（，ｘ）／∂（ln）²＝∂／∂（ln）
∂α（，ｘ）／∂（ln）＝∂／∂（ln）β（
ｘ）α（，ｘ）／（１＋ln） ＝β（ｘ）α（，ｘ）〔（β（ｘ）−１）（１＋
ln）−１〕／²（１＋ln）³……(8) (8)／(7)より β（ｘ）＝（１＋ln）〔∂²α（，ｘ）／∂
（ln）²／∂α（，ｘ）／∂（ln）〕＋２＋l
n／１＋ln……(9) (9)式と(7)式より α（，ｘ）＝∂α（，ｘ）／∂（ln）／
〔∂α（，ｘ）／∂（ln）²／∂α（，ｘ）／
∂（ln）〕＋２＋ln／（１＋ln）²……(10) となり、減衰係数α（、ｘ）およびこの周波数
依存性β（ｘ）が求まる。すなわち(9)式および(4)、
(5)式より

【表】
〓
又(10)式および(4)、(5)式より α（，ｘ）＝〔−１／４∂²lnI（，ｘ）／∂x∂
（ln）＋J₂〕／〔−１／４∂³lnI（，ｘ）／∂x
∂（ln）²＋J₃〕／〔−１／４∂²lnI（，ｘ）／
∂x∂（ln）＋J₂〕＋２＋ln／（１＋ln）²…
…（12） 又(9)式のβ（ｘ）、(10)式のα（、ｘ）よりα₀
（ｘ）が、α₀（ｘ）＝α（、ｘ）／〓^(x)より求め
ら
れる。 （11）式の分母および分子の第１項は測定量で
あるが、第２項以降は非測定量である音響特性不
連続面の反射強度による項であり、β（ｘ）を求
める場合には誤差を生じる原因となる。 反射強度の周波数依存性を示す係数ａ（、ｘ）
が測定する周波数の範囲で一定であれば、すなわ
ち∂a（，ｘ）／∂＝０であれば、ａはｘのみの関
数 ととなり、従つて(4)式および(5)式はそれぞれ次式
のようになる。 ∂α（，ｘ）／∂（ln）＝−１／４∂²lnI
（，ｘ）／∂x∂（ln）＋１／４１／（１＋ln
）∂a（，ｘ）／∂x……（13） ∂²α（，ｘ）／∂（ln）²＝−１／４∂³l
nI（，ｘ）／∂x∂（ln）²＋１／４２＋ln／
²（１＋ln）³∂a（，ｘ）／∂x……（14） 従つて（11）式は（13）および（14）式より、 β（ｘ）＝（１＋ln）〔−１／４∂³lnI（，ｘ
）／∂v∂（ln）²＋１／４２＋ln／²（１＋ln
）³∂a（，ｘ）／∂_x〕／〔−１／４∂²lnI（，
ｘ）／∂x∂（ln）＋１／４ １／（１＋ln）
∂a（，ｘ）／∂_x〕……（15） となり、同様に（12）式は次のようになる。 α（，ｘ）＝ 〔−１／４∂²lnI（，ｘ）／∂x∂（ln）＋１／
４１／（１＋ln）∂a（，ｘ）／∂_x〕／〔−１／
４∂³lnI（，ｘ）／∂x∂（ln）²＋１／４２＋ln
／²（１＋ln）³∂a（，ｘ）／∂_x〕／〔１／４
∂²lnI（，ｘ）／∂x∂（ln）＋１／４1∂a（
，ｘ）／（ln）∂_x〕※ ※ ／＋（２＋ln／（１＋ln）²） ……（16） これらは（11）式および（10）式に比べて誤差
の項が少なくなつていることがわかる。 更に標的強度ａ（ｘ）が距離ｘに依存せず一定
であれば、∂a（ｘ）／∂x＝０より（15）式は次式
のように誤差を含まない式となる。 β（ｘ）＝（１＋ln）〔∂³lnI（，
ｘ）／∂x∂（ln）²／∂²lnI（，ｘ）／∂x∂（
ln）〕……（17） つまり、反射鏡度の周波数依存性を示す係数ａ
（、ｘ）が測定する範囲の周波数および距離
ｘに依存せず一定であれば、^aとなり、（17）式
によつて減衰係数α（、ｘ）の周波数依存性β
（ｘ）が正確に測定されるわけである。この場合
には（16）式より、α（、ｘ）も次式のように
正確に測定することが可能となる。 α（，ｘ）＝〔−１／４∂²lnI（，ｘ）／∂x∂
（ln）〕／〔−１／４∂³lnI（，ｘ）／∂x∂（
ln）²〕／〔−１／４∂²lnI（，ｘ）／∂x∂（
ln）＋〔２＋ln／（１＋ln）²〕……（18） 以上のようにａ（、ｘ）が、特定の条件下で
ない場合にも、本発明によればβ（ｘ）、α（、
ｘ）をより正確に求めることが出来る。 具体的に説明すると、ａ（、ｘ）＝a₀（ｘ）
で与えられる場合、つまり測定周波数の範囲でａ
（、ｘ）が、周波数に比例して変化する場合に
(4)式および(5)式はそれぞれ次式のようになる。 ∂α（，ｘ）／α（ln）＝−１／４
∂²lnI（，ｘ）／∂x∂（ln）＋１／４∂a₀（ｘ
）／∂x……（19） ∂²α（，ｘ）／α（ln）²＝−１／４∂³lnI（
，ｘ）／∂x∂（ln）² ……（20） 従つて（11）式は（19）および（20）式より β（ｘ）＝（１＋ln）〔−１／４∂³lnI（，ｘ
）／∂x∂（ln）²／−１／４∂²lnI（，ｘ）／∂
x∂（ln）＋１／４∂a₀（ｘ）／∂x〕＋２＋ln／
１＋ln……（21） となり、同様に(10)式は次のようになる。 α（，ｘ）＝−１／４∂²lnI（，ｘ）／∂x∂（
ln）＋１／４∂a₀（ｘ）／∂x／〔−１／４∂³lnI
（，ｘ）／∂x∂（ln）²／−１／４∂²lnI（，
ｘ）／∂x∂（ln）＋１／４∂a₀（ｘ）／∂x〕＋２
＋ln／（１＋ln）²
……（22） これは（11）式および（12）式に比べて誤差の
項が少なくなつていることがわかる。すなわち
（21）式では誤差の分母の１／４・∂a₀（ｘ）／∂xの
みで あり、（22）式では分母および分子の１／４ ∂a₀（ｘ）／∂xのみである。 以上の説明ではエコー信号強度Ｉ（、ｘ）が
距離ｘおよび周波数に関して連続的に得られる
場合について述べたが、実際の測定の場合には、
エコー信号は散乱体のある位置ｘにより離散的に
与えられるので、若干の修正を必要とする。 すなわち、(4)式および(5)式で用いた微分の代わ
りに差分を用いなければならない。第２図に示す
ように、散乱体２，３が離散的に存在する位置を
x₁およびx₂として両位置の間で(2)式の差分を求め
変形すると、次式のようになる。 ∫^x2 _x1α（，ｘ）dx＝−１／４ln〔Ｉ（，x₂）／
Ｉ（，x₁）〕＋１／４〔ａ（，x₂）−ａ（，x₁）
〕ln ＋１／４ln〔ｇ（x₂）・x₂ ^b(x2)／ｇ（x₁）・x₁ ^b(x
1)〕……（23） (4)式に対応する式は次のようになる。 K₂（₁，₃）≡１／x₂−x₁∫^x2 _x1α（₁，ｘ）−
α（₃，ｘ）／₁ln₁−₃ln₃dx ＝−１／４（x₂−x₁）（₁ln₁−₃ln₃）ln〔
Ｉ（₁，x₂）／Ｉ（₁，x₁）／Ｉ（₃，x₂）／Ｉ（
₃，x₁）〕＋L₂……（24） L₂（₁，₃）＝１／４｛〔ａ（₁，x₂）−ａ（₁
，x₁）〕ln₁−〔ａ（₃，x₂）−ａ（₃，x₁）〕ln
₃｝／（x₂−x₁）（₁ln₁−₃ln₃）……（25） 又(5)式に対応する式は次のようになる。 K₃（₁，₂，₃）≡K₂（₁，₂）−K₂（₂，
₃）／₁ln₁＋₂ln₂／２−₂ln₂＋₃ln₃／
２ ＝−１／２（x₂−x₁）（₁ln₁−₂ln₂）（₂l
n₂−₃ln₃）（₁ln₁−₃ln₃） ｛ln〔Ｉ（₁，x₂）／（₁，x₁）〕（₂ln₂−
₃ln₃）＋ln〔Ｉ（₂，x₂）／（₂，x₁）〕（₃ln
₃−₁ln₁） ＋ln〔Ｉ（₃，x₂）／（₃，x₁）〕（₁ln₁−
₂ln₂）｝＋L₃（₁，₂，₃）……（26） L₃（₁，₂，₃）＝２｛L₂（₁，₂）−L₂（₂，
₃）｝／₁ln₁−₃ln₃ ＝１／２（〔ａ（₁，x₂）−ａ（₁，x₁）〕ln₁
（₂ln₂−₃ln₃）／（x₂−x₁）（₁ln₁−₂l
n₂）（₂ln₂−₃ln₃）（₁ln₁−₃ln₃）
＊ ＊＋〔ａ（₂，x₂）−ａ（₂，x₁）〕ln₂（₃ln
₃−₁ln₁）／ ＊＊ ＊＊〔ａ（₃，x₂）−ａ（₃，x₁）〕ln₃（₁ln
₁−₂ln₂）／ ……（27） α（、ｘ）＝α₀（ｘ）〓^(x)より（24）式を左辺
を(7)式に次式のように、近似する。 １／x₂−x₁∫^x2 _x1α（₁，ｘ）−α（₃，ｘ）／₁l
n₁−₃ln₃dx≒１／x₂−x₁∫^x2 _x1β（ｘ）α（₃
＋₁／２，ｘ）dx／（₁＋₃／２）・（１＋ln₁＋
₃／２）……（28） 同様（26）式の左辺を(8)式より次式のように近
似する。 １／x₂−x₁∫^x2 _x1（α（₁，ｘ）−α（₂，ｘ）／
₁ln₁−₂ln₂−α（₂，ｘ）−α（₃，ｘ）／
₂ln₂−₃ln₃）２／₁−₃dx ≒１／x₂−x₁∫^x2 _x1β（ｘ）α（₁＋₂＋₃／３
ｘ）〔（β（ｘ）−１）（１＋ln₁＋₂＋₃／３）
−１〕dx／（₁＋₂＋₃／３）²・（１＋ln₁＋₂
＋₃／３）³……（29） （28）式に₁＋₃／２（１＋ln₁＋₃／２）をか
けて 整理し、（24）、（25）式とまとめると次式のよう
になる。 K₂′（₁，₃）≡（₁−₃／２）（１＋ln₁＋
₃／２）K₂ ≒１／x₂−x₁）∫^x2 _x1β（ｘ）α（₁＋₃／２，ｘ）
dx＝₁ ＋₃／８（x₂−x₁）（₁ln₁−₃ln₃）（１＋
ln₁−₃／２）× ln〔Ｉ（₁，x₁）／Ｉ（₁，x₂）／Ｉ（₃，x₁）／
Ｉ（₃，x₂）〕＋L₂′L₂′＝₁ ＋₃／８（x₂−x₁）（₁ln₁−₃ln₃）（１＋
ln₁＋₃／２）× ｛〔ａ（₁，x₂）−ａ（₁，x₁）〕ln₁−〔ａ（₃
，x₂）−
ａ（₃，x₁）〕ln₃｝ ……（30） （29）式に（₁＋₂＋₃／３）²・（１＋ln₁ ＋₂＋₃／３）³をかけて整理し、（26）、（27）
式と まとめると次式のようになる。 K₃′（₁，₂，₃）＝（₁＋₂＋₃／３）²（１
＋ln₁ ＋₂＋₃／３）³K₃（₁，₂，₃）≒１／x₂−x
₁∫^x2 _x1β （ｘ）α（，ｘ）〔（β（ｘ）−１）（１＋ln）−
１〕
dx＝−Ｂ／４（x₂−x₁）｛ln〔Ｉ（₁，x₁）／Ｉ（₁
，x₂）〕（₂ln₂−₃ ln₃）＋ln〔Ｉ（₂，x₁）／Ｉ（₂，x₂）〕（
₃ln₃−₁ln₁）＋ln 〔Ｉ（₃，x₁）／Ｉ（₃，x₂）〕×（₁ln₁−₂l
n₂）｝（₁＋₂＋₃／３）² （１＋ln₁＋₂＋₃／３）³＋L₃′L₃′＝Ｂ／４（x
₂−x₁） ｛〔ａ（₁，x₂）−ａ（₁，x₁）〕ln₁（₂ln₂
−₃ln₃）
＋〔ａ（₂，x₂）−ａ（₂，x₁）〕ln₂）〕ln₂（
₃ln₃
−₁ln₁）＋〔ａ（₃，x₂）−ａ（₃，x₁）〕ln
₃
（₁ln₁−₂ln₂）｝×（₁＋₂＋₃／３）²（
１＋ln₁ ＋₂＋₃／３）³ Ｂ＝２／（₁ln₁−₂ln₂）（₂ln₂−₃ln₃
）（₁ln1−3ln₃）……（31） ここでβ（ｘ）を求めるための(9)式に対応する
式は次式のようになる。 β（ｘ）＝₁＋₂＋₃／３・｛１＋ln（₁＋₂＋
₃／３）｝・〔K₃′／K₂′〕＋２＋ln（₁＋₂＋₃
／３）／１＋ln（₁＋₂＋₃／３）……（32） 又、α（、ｘ）を求めるための(10)式に対応す
る式は次式のようになる。 α（，ｘ）＝K₂′／〔K₃′／K₂′〕＋２＋ln（₁＋
₂＋₃／３）／₁＋₂＋₃／３・｛１＋ln₁＋
₂＋₃／３｝²……（33） ここで、前に述べたと同様に、反射強度の周波
数依存性^a(、^x)が測定する周波数の範囲で一定
であれば、すなわちａ（₁、x₁）＝ａ（₂、x₁）＝ａ
（₃、x₁）≡ａ（x₁）およびａ（₁、x₂）＝ａ（₂
、
x₂）＝ａ（₃、x₂）≡ａ（x₂）であればａはｘのみ
の関数となり、よつて（30）式および（31）式は
次式のようになる。 K₂′≒１／x₂−x₁∫^x2 _x1β（ｘ）・α（₁＋₃／２
，ｘ）dx ＝₁＋₃／８（x₂−x₁）（₁ln₂−₁ln₃）・（
１＋ln₁ ＋₃／２）・×ln〔Ｉ（₁，x₁）／Ｉ（₁，x₂）
／Ｉ（₃，x₁）／Ｉ（₃，x₂）〕＋₁ ＋₃／８（x₂−x₁）（₁ln1−3ln₃）×（１
＋ln₁ ＋₃／２）・〔（ａ（x₂）−ａ（x₁））（ln₁／
₃）〕≡M₂′ ……（34） K₃′≒１／x₂−x₁∫^x2 _x1β（ｘ）α（₁＋₂＋₃
／３，ｘ） 〔（β（ｘ）−１（１＋ln₁＋₂＋₃）／３）−１
〕dx −Ｂ／４（x₂−x₁）｛ln〔Ｉ（₁，x₁）／Ｉ（₁，x₂
）〕（₂ln₂−₃ln₃） ＋ln〔Ｉ（₂，x₁）／Ｉ（₂，x₂）〕（₃ln₃−
₁ln₁）＋ln 〔Ｉ（₃，x₁）／Ｉ（₃，x₂）〕（₁ln₁−₂ln
₂）｝（₁＋₂＋₃／３）² （１＋ln₁＋₂＋₃／３）³＋Ｂ／４（x₂−x₁）〔ａ
（x₂）−ａ （x₁）〕・｛₁ln₁（₃／₂）＋₂ln₂・ln（
₁／₃）＋
₃ln₃ln（₂／₁）｝×（₁＋₂＋₃／３）²（
１＋ln₁ ＋₂＋₃／３）³≡M₃′ ……（35） （32）式と（34）式および（35）式によつてβ
（ｘ）は次式で求められる。 β（ｘ）＝₁＋₂＋₃／３｛１＋ln₁＋₂＋₃／
３｝〔M₃′／M₂′〕＋２＋ln（₁＋₂＋₃／３）／
１＋ln（₁＋₂＋₃／３）……（36） 又α（、ｘ）は（33）、（34）および（35）式
によつて、次式で求められる。 α（，ｘ）＝M₂′／〔M₃′／M₂′〕＋２＋ln（₁＋
₂＋₃／３）／₁＋₂＋₃／３｛１＋ln₁＋₂
＋₃／３｝²……（37） となり、これら２つの式は（32）式および（33）
式に比べると、誤差の項が少なくなつていること
がわかる。 又ａ（ｘ）が、距離ｘに依存せず一定であれば、
ａ（x₂）＝ａ（x₁）より（36）式は次式のように、
さらに誤差の少ない式となる。 β（ｘ）＝₁＋₂＋₃／３｛１＋ln₁＋₂＋₃／
３｝〔N₃′／N₂′〕＋２＋ln（₁＋₂＋₃／３）／
１＋ln（₁＋₂＋₃／３） N₂′＝₁＋₃／８（x₂−x₁）（₁ln₁−₃ln₃）
（１＋ln₁＋₃／２）・×ln〔Ｉ（₁，x₁）／Ｉ（
₁，x₂）／Ｉ（₃，x₁）／Ｉ（₃，x₂）〕 N₃′＝−Ｂ／４（x₂−x₁）｛ln〔Ｉ（₁，x₁）／Ｉ（
₁，x₂）〕（₂ln₂−₃ln₃）＋ln〔Ｉ（₂，x₁
）／Ｉ（₂，x₂）〕（₃ln₃−₁ln₁） ＋ln〔Ｉ（₃，x₁）／Ｉ（₃，x₂）〕（₁ln₁−
₂ln₂）｝（₁＋₂＋₃／３）²（１＋ln₁＋₂
＋₃／３）³……（38） α（、ｘ）も次式より、より誤差を少なくし
て求められる。 α（、ｘ）＝N₂′／〔N₃′／N₂′〕＋２＋ln（₁＋
₂＋₃／３）／₁＋₂＋₃／３｛１＋ln（₁＋
₂＋₃／３）｝² 以上のように、ａ（、ｘ）が、特定の条件下
でない場合にも本発明ではβ（ｘ）、α（、ｘ）
をより正確に求めることが出来る。すなわちａ
（、ｘ）＝a₀（ｘ）・で与えられる場合、つまり
測定周波数の範囲で、ａ（、ｘ）が周波数に比
例して変化する場合に、（30）式、（31）式は次式
のようになる。 K₂′＝₁＋₃／８（x₂−x₁）（₁ln₂−₃ln₃
）（１＋ln₁＋₃／２）・×ln〔Ｉ（₁，x₁）／Ｉ
（₁，x₂）／Ｉ（₃，x₁）／Ｉ（₃，x₂）〕 ₁＋₃／８（x₂−x₁）（１＋ln₁＋₃／２）×
〔a₀（x₂）−a₀（x₁）〕≡P₂′……（40） K₃′＝−Ｂ／４（x₂−x₁）｛ln〔Ｉ（₁，x₁）／Ｉ
（₁，x₂）〕（₂ln₂−₃ln₃）＋ln〔Ｉ（₂，
x₁）／Ｉ（₂，x₂）〕 （₃ln₃−₁ln₁）＋ln〔Ｉ（₃，x₁）／Ｉ
（₃，x₂）〕（₁ln₁−₂ln₂）｝（₁＋₂＋
₃／３）²・（１＋ln₁＋₂＋₃／３）³≡P₃′
……（41） 従つてβ（ｘ）は、（32）式および（40）式よ
り、 β（ｘ）＝₁＋₂＋₃／３｛１＋ln₁＋₂＋₃／
３｝〔P₃′／P₂′〕＋２＋ln（₁＋₂＋₃／３）／
１＋ln（₁＋₂＋₃／３）……（42） で与えられ、又α（、ｘ）は（33）式および
（41）式より、 α（、ｘ）＝P₂′／〔P₃′／P₂′〕＋２＋ln（₁＋
₂＋₃／３）／₁＋₂＋₃／３｛１＋ln（₁＋
₂＋₃／３）｝²……（43） これらは（32）式、（33）式に比べて誤差の項
が少なくなつていることがわかる。 ところで現実には、超音波プローブ１から、例
えば減衰の極めて小さい水中に向けて送信された
超音波ビームは、そのプローブ１の開口、あるい
は中心周波数によつてその音場が変化する。音場
は、近似的には第３図Ａのようになり、また中心
軸上の強度は超音波プローブ１からの距離ｘによ
つて同Ｂに示すように変化する。なお同Ｂの縦軸
は、最大強度I₀に対する距離ｘにおける強度を示
している。 そこで装置の実用上はこれを較正しておかなけ
れば正しい測定が行なえない。すなわち、音圧の
変動をあらかじめ標準媒質で測定し、被測定物体
からのエコー振幅（音圧）を標準音圧で割り算し
て規格化することで、超音波プローブの音場特性
の影響を除き、測定された減衰度を、より普遍的
な値とすることが可能となる。 標準媒質による測定は次のようにして行なうの
が適当である。第４図のように、脱気水１００中
にたとえばステンレスの完全反射体１０２を、設
け、この完全反射体１０２からのエコー振幅を標
準音圧とする。超音波プローブ１と完全反射体１
０２の距離を相対的に変化させ、各距離からのエ
コー振幅を測定すれば、標準音圧曲線が第５図の
ように求められる。 実用に適した装置としては、標準音圧曲線をあ
らかじめ測定し、装置内に記憶させておくことが
有利である。又再度、標準音圧曲線を測定したい
場合には、第６図のように段階状の完全反射体１
０２ａを有する装置を作成すればよい。すなわち
超音波プローブ１を走査機構８で同図の矢印の方
向に水平走査すると、超音波プローブ１と反射体
１０２ａの反斜面との距離が段階的に変化する。
そこで完全反射体１０２ａからのエコー振幅を測
定し、順次記憶すれば標準音圧曲線を得ることが
できる。 以上のように、複数の異なる周波数（₁、₂、
₃）によつて減衰係数α（、ｘ）、α₀（ｘ）およ
びその周波数依存性β（ｘ）を近似的に測定でき
ることが示されたわけである。とくに、反射強度
が観測する範囲の周波数に比例し変化する場合
ａ（、ｘ）＝a₀（ｘ）・によα（、ｘ）、α₀（
ｘ）
およびβ（ｘ）をさらに正確に測定することが可
能である。 次に第１図にブロツク図で示す本発明の実施例
について詳細説明を行う。 第１図に示す実施例は、本発明による超音波測
定方法を実現する装置であり、生体４などの被測
定物体の表面に設定された電気信号と超音波の相
互変換を行なう超音波探触子１に送信回路５が接
続されて送信系が構成され、また受信回路７、対
数増幅回路１８、検波回路１９およびSTC回路
２０によつて受信系が構成されている。探触子１
による走査は走査部８によつて制御回路１７の制
御のもとに行なわれる。 生体４の内部で反射された超音波エコーは表示
部１５に可視像として表示される。表示部１５に
表示されるのは、本発明に従つて算出された、た
とえばエコー振幅の減衰係数および（または）そ
の周波数依存性であるが、これは、受信した超音
波エコーから２つのメモリ回路１０および１２の
使用して演算回路１１によつて得られる。 送信回路５より探触子１に第７図のような急峻
に減衰する広帯域の駆動パルスが印加される。探
触子１は高分子系振動子（ポリフツ化ビニリデ
ン：PVDF）、高分子と無機物の複合系振動子、
あるいは音響整合層を付加したPZT振動子が広
帯域特性をもつ探触子として好ましい。この結果
探触子１より第８図のような広帯域超音波パルス
が被測定物体４の内部へ送波（信）される。 被測定物体４中の音響特性不連続面（たとえば
第２図の２，３）で反射散乱された超音波エコー
は同じ探触子１によつて受波（信）され受信回路
７へ入力される。受信回路７で増幅されたエコー
信号（Ａモード信号）はＡ／Ｄ変換器９でデジタ
ル化され、メモリ１０へ蓄積される。 このメモリ１０内のＡモード信号は、第９Ａ図
および第９Ｂ図に示した所定のアルゴリズムによ
つて演算回路１１で処理され、これによつて前述
のβ（ｘ）、α₀（ｘ）、α（、ｘ）を求め、これら
の値をメモリ１２へ蓄積する。このメモリ１２内
のβ（ｘ）、α₀（ｘ）あるいはα（、ｘ）を必要に
応じて選択し、Ｄ／Ａ変換器１３でアナログ化
し、映像出力増幅回路１４へ入力し、表示部１５
へ可視画像として出力する。 生体４におけるβ（ｘ）、α₀（ｘ）あるいはα
（、ｘ）を測定したい関心領域は、表示部１５
へ展開された生体４のＢモード像２００（第１０
図）上においてたとえば枠２０２で示すように、
関心領域設定回路１６によつて設定される。すな
わち従来行なわれているように、Ｂモード像は、
走査部８によつて探触子１を被測定物体４上に走
査し、Ａモード信号を収集することによつて得ら
れる。走査の方法は、メカニカル・セクタ、コン
パウンドスキヤン、リニア電子スキヤンおよびセ
クタ電子スキヤン等、多くの方式がある。しかし
それらの詳細については本発明に直接関係ないの
で、ここでは説明を省略する。このような走査で
得られたエコー信号は、受信回路７を通り、対数
増幅回路１８によつて対数増幅され、検波回路１
９およびSTC回路２０によつてSTC補正を受け
る。STC補正されたエコー信号は、メモリ２１
に蓄積され、映像出力増幅回路１４によつて表示
部１５へＢモード像として展開される。このＢモ
ード像を出力する各回路の動作は公知であるの
で、説明を省略する。 ところでメモリ２１への書き込み（入力）を停
止し、このメモリ２１内の情報を繰返し出力（再
生すれば）、フリーズＢモード像が表示部１５へ
展開される。このフリーズ像２００上において、
第１０図のように、α（、ｘ）、α₀（ｘ）、β（ｘ
）
を測定したい関心領域２０２を関心領域設定回路
１６で指定する。指定法は、設定回路１６の操作
部（図示せず）を操作して例えば第１０図のよう
に関心領域の区間を枠２０２で指定すればよく、
この方法は、いわゆるキヤリパ計測による２点間
の距離計測法に類似したものである。 次に、演算回路１１で実行されるα（、ｘ）、
α₀（ｘ）、β（ｘ）を測定するアルゴリズムについ
て、第９Ａ図、第９Ｂ図および第９Ｃ図のフロー
チヤートに従つて詳細説明を行う。第９Ａ図、第
９Ｂ図および第９Ｃ図の左側のフローチヤートに
対応して、右側にその処理結果がグラフなどで示
されている。また演算回路１１の機能ブロツクが
第１１図に示されている。 まず加算平均処理部５０では、メモリ１０から
読み出した受信信号を指定区間について加算平均
を行う３００，３０２。次に検波処理部５１で上
記加算平均された受信信号１０００を検波する３
０４。次に高周波成分抽出部５２では、受信信号
１００１の高周波成分のみを抽出した受信信号１
００２を作り、所定レベルL₀以上の受信信号１
００３に変換する３０６。次に極大値検出部５３
では、受信信号１００３の極大値の位置２０００
〜２００３を検出する３０８。 次にエコー強度算出部５４では、検出された極
大値２０００〜２００３のうちで指定区間の始め
と終りの位置に最も近い極大値２０００および２
００３を決定し、この位置に対応した受信信号１
０００の波形について所定幅のハミングウインド
ウをかけて高速フリーエ変換（FFT）を実行し、
Ｉ（、ｘ）を求める３１０。 そこで探触子１による音場特性を較正するため
に、規格部５６では、前述のようにして標準媒質
１００または１００ａなどで求めたＪ（、ｘ）
でＩ（、ｘ）を規格化し、I^（、ｘ）を与える
３１２。生体４の音速C₀を一定、たとえばC₀＝
1530ｍ／ｓと仮定しているので、x₁＝ct₁、x₂＝
ct₂でx₁、x₂が得られる。 次に、前述したように周波数₁、₃における指
定区間での超音波減衰係数の平均値K^₂′及び₁、
₂、₃における同平均値に比例する量K^₃′が減衰
係数算出部５７で求められる３１４，３１６。こ
のK^₂′およびK^₃′よりさらに、目的のβ（ｘ）、α
（、ｘ）およびα₀（ｘ）が次に求められる３１
８，３２０，３２２。これらのK^₂′、K^₃′、β
（ｘ）、α₀（ｘ）およびα（、ｘ）の値がメモリ１
２へ入力され、映像出力増幅器１４を通つて表示
部１５へ出力され、メモリ２１から読み出された
Ｂモード画像と共に表示部１５に併示される。
K^₂′、K^₃′、およびα（、ｘ）は〔Neper／cm〕
の単位で表現されるが、必要に応じて演算回路１
１内で〔dB／cm〕の単位に変換してもよい。又
は、β（ｘ）も表示可能である。α₀（ｘ）も
〔dB／cm・ＭHz〕の単位で表示できる。 ところで本装置による第２の表示モードでは、
指定区間のみの減衰係数を測定するのではなく、
全画面に対して単位画素（区間）当りのK^₂′、
K^₃′、β（ｘ）、α₀（ｘ）およびα（、ｘ）を測定
し、その分布を画像として表示することができ
る。これは、第１２図のように画像の各単位画素
Δsに対して前記と同様のアルゴリズムを適用す
れば、可能となる。この分布像は当然、Ｂモード
像に比べて空間分解能は劣化するが、表示してい
る情報としては、従来全く測定できなかつた減衰
係数に関する情報が得られることが根本的な相違
点となつている。 本実施例では探触子１として広帯域探触子を使
用したが、特開昭56−147082に記載してあるよう
に、複数の異なる周波数帯域を有する探触子を使
用してもよい。 発明の具体的効果 以上のように本発明によれば、複数の異なる周
波数で被測定物体からのエコー強度を測定し、こ
れによつて被測定物体の減衰係数および減衰係数
の周波数依存性をより正確に測定することがで
き、またこれらの値の分布像を得ることができ
る。よつて、従来の形態学的な情報しか得られな
かつた超音波断層測定方法および装置と異なり、
被測定物体の減衰に関する定量的情報を得ること
ができる。

【図面の簡単な説明】

第１図は本発明による超音波測定方法を実現す
る装置の実施例を示すブロツク図、第２図および
第３図は本発明の基本的な原理を説明するための
説明図、第４図、第５図および第６図は本発明の
実施例に使用する超音波探触子の較正を説明する
ための説明図、第７図および第８図は、第１図に
示す実施例の動作説明に使用するパルス波形を示
す波形図、第９Ａ図、第９Ｂ図および第９Ｃ図
は、第１図に示す実施例の装置の動作アルゴリズ
ムを示す説明フロー図、第１０図は、第１図に示
す実施例の装置で表示される超音波段層映像の例
を示す図、第１１図は、第１図に示す実施例の装
置における演算回路の機能的構成を示す機能ブロ
ツク図、第１２図は、第１図に示す実施例の装置
の他の表示モードを説明するための説明図であ
る。 主要部分の符号の説明、１……超音波探触子、
１０，１２，２１……メモリ、１１……演算回
路、１５……表示部、１６……関心領域設定回
路、１８……対数増幅回路、２０……STC回路、
５０……加算平均処理部、５２……高周波成分抽
出部、５３……極大値検出部、５４……エコー強
度算出部、５５……指定区間距離算出部、５６…
…規格化部、５７……減衰係数算出部。

Claims (1)

1. 【特許請求の範囲】 １ 複数の異なる周波数の超音波パルスを被測定
   物体に送信し、該物体内より反射された超音波パ
   ルスのエコーを検出し、該検出された複数の周波
   数のエコーを情報処理することによつて該物体の
   超音波特性に関する定量的な情報を得る超音波測
   定方法において、 (a) 所望の指定区間にて物体内より反射された少
   なくとも３つの周波数f₁、f₂、f₃のエコー信号
   を走査した方位方向にわたつて加算平均する工
   程と、 (b) 前記工程(a)にて加算平均された信号を検波す
   る工程と、 (c) 前記工程(b)にて検波された信号の低周波成分
   を除去するために高周波通過フイルタを通して
   所定レベル以上の信号を抽出する工程と、 (d) 前記工程(c)にて抽出された信号の極大値の位
   置を検出する工程と、 (e) 前記工程(d)にて検出された極大値のうち、前
   記指定区間の始めの位置（x₁）と終りの位置
   （x₂）にそれぞれ最も近い極大値に対応した信
   号の波形について所定幅のハミングウインドウ
   を乗算して高速フーリエ変換を行ない、その信
   号のそれぞれの周波数f₁、f₂、f₃における反射
   強度Ｉ（f₁、x₁）、Ｉ（f₂、x₁）、Ｉ（f₃、x₁）、Ｉ
   （f₁、x₂）、Ｉ（f₂、x₂）、Ｉ（f₃、x₂）を求める工
   程と、 (f) 前記工程(e)にて求めたそれぞれの反射強度Ｉ
   （f₁、x₁）…Ｉ（f₃、x₂）を標準媒質で求めた信
   号の反射強度Ｊ（ｆ、ｘ）にて規格化して規格
   化反射強度I^（f₁，x₁）…I^（f₃，x₂）をそれぞれ
   求める工程と、 (g) 前記工程(f)にて求めた規格化反射強度のう
   ち、２つの周波数f₁、f₃について、次式(1)によ
   り平均値k₂を求める工程と、 k₂＝［（f₁＋f₃）／｛８・（x₂−x₁）・（f₁・1nf₁−
   f₃・1nf₃）｝］Ｘ［１＋1n｛（f₁ ＋f₃）／２｝］X1n［｛I^（f₁，x₁）／I^（f₁，x₂
   ）｝／｛I^（f₃，x₁）／I^（f₃，x₂）｝］……(1) ただし、1nは自然対数 (h) 前記工程(f)にて求めた規格化反射強度のう
   ち、３つの周波数f₁、f₂、f₃について、次式(2)
   により平均値k₃を求める工程と、 k₃＝［−Ｂ／｛４・（x₂−x₁）｝］・｛（f₁＋f₂＋f₃
   ）／３｝²・｛１＋1n（f₁＋f₂＋f₃）｝³ ・［1n｛I^（f₁，x₁）／I^（f₁、x₂）｝・（f₂・1n
   f₂−f₃・1nf₃） ＋1n｛I^（f₃，x₁）／I^（f₃，x₂）｝・（f₁・1nf₁
   −f₂・1nf₂）］……(2) ただし、 Ｂ＝２／｛（f₁・1nf₁−f₂・1nf₂）・（f₂・1nf₂−f₃
   ・1nf₃）・（f₁・1nf₁−f₃・1nf₃）｝ (i) 前記工程(g)、(h)にて求めた平均値k₂、k₃から
   減衰係数α（ｆ、ｘ）、α₀（ｘ）およびその周波
   数依存性β（ｘ）をそれぞれ下式(3)〜(5)にて求
   める工程とを備えたことを特徴とする超音波測
   定方法。 β（ｘ）＝（f₁＋f₂＋f₃）／３・〔［１＋1n（f₁＋f₂
   ＋f₃）／３｝・k₃／k₂ ＋［２＋1n｛（f₁＋f₂＋f₃）／３｝］／［１＋1n（
   f₁＋f₂＋f₃）／３｝］〕……(3) α（ｆ，ｘ）＝k₂〔（k₃／k₂）＋［２＋1n｛（f₁＋f₂
   ＋f₃）／３｝］／［｛（f₁ ＋f₂＋f₃）／３｝・｛１＋1n（f₁＋f₂＋f₃）／３｝
   ｝²〕……(4) α₀（ｘ）＝α（ｆ，ｘ）／f〓^(x) ……(5) ２ 複数の異なる周波数の超音波パルスを被測定
   物体に送信する手段と、該物体内より反射された
   超音波パルスのエコーを検出する手段と、検出さ
   れた複数の周波数のエコーを情報処理し、該物体
   の超音波特性に関する定量的な情報を表示する手
   段とを有する超音波測定装置において、 (a) 所望の指定区間にて物体内より反射された少
   なくとも３つの周波数f₁、f₂、f₃のエコーを走
   査した方位方向にわたつて加算平均する加算平
   均手段と、 (b) 加算平均手段にて加算平均された信号を検波
   する検波手段と、 (c) 検波手段にて検波された信号の低周波成分を
   除去するために高周波通過フイルタを通して所
   定レベル以上の信号を抽出する高周波成分抽出
   手段と、 (d) 高周波成分抽出手段にて抽出された信号の極
   大値の位置を検出する極大値検出手段と、 (e) 極大値検出手段にてケーシング極大値のうち
   指定区間の始め（x₁）と終り（x₂）の位置にそ
   れぞれ最も近い極大値に対応した信号の波形に
   ついて所定幅のハミングウインドウを乗算して
   高速フーリエ変換を行ない、その信号のそれぞ
   れの周波数f₁、f₂、f₃における反射強度Ｉ（f₁、
   x₁）、Ｉ（f₂、x₁）、Ｉ（f₃、x₁）、Ｉ（f₁、x₂）、Ｉ
   （f₂、x₂）、Ｉ（f₃、x₂）を求めるエコー強度算
   出手段と、 (f) エコー強度算出手段にて求めた前記反射強度
   Ｉ（f₁、x₁）…Ｉ（f₃、x₂）を標準媒質で求めた
   反射強度Ｊ（ｆ、ｘ）にてそれぞれ規格化して
   規格化反射強度I^（f₁、x₁）…I^（f₃、x₂）をそれ
   ぞれ求める規格化手段と、 (g) 規格化手段にて求めた規格化反射強度のうち
   ２つの周波数について、平均値k₂を求め、さら
   に３つの平均値k₃を求めて、これら平均値k₂、
   k₃より、減衰係数α（ｆ、ｘ）、α₀（ｘ）および
   周波数依存性β（ｘ）を求める減衰係数算出手
   段とを備えたことを特徴とする超音波測定装
   置。 ３ 前記減衰係数算出手段は、求めた規格化反射
   強度のうち２つの周波数f₁、f₃について、次式(1)
   により平均値k₂を求めることを特徴とする特許請
   求の範囲第２項に記載の超音波測定装置。 k₂＝［（f₁＋f₃）／｛８・（x₂−x₁）・（f₁・1nf₁−
   f₃・1nf₃）｝］・［１＋1n｛f₁ ＋f₃）／２］・1n［｛I^（f₁，x₁）／I^（f₁，x₂）
   ｝／｛I^（f₃，x₁）／I^（f₃，x₂）｝］……(1) ４ 前記減衰係数算出手段は、求めた規格化反射
   強度のうち３つの周波数f₁、f₂、f₃について次式
   (2)により平均値k₃を求めることを特徴とする特許
   請求の範囲第２項に記載の超音波測定装置。 k₃＝｛−Ｂ／４・（x₂−x₁）｝・｛（f₁＋f₂＋f₃）／
   ３｝²・［１＋1n｛（f₁＋f₂ ＋f₃）／３｝³］・［1n｛I^（f₁，x₁）／I^（f₁，x
   ₂）｝・（f₂・1nf₂−f₃・1nf₃） ＋1n｛I^（f₂，x₁）／I^（f₂，x₂）｝・（f₃・1nf₃
   −f₁・1nf₁） ＋1n｛I^（f₃，x₁）／I^（f₃，x₂）｝・（f₁・1nf₁
   −f₂・1nf₂）］……(2) ただし、 Ｂ＝２／（f₁・1nf₁−f₂・1nf₂）・（f₂・1nf₂−f₃・
   1nf₃）・（f₁・1nf₁−f₃・1nf₃） ５ 前記減衰係数演算手段は、求めた平均値k₂、
   k₃から減衰係数α（ｆ、ｘ）、α₀（ｘ）およびその
   周波数依存性β（ｘ）をそれぞれ下式(3)〜(5)にて
   求めることを特徴とする特許請求の範囲第２項に
   記載の超音波測定装置。 β（ｘ）＝［（f₁＋f₂＋f₃）／３］・［１＋1n｛（f₁
   ＋f₂＋f₃）／３｝］・k₃／k₂ ＋［２＋1n｛（f₁＋f₂＋f₃）／３｝］／［１＋1n｛
   （f₁＋f₂＋f₃）／３｝］……(3) α（ｆ，ｘ）＝k₂／〔（k₃／k₂）＋［２＋1n｛（f₁＋
   f₂＋f₃）／３｝］／［｛（f₁ ＋f₂＋f₃）／３｝・｛１＋1n｛（f₁＋f₂＋f₃）／３
   ｝｝²〕……(4) α₀（ｘ）＝α（ｆ，ｘ）／f〓^(x) ……(5)