JP3718856B2 - 物理量解析装置 - Google Patents
物理量解析装置 Download PDFInfo
- Publication number
- JP3718856B2 JP3718856B2 JP30557193A JP30557193A JP3718856B2 JP 3718856 B2 JP3718856 B2 JP 3718856B2 JP 30557193 A JP30557193 A JP 30557193A JP 30557193 A JP30557193 A JP 30557193A JP 3718856 B2 JP3718856 B2 JP 3718856B2
- Authority
- JP
- Japan
- Prior art keywords
- magnetic field
- physical quantity
- value
- equation
- calculation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F8/00—Arrangements for software engineering
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/12—Simultaneous equations, e.g. systems of linear equations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q50/00—Systems or methods specially adapted for specific business sectors, e.g. utilities or tourism
- G06Q50/10—Services
Description
【０００１】
【産業上の利用分野】
この発明は物理量解析方法およびその装置に関し、個々の物理源の物理量ｕiと、個々の物理源に起因して生ずる任意箇所における物理量Ｏjとが数３の関係を有する物理系における未知の物理源の物理量ｕiを解析するための新規な方法およびその装置に関する。
【０００２】
【従来の技術】
従来から、物理量解析方法として、（１）複数個の磁束計の出力に基づいて、スーパーコンピュータを用い、かつモンテカルロ法に基づいて磁場源の解析を行なう方法、（２）適応ノイズキャンセラに適用された方法、（３）インパルス応答推定装置に組み込まれた方法、および（４）ニューラルネットワークを構成するニューロン素子の閾値、結合係数等を系の物理法則に基づいて設定することにより、これらの学習所要時間を不要にした方法（特開平５−９４５４３号公報参照）が提案されている。また、（５）有限要素法も提案されている。
【０００３】
上記（１）の方法は、
ａ）複数個の磁束計による探査空間に乱数を用いてｍ個の電流素片をばらまく。
ｂ）推定誤差演算プロセスにより全推定誤差（推定誤差の総和）を算出する。
ｃ）以下のｄ）からｇ）の処理を反復する。
ｄ）任意に電流素片ｋを選択し、該当する電流素片のパラメータおよび全推定誤差を退避する。
ｅ）電流素片ｋのパラメータを乱数を用いて微小な量だけ変化させる。
ｆ）推定誤差演算プロセスにより全推定誤差を算出する。
ｇ）退避した全推定誤差とｆ）で算出された全推定誤差とを比較し、退避した全推定誤差の方が小さければｄ）で退避した情報を復帰させる。
【０００４】
上記（２）の方法を組み込んだ適応ノイズキャンセラは、図２７に示すように、信号源７１からの情報にノイズ源７２からのノイズが混入した入力Ｓjを誤差演算器７３の非反転入力端子に供給し、ノイズ源７２からのノイズのみをＦＩＲフィルタ７４を介して誤差演算器７３の反転入力端子に供給し、しかも、誤差演算器７３からの出力信号をＦＩＲフィルタ７４にフィードバックしている。また、ＦＩＲフィルタ７４には、ＬＭＳ（Least Mean Square ）アルゴリズムが採用されている。
【０００５】
この適応ノイズキャンセラでは、推定ゲインを適正に設定することにより、情報に混入したノイズのみを除去でき、空調ダクト騒音の消去、自動車の車内の静音化等を達成できる。即ち、除去すべきノイズを高精度に推定できる。
上記（３）の方法を組み込んだインパルス応答推定装置は、高速フーリエ変換（以下、ＦＦＴと略称する）を採用して周波数成分の解析を行ない、解析結果に基づいてインパルス応答の推定を行なう方法である。
【０００６】
上記（４）の方法は、上記（１）の方法の不都合を解消すべく本件特許出願人が既に特許出願を行なった方法であり、図２８に示すように、入力パターンを複数の物理公式演算ユニット８１１，８１２，・・・，８１ｍに供給してそれぞれ既知の物理公式に基づく演算を行ない、全ての物理公式演算ユニット８１１，８１２，・・・，８１ｍからの出力をシグマユニット８２に供給して総和を得、総和および実際の計測値を誤差演算器８３に供給して誤差（両者の差）を得、得られた誤差を全ての物理公式演算ユニット８１１，８１２，・・・，８１ｍの補正部８１１ａ，８１２ａ，・・・，８１ｍａにフィードバックしている。また、物理公式演算ユニット８１１，８１２，・・・，８１ｍにおいて推定されている変数の値を情報収集ユニット８４により収集して解析結果として出力する。
【０００７】
したがって、物理公式自体に関しては何ら学習を行なう必要がなく、物理公式に含まれる変数のみについて誤差および推定ゲインに基づく補正を反復するだけで高精度の解析結果を得ることができる。
上記（５）の方法は、建築、機械の構造解析、気象計算、計算天文学、電磁解析等、多種多様な分野で広く利用されている方法であり、解析対象を有限個のメッシュに分割してモデルを構築し、このモデルに基づいて与えられる連立一次方程式を解くことにより物理源解析結果を得ることができる。そして、連立一次方程式を解くための方法として、ガウスの消去法、ＳＯＲ法、ＣＧ法、ＳＤ法等の種々の方法が知られている。
【０００８】
【発明が解決しようとする課題】
上記（１）の方法は、全推定誤差が小さくなるように電流素片ｋのパラメータを微小量ずつ変化させるのであるから、最終的に正しい解析結果が得られるように思われる。しかし、１回の処理を行なっても全推定誤差が小さくなるという保証が全くなく、しかも上記処理の一部についてのみ並列処理が可能であり、他の部分（ｍ個の電流素片が測定点につくる磁場の演算処理、および全推定誤差を求める処理）については並列処理が不可能であるから、スーパーコンピュータを用いても著しく長時間がかかり、しかも処理回数を増加させても最終的な解を得ることができないという不都合がある。
【０００９】
上記適応ノイズキャンセラは、ＬＭＳアルゴリズムを採用しているのであるから推定ゲインの設定が必要であり、適切な推定ゲインの設定が困難であるという不都合がある。さらに詳細に説明すると、ノイズの伝播経路は著しく多いのであり、これらの全ての伝播経路を通ってノイズキャンセル対象位置に到達する実際のノイズを高精度に推定するためには、ノイズ源７２のノイズを微小時間ずつ遅延させ、それぞれに対して適切な推定ゲインを設定しなければならなくなる。この結果、設定すべき推定ゲインの数が著しく多くなるだけでなく、全ての推定ゲインを適切に設定しなければならないことになる。もし、全ての推定ゲインを同一の値に設定するのであれば、推定処理が発散しないように推定ゲインを小さな値に設定せざるを得ず、これに伴なって、解の収束が遅くなってしまうという不都合もある。
【００１０】
上記（３）の方法は、ＦＦＴを採用しているのであるから次の不都合が生じる。ＦＦＴはサンプリング定理に基づく処理を行なう方法であるから、通常測定信号に含まれている不要な高調波を除去するためにアンチエリアシングフィルタと呼ばれるローパスフィルタを設けることが必須であり、構成が複雑化するという不都合がある。また、サンプリング区間のデータが周期的に連続するという補償がなければ周波数解析結果の精度が著しく低下してしまうのであるから、適用可能な信号の種類が制限されてしまうという不都合もある。このような不都合を解消するために、ハミング、ハニング等の窓関数を用いることが提案されているが、逆フィルタの演算時に窓関数を用いると演算後にサンプリング区間全体にわたって波形が歪み、解析精度が低下するという新たな不都合が生じる。さらに、サンプリング間隔に基づいて定まる周波数間隔の出力しか得られないのであるから、広帯域の解析が必要な場合には必然的にサンプル数を増加させなければならないという不都合もある。さらにまた、周波数軸が対数目盛であっても等間隔にサンプリングを行なわなければならないので著しく多量のメモリが必要になり、また、サンプル数も２Nでなければならないという制約があるという不都合もある。
【００１１】
上記（４）の方法は、物理公式演算結果の偏微分値と誤差との積の総和に推定ゲインを乗算して得た値で補正を行なうのであるから、推定ゲインの設定が必要である。そして、設定すべき推定ゲインが多くなるとともに、適切な推定ゲインを設定しなければならないのであるから、推定ゲインの設定が著しく繁雑になってしまうという不都合がある。簡単化のために、もし、全ての推定ゲインを同一の値に設定するのであれば、推定処理が発散しないように推定ゲインを小さな値に設定せざるを得ず、これに伴なって解の収束が遅くなってしまうという不都合もある。
【００１２】
上記（５）の方法は、メッシュ分割数を増加させることにより物理量の解析精度を高めることができる。
しかし、連立一次方程式を解くためにガウスの消去法を採用すると、未知数の数がｎである場合に、ｎ3／３のオーダーの演算回数が必要であり、著しく演算負荷が大きくなるので、未知数の数が多い場合には、スーパーコンピュータを用いても実用上十分な速度での物理源解析を達成することができない。これに対してコレスキー法を採用すれば、演算回数をｎ3／６のオーダーにすることができるが、サイズが大きい連立一次方程式を解くための方法としては演算負荷が大きすぎる。そして、何れの方法においても誤差が累積されることになるので、物理源解析精度を余り高めることができない。
【００１３】
上記ＳＯＲ法は、ガウス＝ザイデル法の収束を加速するために、加速パラメータを導入して過大修正を行なう方法であるから、所定の条件を満足する係数マトリクスに関しては解の収束が保証されているが、所定の条件を満足しない係数マトリクスに関しては解の収束が保証されないという不都合がある。そして、解の収束が保証される係数マトリクスであっても、加速パラメータの最適値を簡単にかつ確実に求める方法がなく、加速パラメータの設定が不適当な場合には、収束を十分に加速することができない。
【００１４】
上記ＳＤ法は、最急降下法と呼ばれる方法であり、ＳＯＲ法のような係数マトリクスの制約はないが、反復回数が著しく多く、しかも最急降下方向の算出、未知数を補正するためのゲインの算出を毎回行なう必要があるとともに、これらの演算はベクトル演算および除算を含むのであるから、全体として演算負荷が著しく大きくなってしまう。
【００１５】
上記ＣＧ法は、共役傾斜法と呼ばれる方法であり、既に修正した方向の全てと直交する方向を修正方向として順次未知数を補正するのであるから、理論的には反復回数を著しく少なくすることができ（ｎ回以下にすることができ）、しかも物理源解析精度を著しく高くすることができる。しかし、実際には、殆ど全ての補正処理において誤差が発生するのであるから、係数マトリクスによって反復回数が大きくばらついてしまう。また、各回の処理を行なうに当ってベクトル演算、マトリクス演算が必須になるのであるから、全体として演算負荷が著しく大きくなってしまう。
【００１６】
【発明の目的】
この発明は上記の問題点に鑑みてなされたものであり、推定ゲインの設定が不要であり、しかも係数マトリクス、サンプル数の制約を排除して高速かつ高精度に物理量の解析を行なうことができる物理量解析方法およびその装置を提供することを目的としている。
【００１７】
【課題を解決するための手段】
上記の目的を達成するための、請求項１の物理量解析装置は、測定対象物が発生する磁場を複数の磁場センサで計測し、前記測定対象物内部の磁場源における物理量を解析する装置であって、
（１） 個々の磁場源ｉに対応する単位量の磁場源が各所定箇所に発生する磁場を示す、数１で与えられる定数αijを予め定めておく手段と、
（２） 物の内部に設定された格子点ｋ（ｋ＝１，２・・・，ｐ）に位置する各磁場源の、数２で示すx方向電流成分およびｙ方向電流成分からなる物理量ｕiを仮に定めておく手段と、
（３） 数３の演算を行なって観測点ｊに発生すると予測される磁場Ｏjを算出する手段と、
（４） 観測点ｊにおける磁場の実測値Ｓjを測定する手段と、
（５） 観測点ｊにおける磁場の実測値Ｓjと上記磁場Ｏjとの差Ｓj−Ｏjを算出する手段と、
（６） 数４の演算を行なって仮の解としての物理量ｕiを得る手段と、
（７） （３）から（６）までの各手段による処理を全ての磁場源ｉについて行なわせる手段と、
（８） 上記差Ｓj−Ｏjの二乗の総和が所定の閾値よりも小さくなるまで、（３）から（６）までの各手段による処理の反復処理を反復させる手段と、
（９） 最終的に得られた各磁場源の物理量ｕiを磁場源の物理量ｕiとして採用する手段とを含むものである。
【００４６】
請求項１の物理量解析装置であれば、測定対象物が発生する磁場を複数の磁場センサで計測し、前記測定対象物内部の磁場源における物理量を解析する装置であって、
（１） 個々の磁場源ｉに対応する単位量の磁場源が各所定箇所に発生する磁場を示す、数１で与えられる定数αijを予め定めておく手段と、
（２） 物の内部に設定された格子点ｋ（ｋ＝１，２・・・，ｐ）に位置する各磁場源の、数２で示すx方向電流成分およびｙ方向電流成分からなる物理量ｕiを仮に定めておく手段と、
（３） 数３の演算を行なって観測点ｊに発生すると予測される磁場Ｏjを算出する手段と、
（４） 観測点ｊにおける磁場の実測値Ｓjを測定する手段と、
（５） 観測点ｊにおける磁場の実測値Ｓjと上記磁場Ｏjとの差Ｓj−Ｏjを算出する手段と、
（６） 数４の演算を行なって仮の解としての物理量ｕiを得る手段と、
（７） （３）から（６）までの各手段による処理を全ての磁場源ｉについて行なわせる手段と、
（８） 上記差Ｓj−Ｏjの二乗の総和が所定の閾値よりも小さくなるまで、（３）から（６）までの各手段による処理の反復処理を反復させる手段と、
（９） 最終的に得られた各磁場源の物理量ｕiを磁場源の物理量ｕiとして採用する手段とを含むものであるから、推定ゲインの設定が不要になり、推定ゲインがないことに伴って解の収束を高速化できる。即ち、磁場源の解析を高速に達成することができる。さらに、磁場のｚ方向成分のみを計測することにより、各格子点におけるx方向電流成分、ｙ方向電流成分を解析することができる。
【００７３】
【実施例】
以下、実施例を示す添付図面によって詳細に説明する。
図１はこの発明の物理量解析方法の一実施例を説明するフローチャートであり、複数（ｍ個）の微小物理源の位置および複数（ｎ個）の観測点の位置がそれぞれ設定され、かつ系に線形加算性が成立する場合に、ステップＳＰ１において各微小物理源ｉ（ｉ＝１，２，・・・，ｍ）と各観測点ｊ（ｊ＝１，２，・・・，ｎ）の位置的関係から物理法則を用いて定数αｉｊを算出し、ステップＳＰ２において各微小物理源ｉの物理量ｕｉを仮に設定する（例えば、全ての物理量を０に設定する）。そして、ステップＳＰ３において数１の演算を行なうことにより観測点ｊに発生すると予測される物理場Ｏｊを算出し、ステップＳＰ４において観測点ｊにおける物理場の実測値Ｓｊを得る。
【００７４】
次いで、ステップＳＰ５において、実測値Ｓｊと物理場Ｏｊとの差および定数αｉ jを用いて数４の演算を行なうことにより仮の解としての物理量ｕｉを得る。
尚、ステップＳＰ３，ＳＰ４，ＳＰ５の処理は全ての微小物理源ｉについて順次行なわれる。
その後、ステップＳＰ６において実測値Ｓｊと物理場Ｏｊとの差の二乗の総和が十分に小さくなったか否かを判別し、十分には小さくなっていないと判別された場合には、再びステップＳＰ３の処理を行なう。逆に、十分に小さくなったと判別された場合には、ステップＳＰ７において、最終的に得られている物理量ｕｉを解析結果として取り出し、そのまま一連の処理を終了する。
【００７５】
さらに詳細に説明すると、誤差関数Ｅとして数５に示すように、実測値Ｓｊと物理場Ｏｊとの差の二乗の総和を採用する。
【００７６】
【数５】
【００７７】
ここで、
Ｓｊ−Ｏｊ＝（Ｓｊ−Ｏｊ＋αｉｊｕｉ）−αｉｊｕｉ
であるから、未知数の物理量ｕｉを含まない項と未知の物理量ｕｉを含む項とに分離できる。したがって、数５は数６と等価である。
【００７８】
【数６】
【００７９】
但し、Ａｉ ２，ＡｉＢｉ，Ｃiはそれぞれ数７，数８，数９で与えられる。
【００８０】
【数７】
【００８１】
【数８】
【００８２】
【数９】
【００８３】
数６から明らかなように、誤差関数Ｅは，図２に示すように、ｕｉに関して軸をＢｉ/Ａｉ、最小値をＣｉ−Ｂｉ ２とする放物線になる。即ち、誤差関数Ｅは，全体としてｍ次元放物面体であり、一意収束性がある。したがって、ｕｉの軸値への修正（ｕｉ←Ｂｉ/Ａｉ）は解への漸近を意味することになる。また、上記修正式に数７、数８を代入して整理することにより数４が得られる。
【００８４】
図３は分布電流を格子近似した磁場源モデルと磁場センサ配列面との関係を示す概略図であり、磁場センサにより対象物から発生する磁場を検出し、検出された磁場に基づいて対象物の内部の電流分布を測定する場合における磁場源モデルと磁場センサ配列面との関係を示している。図３においては、ｎ個の磁場センサが配列されているとともに、ｍ個の磁場源が想定されている。そして、磁場のｚ方向成分のみの計測を行なうと仮定すれば、磁場源はｘ方向とｙ方向の電流成分を有することになる。また、ｐ個の３次元格子で関心領域を近似すれば、各格子座標にｘ方向とｙ方向の２個の独立した電流双極子を割り当てることになる。即ち、ｍ＝２ｐとなる。
【００８５】
格子点ｋ（ｋ＝１，２，・・・，Ｐ）のｘ方向電流成分をＰｘｋ、ｙ方向電流成分をＰｙｋとし、その座標を（Ｘｋ，Ｙｋ，Ｚｋ）とし、観測点ｊの座標を（ｘｊ，ｙｊ，ｚｊ）とする。そして、磁場センサのピックアップコイル径を無視し、コイル中心の磁束密度を計測していると仮定すれば、数２の下で、ビオサバールの法則により、αｉｊは数１となる。そして、このように設定することにより、数３と同様に線形加算が成立する系になるので、前記と同様の処理を行なうことができ、磁場源の内部の電流分布を算出することができる。
【００８６】
そして、数４による推定を安定化するために、事前に正規化（αｉｊ’＝αｉｊ/Ａｉ）を行なって各放物線の開口度を１に統一する。また、正規化後の定数αｉｊ’による推定値を補正するためにｕｉ＝ｕｉ’/Ａｉの演算を行なう。
以上の一連の解析処理を行なうことにより、各格子点における電流成分を高速かつ正確に推定できた。
【００８７】
この実施例は上記具体例に限定されるものではなく、例えば、集光特性が既知の光量センサを採用することにより光量解析方法に適用することが可能であり、その他、指向性を持ち、かつ指向性が既知のセンサを採用することにより種々の解析方法に適用することが可能である。
【００８８】
【実施例２】
図４はこの発明の物理量解析装置の一実施例を示すブロック図であり、解析対象となる物理源の数と等しい数の格子ユニット１１，１２，・・・，１ｍと、各格子ユニット１１，１２，・・・，１ｍに対応する修正ユニット１１ａ，１２ａ，・・・，１ｍａとからなる格子ユニット層１と、観測点の数と等しい数の観測ユニット２１，２２，・・・，２ｎと、各観測ユニット２１，２２，・・・，２ｎに対応する物理場演算ユニット２１ａ，２２ａ，・・・，２ｎａとからなる観測ユニット層２と、格子ユニット層１、観測ユニット層２を制御する制御回路３とを有している。ここで、制御回路３は、全ての観測ユニット２１，２２，・・・，２ｎに対して互に同じタイミングで格納指示信号を供給し、各格子ユニット１１，１２，・・・，１ｍに対して順次修正指示信号を供給するものである。
【００８９】
図５は格子ユニット１ｉおよび修正ユニット１ｉａの構成を詳細に示すブロック図であり、修正量基本値演算セル３１ｉと、修正値演算セル３２ｉと、推定値保持セル３３ｉとを有している。
上記修正量基本値演算セル３１ｉは、ｎ個の入力端子を有しているとともに、各入力端子にそれぞれ定数乗算機能を持たせており、各観測ユニット２１，２２、・・・、２ｎからの出力値（Ｓｊ−Ｏｊ）を対応する入力端子に供給している。そして、定数αｉｊおよび出力値（Ｓｊ−Ｏｊ）に基づいて数１０で示す修正量基本値を出力する。
【００９０】
【数１０】
【００９１】
上記修正値演算セル３２ｉは、２個の入力端子を有しているとともに、各入力端子にそれぞれ１、数１１の定数を乗算する定数乗算機能を持たせており、１の定数を乗算する入力端子に推定値保持セル３３ｉから出力される今回の推定前の推定値ｕｉを供給しているとともに、他方の入力端子に修正量基本値を供給している。そして、定数および供給された値に基づいて数４で示す修正値を得、推定値保持セル３３ｉに供給する。
【００９２】
【数１１】
【００９３】
上記推定値保持セル３３ｉは、制御回路３から修正指示信号ｉが供給されたことに応答して、修正値演算セル３２ｉから供給される修正値を新たに保持するものである。
尚、例えば、上記修正量基本値演算セル３１ｉとしては、ニューロンデバイスを採用することが可能であり、しかもニューロンデバイスの飽和関数のうちリニアな部分のみを使用すれば足りるのであるから、ニューロンデバイスに代えて、乗算器、加算器のみを用いて修正量基本値演算セル３１ｉを構成することが可能である。
【００９４】
図６は観測ユニット２ｉおよび物理場演算ユニット２ｉａの構成を詳細に示すブロック図であり、物理場演算セル４１ｉと、誤差演算セル４２ｉと、観測値保持セル４３ｉとを有している。
上記物理場演算セル４１ｉは、ｍ個の入力端子を有しているとともに、各入力端子にそれぞれ定数乗算機能を持たせており、各格子ユニット１１，１２，・・・１ｍからの出力値ｕｉを対応する入力端子にを供給している。そして、定数αｉｊおよび出力値ｕｉに基づいて数３で示す物理場演算値を出力する。
【００９５】
上記観測値保持セル４３ｉは、制御回路３から格納指示信号ｊが供給されたことに応答して、観測値を新たに保持するものである。
上記誤差演算セル４２ｉは、所定の物理量を検出するセンサが検出した観測値から物理場演算値を減算して得た値を誤差として出力するものである。
したがって、制御回路３により修正指示信号ｉおよび格納指示信号ｊを出力することにより図１のフローチャートと同様に高速かつ高精度に物理量の解析を行なうことができる。
【００９６】
図７は、格子ユニット１ｉおよび修正ユニット１ｉａの構成を詳細に示すブロック図であり、図５の構成例と異なる点は、修正値演算セル３２ｉを省略し、修正量基本値演算セル３１ｉに、推定値を入力とし、かつ１の定数に基づく乗算機能を有する入力端子を追加するとともに、他の各入力端子に、図５の各入力端子が定数αｉｊに対して数１１を乗算して得た定数を乗算機能として持たせて新たな修正値演算セルを構成した点のみである。
【００９７】
したがって、この場合には、構成を簡素化でき、しかも図５と同様の作用を達成できる。
図８は観測ユニット２ｉおよび物理場演算ユニット２ｉａの構成を詳細に示すブロック図であり、図６の構成例と異なる点は、誤差演算セル４２ｉを省略し、物理場演算セル４１ｉの各入力端子に、図６の各入力端子が定数αijに対して−１を乗算して得た定数を乗算機能として持たせるとともに、観測値を入力とし、かつ１の定数に基づく乗算機能を有する入力端子を追加して物理場誤差演算セルを構成した点のみである。
【００９８】
したがって、この場合にも、構成を簡素化でき、しかも図６と同様の作用を達成できる。
図９は格子ユニット１ｉおよび修正ユニット１ｉａの構成を詳細に示すブロック図であり、図７の構成例と異なる点は、順次供給される誤差を保持する誤差格納セル（例えば、サンプルホールド回路）３４ｉ１，３４ｉ２，・・・，３４ｉｎをさらに有している点および推定値ｕiに対して定数αijを乗算する乗算セル３５ｉをさらに有している点のみである。
【００９９】
したがって、この場合には、観測ユニットを１つのみにすることにより構成を簡素化でき、しかも図７と同様の作用を達成できる。
【０１００】
【実施例３】
図１０はこの発明の他の実施例としての適応ノイズキャンセル方法を説明するフローチャートであり、ステップＳＰ１において、現時点を基準とする過去のノイズのリファレンスαj-kを得、ステップＳＰ２において各ノイズの実測値に対する相関係数（混ざり具合を示す値）ｕｉを仮に設定する。そしてステップＳＰ３において数１２の演算を行なうことにより観測点ｊに発生すると予測されるノイズＯｊを算出し、ステップＳＰ４において観測点ｊにおける実測値Ｓｊを得る。
【０１０１】
【数１２】
【０１０２】
次いで、ステップＳＰ５において、実測値Ｓｊと予測ノイズＯｊとの差およびリファレンスαｋ−ｊを用いて数１３の演算を行なうことにより仮の解としての相関係数ｕｉを得る。
【０１０３】
【数１３】
【０１０４】
尚、ステップＳＰ３，ＳＰ４，ＳＰ５の処理は全ての時点のノイズｉについて順次行なわれる。
その後、ステップＳＰ６において実測値Ｓjと予測ノイズＯjとの差の二乗の総和が十分に小さくなったか否かを判別し、十分には小さくなっていないと判別された場合には、再びステップＳＰ３の処理を行なう。逆に、十分に小さくなったと判別された場合には、ステップＳＰ７において、最終的に得られている相関係数ｕiを解析結果として採用し、ステップＳＰ８において、情報とノイズが混在した信号からリファレンスαijおよび相関係数ｕiに基づいて得られる予測ノイズＯj を減算して情報のみを抽出し、そのまま一連の処理を終了する。
【０１０５】
図１１はこの適応ノイズキャンセル方法の具体的適用例を示す波形図であり、５０Ｈｚの交流成分が重畳された心電図｛図１１（Ａ）参照｝から交流成分を除去した場合｛図１１（Ｂ）参照｝を示している。
図１１から明らかなように、従来の適応ノイズキャンセラのように推定ゲイン（ステップサイズとも呼ばれる）を設定する必要がなく、しかも高精度にノイズを除去して正確な心電図が得られる。
【０１０６】
【実施例４】
図１２はこの発明のさらに他の実施例としての音響探査方法を説明するフローチャートであり、ステップＳＰ１において、各サンプリング時点における送信音波αｋを得、ステップＳＰ２において探査対象となる物理源のインパルス応答ｕｋを仮に設定する。そして、ステップＳＰ３において数１４の演算を行なうことにより観測点ｊに発生すると予測される予測受信音波Ｏｊを算出し、ステップＳＰ４において観測点ｊにおける実測値Ｓｊを得る。但し、ｍは送信音波の１周期におけるサンプル数、ｎは受信音波の総サンプル数である。
【０１０７】
【数１４】
【０１０８】
次いで、ステップＳＰ５において、実測値Ｓｊと予測受信音波Ｏｊとの差および送信音波αｋを用いて数１５の演算を行なうことにより仮の解としての相関係数ｕｋを得る。
【０１０９】
【数１５】
【０１１０】
尚、ステップＳＰ３，ＳＰ４，ＳＰ５の処理は該当する全ての時点の送信音波αkについて順次行なわれる。
その後、ステップＳＰ６において実測値Ｓjと予測ノイズＯjとの差の二乗の総和が十分に小さくなったか否かを判別し、十分には小さくなっていないと判別された場合には、再びステップＳＰ３の処理を行なう。逆に、十分に小さくなったと判別された場合には、ステップＳＰ７において、最終的に得られているインパルス応答ｕk解析結果として採用し、そのまま一連の処理を終了する。
【０１１１】
尚、数１４、数１５の上段においてｍ個の総和を算出しているのは、送信音波が１周期分のみであると仮定しているからであり、数１５の下段において（ｎ-ｋ）個の総和を算出しているのは、受信音波の観測時間が有限であり、送信音波の一部のみに対応する受信音波しか得られていないからである。
図１３は音響探査の具体例を示す波形図であり、図１３（Ａ）に示す送信音波に基づいて図１３（Ｂ）に示す受信音波が得られた場合に、ノイズが存在していなければ図１３（Ｃ）に示すインパルス応答推定結果を得ることができ、ノイズ（１％）が存在している場合に図１３（Ｄ）に示すインパルス応答推定結果を得ることができた。
【０１１２】
図１３から明らかなように、ノイズが存在していなければ著しく高精度にインパルス応答を推定でき、ある程度のノイズが存在していてもかなり高い精度でインパルス応答を推定できることが分る。
尚、この実施例は音波のみについて説明したが、電磁波等にも適用できることはもちろんである。
【０１１３】
実施例３、実施例４の各方法は、実施例１の定数に代えて、解析目的に適合する値を採用し、推定値も実施例１の推定値と異なるのであるが、採用された値、推定値に基づく処理は実施例１とほぼ同様であるから、実施例２の装置に適用することにより簡単に実施例３、実施例４に対応する装置を得ることができる。但し、適応ノイズキャンセラに適用する場合には、各観測ユニット２１，２２，・・・，２ｎにそれぞれ別個の観測値を供給する代わりに、観測値が１つの観測ユニットに供給され、各観測ユニットの値が順次隣合う観測ユニットに供給されるよう構成すればよい。
【０１１４】
【実施例５】
図１４はこの発明の物理量解析方法のさらに他の実施例を説明するフローチャートであり、個々の物理源を含む、線形連立一次方程式が成立する領域内において、個々の物理源の物理量ｕiおよび上記領域に基づいて定まる比例定数αijに基づいて任意箇所で得られる物理量Ｏjを解析する方法を説明している。
【０１１５】
ステップＳＰ１において、領域の性質等に基づいて比例定数αijを得、ステップＳＰ２において各物理源の物理量ｕｉを仮に設定する（例えば、全ての物理量を０に設定する）。そして、ステップＳＰ３において数３の演算を行なうことにより任意箇所で得られる物理量Ｏｊを算出し、ステップＳＰ４において成立条件（領域内に与えた既知の物理量）Ｓｊと物理量Ｏｊとの差を算出する。
【０１１６】
次いで、ステップＳＰ５において上記差および比例定数αijを用いて数４の演算を行なうことにより仮の解としての物理量ｕｉを得る。
尚、ステップＳＰ３，ＳＰ４，ＳＰ５の処理は全ての物理源について順次行なわれる。
その後、ステップＳＰ６において成立条件Ｓｊと物理量Ｏｊとの差の二乗の総和が十分に小さくなったか否かを判別し、十分には小さくなっていないと判別された場合には、再びステップＳＰ３の処理を行なう。逆に、十分に小さくなったと判別された場合には、ステップＳＰ７において、最終的に得られている物理量ｕｉを解析結果として取り出し、そのまま一連の処理を終了する。
【０１１７】
さらに詳細に説明すると、誤差関数Ｅとして数５に示すように成立条件Ｓｊと物理量Ｏｊとの差の二乗の総和を採用すれば、数６から数９の関係が得られる。したがって、ｕｉの軸値への修正により解への漸近を達成できるとともに、収束を高速化できる。また、ｕｉの軸値への修正に当って、数７に基づく除算が必要であるが、数７の値は事前に算出可能な定数であるから、予め数７の逆数を算出しておくことにより、上記一連の処理における除算回数を著しく低減することができ、ひいては解析所要時間を一層短縮することができる。
【０１１８】
しかし、一連の処理を行なう毎に数３の演算を行なう必要があり、演算負荷を十分には低減できていない。この点に関連して、本件発明者は、物理量ｕｉを軸値に修正するための修正量△ｕｉ、対応する比例定数αij、および先行する物理量Ｏｊに基づいて数１６の演算を行なうことにより、物理量Ｏｊを修正することができることを見出した。
【０１１９】
【数１６】
【０１２０】
したがって、一連の処理を行なう毎に数３の演算を行なう代わりに数１６の演算を行なうことにより、演算負荷を大幅に低減し、ひいては解析所要時間をさらに短縮することができる。
さらに、予めＤｊ＝Ｓｊ-Ｏｊを算出しておき、数１６の演算を行なう代わりに数１７の演算を行なうようにすれば、Ｓｊ、Ｏｊをそれぞれ保持しておく必要がなくなり、メモリ容量を低減することができる。
【０１２１】
【数１７】
【０１２２】
この実施例を有限要素法に適用すれば、一連の処理を行なう毎に必要な演算負荷についても上述のように著しく低減することができる。したがって、係数マトリクスの次数が著しく大きい解析対象領域についても簡単に、かつ高速に物理量の解析を行なうことができる。
図１５は多数の電極により収集した体表面電位分布データに基づいて体内電気活動状態を可視化するための装置を概略的に示す図であり、図１６は図１５の装置を実現する処理の一部として電位パターンを得るための物理量解析方法の具体例を説明するフローチャートである。
【０１２３】
図１５の装置は、人体の頭部、胸部等に被着される基部部材２０１の所定位置に所定の配置パターンで電位計２０２が装着されてあり、電位計２０２からの出力信号が生体計測用絶縁アンプ２０３を介してＡ／Ｄコンバータ２０４に導かれてディジタル信号に変換され、ディジタル信号が処理部２０５に供給される。この処理部２０５は電位パターン保持部２０６を有しており、処理部２０５において図１６のフローチャートの処理を行なって得た電位パターンを電位パターン保持部２０６に保持させる。また、電位パターンが得られた後に、電位パターンに基づいて人体内部の電流の解析を行なう電流解析部２０７と、電流解析部２０７による解析結果を表示するための表示部２０８とを有している。
【０１２４】
上記電位計２０２は、例えば人体の表面に圧接できるように付勢機構を有しており、しかも、人体に対する圧接状態において、基部部材２０１に対する相対位置を検出できるように構成してあることが好ましい。
上記電流解析部２０７は、得られた電位パターンを比例定数αijとして採用し、例えば、実施例１と同様の処理を行なって人体内部の電流の解析を行なうものである。
【０１２５】
また、図１６のフローチャートに基づく処理を行なう前に、予めＭＲＩ、Ｘ線ＣＴ等により体内臓器の３次元形状を抽出し、別な計測法により得られ、または医学上の既知情報に基づいて得られた電気抵抗率を用いて、解析対象となる３次元領域を所定ピッチ（例えば５ｍｍピッチ）の３次元格子状抵抗ネットワークでモデル化しておく。また、上記装置は体表面電位分布データに基づいて体内電気活動状態を可視化するものであるから、その前提として、上記３次元領域内の任意の抵抗の一方の端部に単位電流を与えた状態における表面電位分布を、上記３次元領域内に含まれる全ての抵抗について算出しておくことが必要になる。ここで、３次元格子状抵抗ネットモデルのx、y、z方向のノード数をそれぞれＮｘ、Ｎｙ、Ｎｚで表し、図１７に示すように、ノード[ｋx，ｋy，ｋz]の電位をｕ[ｋx，ｋy，ｋz]で表し、ノード[ｋx，ｋy，ｋz]のｘ方向に接続されている抵抗の値をＲx[ｋx−１，ｋy，ｋz]、Ｒx[ｋx，ｋy，ｋz]で表し、ノード[ｋx，ｋy，ｋz]のｙ方向に接続されている抵抗の値をＲy[ｋx，ｋy-１，ｋz]、Ｒy[ｋx，ｋy，ｋz]で表し、ノード[ｋx，ｋy，ｋz]のｚ方向に接続されている抵抗の値をＲz[ｋx，ｋy，ｋz−１]、Ｒz[ｋx，ｋy，ｋz]で表し、外部からノード[ｋx，ｋy，ｋz]に強制的に供給される電流をＳ[ｋx，ｋy，ｋz]で表し、ノード[ｋx，ｋy，ｋz]から接続されている抵抗を通じて流出する電流をＯ[ｋx，ｋy，ｋz]で表す。そして、各抵抗値Ｒx[ｋx−１，ｋy，ｋz]、Ｒx[ｋx，ｋy，ｋz]、Ｒy[ｋx，ｋy−１，ｋz]、Ｒy[ｋx，ｋy，ｋz]、Ｒz[ｋx，ｋy，ｋz−１]、Ｒz[ｋx，ｋy，ｋz]の逆数をそれぞれＧx[ｋx−１，ｋy，ｋz]、Ｇx[ｋx，ｋy，ｋz]、Ｇy[ｋx，ｋy−１，ｋz]、Ｇy[ｋx，ｋy，ｋz]、Ｇz[ｋx，ｋy，ｋz−１]、Ｇz[ｋx，ｋy，ｋz]で表し、実際に接続されている抵抗に対応する総和ＧN[ｋx，ｋy，ｋz]を数１８で表す。
【数１８】
【０１２６】
但し、数１８の右辺においては実際に接続されている抵抗に対応する項のみを加算の対象とする。
したがって、オームの法則に基づいて、流出電流Ｏ[ｋx，ｋy，ｋz]は数１９により算出することができる。
【０１２７】
【数１９】
【０１２８】
また、キルヒホッフの第１法則を成立させるためには、全てのノードにおいて供給電流Ｓ[ｋx，ｋy，ｋz]と流出電流Ｏ[ｋx，ｋy，ｋz]とが等しくなるノード電位ｕ[ｋx，ｋy，ｋz]を求めればよい。尚、キルヒホッフの第２法則に関しては、未知の物理量をノード電位ｕ[ｋx，ｋy，ｋz]とすることにより自動的に満足される。即ち、数２０で示される評価関数Ｅの値を最小にするノード電位ｕ[ｋx，ｋy，ｋz]を算出すればよいことになる。
【０１２９】
【数２０】
【０１３０】
上記の条件下において、ステップＳＰ１において（Ｎx−１）×Ｎy×Ｎz個のＧx[ｋx，ｋy，ｋz]、Ｎx×（Ｎy−１）×Ｎz個のＧy[ｋx，ｋy，ｋz]、Ｎx×Ｎy×（Ｎz−１）個のＧｚ[ｋx，ｋy，ｋz]、Ｎx×Ｎy×Ｎz個のＧN[ｋx，ｋy，ｋz]、およびＮx×Ｎy×Ｎz個の、数２１の関係を有するβ[ｋx，ｋy，ｋz]をそれぞれ算出して格納する。但し、数２１の右辺の各項は実際の接続を反映するように求められている。したがって、解析対象となるノードの数が著しく多くなっても、マトリクスの列方向の非０の比例定数の数は最大でも７にしかならない。したがって、係数マトリクスを保持しておくためのメモリ容量を、マトリクスの次数に比して十分に少なくすることができる。
【数２１】
【０１３１】
そして、ステップＳＰ２においてＮx×Ｎy×Ｎz個の推定誤差Ｄ[ｋx，ｋy，ｋz]＝Ｓ[ｋx，ｋy，ｋz]−Ｏ[ｋx，ｋy，ｋz]を算出して格納する。
次いで、ステップＳＰ３において数２２に基づいて未知のノード電位ｕ[ｋx，ｋy，ｋz]を軸値に修正するための修正量Δｕを算出してノード電位ｕ[ｋx，ｋy，ｋz]を修正する。
【０１３２】
【数２２】
【０１３３】
そして、ステップＳＰ４において数２３の演算を行なって推定誤差を修正する。
その後、ステップＳＰ５においてｋzがＮzに達したか否かを判別し、達していなければステップＳＰ６においてｋzをインクリメントし、再び、ステップＳＰ３の処理を行なう。逆に、ステップＳＰ５においてｋz＝Ｎzであると判別された場合には、ステップＳＰ７においてｋyがＮyに達したか否かを判別し、達していなければステップＳＰ８においてｋyをインクリメントし、再びステップＳＰ３の処理を行なう。逆に、ステップＳＰ７においてｋy＝Ｎyであると判別された場合には、ステップＳＰ９においてｋxがＮxに達したか否かを判別し、達していなければステップＳＰ１０においてｋxをインクリメントし、再びステップＳＰ３の処理を行なう。逆に、ステップＳＰ９においてｋx＝Ｎxであると判別された場合には、ステップＳＰ１１において評価関数Ｅの値が所定の閾値よりも小さいか否かを判別し、評価関数Ｅの値が所定の閾値よりも小さければステップＳＰ１２において、最終的に得られたノード電位ｕ[ｋx，ｋy，ｋz]を解析結果として採用し、そのまま一連の処理を終了する。逆に、評価関数Ｅの値が所定の閾値以上であれば、ｋx，ｋy，ｋzを初期値にリセットして再びステップＳＰ３の処理を行なう。
【数２３】
【０１３４】
以上の一連の処理を行なうことにより、人体の頭部、胴体を解析対象領域とする未知の電位の解析を達成することができた。具体的には、人体の頭部を解析対象とする場合には未知の電位の数が約４００００であり、胴体を解析対象とする場合には未知の電位の数が約４８００００であるが、主記憶が６４Ｍバイトのワークステーションを用いて十分な精度の解析結果を得ることができた。また、所要時間に関しても、前者の場合に３０００回の反復処理を行なって１４４秒、後者の場合に３０００回の反復処理を行なって１９８６秒であった。
【０１３５】
以上のようにして得られた電位パターンを比例定数αijとして採用し、体表面電位分布データを既知の物理量Ｓjとして採用し、実施例１と同様の処理を行なって生体内部に設定された格子点の電流の流入流出状況を解析することができる。
さらに詳細に説明する。
【０１３６】
格子点［１，１，１］に単位電流が流入し、格子点［２，１，１］から電流が流出すると仮定して得た体表面電位分布データを、（格子点数−１）×３列の比例定数αijのうち、第１列の比例定数として採用し、格子点［１，１，１］に単位電流が流入し、格子点［１，２，１］から電流が流出すると仮定して得た体表面電位分布データを第２列の比例定数として採用し、第１列の比例定数として採用し、格子点［１，１，１］に単位電流が流入し、格子点［１，１，２］から電流が流出すると仮定して得た体表面電位分布データを第３列の比例定数として採用し、以下、同様にして他の全ての列の比例定数を採用する。
【０１３７】
このようにして全ての列の比例定数αijを採用すれば、数３の関係を満足することになるので、実施例１と同様の解析処理を行なうことにより、生体内部に設定された格子点の電流の流入流出状況を解析することができる。
具体的には、図１８中の（Ａ）の２次元等電位線図に示される電流源を仮定して図１６の処理を行なって電位分布データを得、この電位分布データを比例定数として採用して実施例１と同様の解析処理を行なった結果、図１８の（Ｂ）に示す電流源解析結果を得ることができた。
【０１３８】
尚、生体内部の電流分布を解析するためには、生体内部に設定された格子点の電流の流入流出状況の解析結果に基づいて図１６のフローチャートと同様の解析処理を行なう必要がある。但し、当初に既知の物理量として生体内部の供給電流Ｓ［ｋx，ｋy，ｋz］および流出電流Ｏ［ｋx，ｋy，ｋz］を採用する代わりに、抵抗を流れる電流を採用しておけば、図１６のフローチャートと同様の解析処理を行なって得た電位分布データを比例定数として採用して実施例１と同様の解析処理を行なうことにより、生体内部の電流分布を得ることができる。
【０１３９】
【実施例６】
図１９はこの発明の物理量解析方法のさらに他の実施例を説明するフローチャートであり、ステップＳＰ１において、解析対象となる多数の物理量の中から任意の２つの物理量を選択し、ステップＳＰ２において、選択した２つの物理量の一方（以下、物理量ｕkと称する）を対象として、実施例５と同様の方法により物理量ｕkを１回だけ修正し、ステップＳＰ３において他方の物理量（以下、物理量ｕiと称する）を対象として、実施例５と同様の方法により物理量ｕiを１回だけ修正する。次いで、ステップＳＰ４，ＳＰ５において、ステップＳＰ２，ＳＰ３と同様にして再び物理量ｕk，ｕiを修正する。その後、ステップＳＰ６において、ステップＳＰ３，ＳＰ４，ＳＰ５の処理により得られた修正量に基づいて、両物理量ｕk，ｕiに基づいて定まる評価関数Ｅの値を最小にするために必要な両物理量ｕk，ｕiの修正量Δｕk，Δｕiを算出し、ステップＳＰ７において、算出された修正量Δｕk，Δｕiに基づいて各物理量ｕk，ｕiを修正する。そして、ステップＳＰ８において推定誤差Ｄを上記修正量および対応する比例定数に基づいて修正し、ステップＳＰ９において全ての物理量に対する処理が行なわれたか否かを判別し、処理が行なわれていない物理量が存在すれば、ステップＳＰ１０において、他の２つの物理量を選択し、再びステップＳＰ２の処理を行なう。逆に、ステップＳＰ９において全ての物理量に対する処理が行なわれたと判別された場合には、ステップＳＰ１１において評価関数Ｅの値が所定の閾値よりも小さいか否かを判別し、評価関数Ｅの値が所定の閾値よりも小さければステップＳＰ１２において、最終的に得られた物理量を解析結果として採用し、そのまま一連の処理を終了する。逆に、評価関数Ｅの値が所定の閾値以上であれば、再びステップＳＰ１の処理を行なう。
【０１４０】
さらに詳細に説明する。
実施例５においては、修正すべき物理量が１つであるから、評価関数Ｅは放物線である。しかし、この実施例においては、修正すべき物理量が２つであるから、評価関数Ｅは図２０に示すように楕円になる。そして、ステップＳＰ２の処理を行なうことにより、図２０中にｐ1で示すように楕円状の点が得られ、次いで、ステップＳＰ３，ＳＰ４，ＳＰ５の処理を順次行なうことにより、楕円の中心により近い楕円状の点、ｐ2，ｐ3，ｐ4が得られる。これらの点ｐ1，ｐ2、ｐ3，ｐ4の関係は、図２１に示すとおりであり、ｐ1とｐ2との間隔をΔｕi´´、ｐ3とｐ4との間隔をΔｕi´、ｐ2とｐ3との間隔をΔｕk´、ｐ1とｐ3との間隔をｒ1、ｐ3と楕円中心との間隔をｒ2、ｐ2と楕円中心とのｕｉ軸方向の間隔をΔｕi、ｐ4と楕円中心とのｕk軸方向の間隔をΔｕkとすれば、数２４から数２６の比例関係式が得られる。
【０１４１】
【数２４】
【０１４２】
【数２５】
【０１４３】
【数２６】
【０１４４】
したがって、数２４から数２７が得られる。
【０１４５】
【数２７】
【０１４６】
そして、数２７を数２５，数２６に代入することにより、数２８のとおり、Δｕk，Δｕiを得ることができる。
【０１４７】
【数２８】
【０１４８】
また、上記ステップＳＰ８における推定誤差Ｄの修正は、Ｄj←Ｄj−αijΔｕk−αkjΔｕiの演算により達成される。
以上から明らかなように、任意の２つの物理量を選択し、各物理量について交互に２回ずつ修正処理を行なった後は、両物理量を楕円中心に遷移させるために必要な修正量を簡単に算出することができる。
【０１４９】
【実施例７】
図２２はこの発明の物理量解析装置の他の実施例を示すブロック図であり、解析対象となる物理源の数と等しい数の格子ユニット１１１，１１２，・・・，１１ｍと、成立条件の数と等しい数の誤差ユニット１２１，１２２，・・・，１２ｎと、格子ユニット、誤差ユニットを制御する制御回路１３０とを有している。
【０１５０】
図２３は格子ユニット１１ｉの構成を詳細に示すブロック図であり、誤差ユニット１２ｊから供給される誤差値ＤJを入力として該当する物理量ｕiの修正量Δｕiを算出する修正量算出部１１ｉａと、物理量推定値ｕiを保持する推定値保持部１１ｉｂと、修正量算出部１１ａからの修正量Δｕiの出力を許容する状態と禁止する状態とを選択する選択部１１ｉｃと、選択部１１ｉｃにより出力が許容された修正量Δｕiと推定値保持部１１ｉｂに保持されている推定値ｕiとを加算して再び推定値保持部１１ｉｂにの供給する推定値修正部１１ｉｄとを有している。尚、選択部１１ｉｃにより出力が許容された修正量Δｕiはそのまま誤差ユニットに対する出力として使用される。また、修正量算出部１１ｉａは、誤差信号ＤJに対応して定数δij（＝βiαij）が設定されてあり、数２９の演算を行なって修正量Δｕiを算出する。
【０１５１】
【数２９】
【０１５２】
図２４は誤差ユニット１２ｊの構成を詳細に示すブロック図であり、格子ユニット１１ｉから出力される修正量Δｕiを入力として該当する誤差Ｄjの修正量αijΔｕiを出力する誤差修正量演算部１２ｊａと、誤差値Ｄjを保持する誤差値保持部１２ｊｂと、誤差修正量演算部１２ｊａから出力される誤差修正量αijΔｕiと誤差値保持部１２ｊｂに保持されている誤差値Ｄjとを減算して再び誤差値保持部１２ｊｂに供給する誤差値修正部１２ｊｃとを有している。尚、誤差値保持部１２ｊｂに保持されている誤差値Ｄjはそのまま格子ユニットに対する出力として使用される。上記誤差修正量演算部１２ｊａは、修正量Δｕiに対応して定数αijが設定されてあり、αijΔｕiの演算を行なって誤差修正量を算出する。
【０１５３】
上記制御回路１３０は、選択部１１ｉｃを選択的に動作させる選択制御信号を出力するとともに、推定値保持部１１ｉｂに対して推定値修正部１１ｉｄからの出力値を保持すべきことを示す指示信号を出力し、また、誤差値保持部１２ｊｂに対して誤差値修正部１２ｊｃからの出力値を保持すべきことを示す指示信号を出力するものである。
【０１５４】
したがって、この実施例の場合には、実施例５と同様の作用を達成することができる。そして、上記の説明から明らかなように、メモリとして推定値保持部１１ｉｂと誤差値保持部１２ｊｂとが必要なだけであるから、メモリ容量を著しく低減することができる。即ち、メモリ容量の大幅な低減および物理源解析所要時間の大幅な短縮を両立させることができる。
【０１５５】
【実施例８】
図２５は格子ユニット２１ｉの他の構成例を示すブロック図であり、図２６は楕円中心に対する修正量を算出するための修正量算出ユニット２２の構成を示すブロック図である。尚、誤差修正量演算ユニットの構成は図２３と同一であるから図示および説明を省略する。
【０１５６】
上記格子ユニット２１ｉは、図２５に示すように、図２３の構成とほぼ同一であり、異なる点は、選択部１１ｉｃに代えて、修正量算出ユニット２２から出力される２つの修正量Δｕi，Δｕkおよび修正量算出部１１ｉａから出力される修正量Δｕiの何れかを選択するか、または何れも選択しないことを選択する選択部２１ｉｃを採用した点、および選択部２１ｉｃにより選択された修正量をそのまま修正量算出ユニット２２に対する出力としても採用した点のみである。
【０１５７】
修正量算出ユニット２２は、図２６に示すように、全ての格子ユニット２１ｉから出力される修正量Δｕiを入力とするマルチプレクサ２２ａと、マルチプレクサ２２ａにより選択された２つの修正量Δｕi，Δｕkをそれぞれ保持する１対の修正量保持部２２ｂ，２２ｃと、両修正量保持部２２ｂ，２２ｃにより保持されている修正量が修正方向を規定する基準値として供給され、しかも何れか一方の修正量を入力として数３０の演算を行なって２つの修正量Δｕi，Δｕkを出力する修正量算出部２２ｄとを有している。
【０１５８】
尚、この実施例においては、特には図示していないが、制御回路１３０によりマルチプレクサ２２ａに対して選択制御信号を供給するとともに、両修正量保持部２２ｂ，２２ｃに対して選択的に、マルチプレクサ２２ａから出力される修正値を保持すべきことを示す指示信号を供給するよう構成されている。したがって、マルチプレクサ２２ａにより１つの修正量を選択して一方の修正量保持部に保持させた後、マルチプレクサ２２ａにより他の修正量を選択して他方の修正量保持部に保持させ、次いで、マルチプレクサ２２ａにより交互に上記２つの修正量Δｕi，Δｕkを選択して、実施例６と同様の処理を行なって楕円中心に対応する２つの修正量Δｕi，Δｕkを算出することができる。尚、他の構成部分の作用は実施例７と同様であるから、詳細な説明を省略する。
【０１５９】
【発明の効果】
以上のように請求項１の発明は、磁場源の解析を高速に達成することができ、しかも磁場のｚ方向成分のみを計測することにより、各格子点におけるx方向電流成分、ｙ方向電流成分を解析することができるという特有の効果を奏する。
【図面の簡単な説明】
【図１】この発明の物理量解析方法の一実施例を説明するフローチャートである。
【図２】誤差関数の一例を示す図である。
【図３】分布電流を格子近似した磁場源モデルと磁場センサ配列面との関係を示す概略図である。
【図４】この発明の物理量解析装置の一実施例を示すブロック図である。
【図５】格子ユニットおよび修正ユニットの構成を詳細に示すブロック図である。
【図６】観測ユニットおよび物理場演算ユニットの構成を詳細に示すブロック図である。
【図７】格子ユニットおよび修正ユニットの構成を詳細に示すブロック図である。
【図８】観測ユニットおよび物理場演算ユニットの構成を詳細に示すブロック図である。
【図９】格子ユニットおよび修正ユニットの構成を詳細に示すブロック図である。
【図１０】この発明の他の実施例としての適応ノイズキャンセル方法を説明するフローチャートである。
【図１１】適応ノイズキャンセル方法の具体的適用例を示す波形図である。
【図１２】この発明のさらに他の実施例としての音響探査方法を説明するフローチャートである。
【図１３】音響探査の具体例を示す波形図である。
【図１４】この発明の物理量解析方法のさらに他の実施例を説明するフローチャートである。
【図１５】多数の電極により収集した体表面電位分布データに基づいて体内電気活動状態を可視化するための装置を示す概略図である。
【図１６】図１５の装置を具体化するための処理の一部として電位パターンを得るための物理量解析方法の具体例を説明するフローチャートである。
【図１７】３次元格子状抵抗ネットモデルの一部を拡大して示す図である。
【図１８】仮定された電流源と解析結果とを示す図である。
【図１９】この発明の物理量解析方法のさらに他の実施例を説明するフローチャートである。
【図２０】修正すべき物理量を２つ選択した場合における誤差評価関数の一例を示す図である。
【図２１】図１９のフローチャートに基づく修正量算出の原理を説明する概略図である。
【図２２】この発明の物理量解析装置の他の実施例を示すブロック図である。
【図２３】格子ユニットの構成を詳細に示すブロック図である。
【図２４】誤差修正量演算ユニットの構成を詳細に示すブロック図である。
【図２５】格子ユニットの他の構成例を詳細に示すブロック図である。
【図２６】楕円中心に対する修正量を算出するための修正量算出ユニットの構成を示すブロック図である。
【図２７】従来の適応ノイズキャンセラを示すブロック図である。
【図２８】従来の物理量解析装置を示すブロック図である。
【符号の説明】
１ｉａ 修正ユニット ２ｉ 観測ユニット
３１ｉ 修正量基本値演算セル ３２ｉ 修正値演算セル
３３ｉ 推定値保持セル ４２ｉ 誤差演算セル
１１ｉａ 修正量算出部 １１ｉｂ 推定値保持部
１１ｉｃ 選択部 １１ｉｄ 推定値修正部
１２ｊａ 誤差修正量演算部 １２ｊｂ 誤差値保持部
１２ｊｃ 誤差値修正部
【産業上の利用分野】
この発明は物理量解析方法およびその装置に関し、個々の物理源の物理量ｕiと、個々の物理源に起因して生ずる任意箇所における物理量Ｏjとが数３の関係を有する物理系における未知の物理源の物理量ｕiを解析するための新規な方法およびその装置に関する。
【０００２】
【従来の技術】
従来から、物理量解析方法として、（１）複数個の磁束計の出力に基づいて、スーパーコンピュータを用い、かつモンテカルロ法に基づいて磁場源の解析を行なう方法、（２）適応ノイズキャンセラに適用された方法、（３）インパルス応答推定装置に組み込まれた方法、および（４）ニューラルネットワークを構成するニューロン素子の閾値、結合係数等を系の物理法則に基づいて設定することにより、これらの学習所要時間を不要にした方法（特開平５−９４５４３号公報参照）が提案されている。また、（５）有限要素法も提案されている。
【０００３】
上記（１）の方法は、
ａ）複数個の磁束計による探査空間に乱数を用いてｍ個の電流素片をばらまく。
ｂ）推定誤差演算プロセスにより全推定誤差（推定誤差の総和）を算出する。
ｃ）以下のｄ）からｇ）の処理を反復する。
ｄ）任意に電流素片ｋを選択し、該当する電流素片のパラメータおよび全推定誤差を退避する。
ｅ）電流素片ｋのパラメータを乱数を用いて微小な量だけ変化させる。
ｆ）推定誤差演算プロセスにより全推定誤差を算出する。
ｇ）退避した全推定誤差とｆ）で算出された全推定誤差とを比較し、退避した全推定誤差の方が小さければｄ）で退避した情報を復帰させる。
【０００４】
上記（２）の方法を組み込んだ適応ノイズキャンセラは、図２７に示すように、信号源７１からの情報にノイズ源７２からのノイズが混入した入力Ｓjを誤差演算器７３の非反転入力端子に供給し、ノイズ源７２からのノイズのみをＦＩＲフィルタ７４を介して誤差演算器７３の反転入力端子に供給し、しかも、誤差演算器７３からの出力信号をＦＩＲフィルタ７４にフィードバックしている。また、ＦＩＲフィルタ７４には、ＬＭＳ（Least Mean Square ）アルゴリズムが採用されている。
【０００５】
この適応ノイズキャンセラでは、推定ゲインを適正に設定することにより、情報に混入したノイズのみを除去でき、空調ダクト騒音の消去、自動車の車内の静音化等を達成できる。即ち、除去すべきノイズを高精度に推定できる。
上記（３）の方法を組み込んだインパルス応答推定装置は、高速フーリエ変換（以下、ＦＦＴと略称する）を採用して周波数成分の解析を行ない、解析結果に基づいてインパルス応答の推定を行なう方法である。
【０００６】
上記（４）の方法は、上記（１）の方法の不都合を解消すべく本件特許出願人が既に特許出願を行なった方法であり、図２８に示すように、入力パターンを複数の物理公式演算ユニット８１１，８１２，・・・，８１ｍに供給してそれぞれ既知の物理公式に基づく演算を行ない、全ての物理公式演算ユニット８１１，８１２，・・・，８１ｍからの出力をシグマユニット８２に供給して総和を得、総和および実際の計測値を誤差演算器８３に供給して誤差（両者の差）を得、得られた誤差を全ての物理公式演算ユニット８１１，８１２，・・・，８１ｍの補正部８１１ａ，８１２ａ，・・・，８１ｍａにフィードバックしている。また、物理公式演算ユニット８１１，８１２，・・・，８１ｍにおいて推定されている変数の値を情報収集ユニット８４により収集して解析結果として出力する。
【０００７】
したがって、物理公式自体に関しては何ら学習を行なう必要がなく、物理公式に含まれる変数のみについて誤差および推定ゲインに基づく補正を反復するだけで高精度の解析結果を得ることができる。
上記（５）の方法は、建築、機械の構造解析、気象計算、計算天文学、電磁解析等、多種多様な分野で広く利用されている方法であり、解析対象を有限個のメッシュに分割してモデルを構築し、このモデルに基づいて与えられる連立一次方程式を解くことにより物理源解析結果を得ることができる。そして、連立一次方程式を解くための方法として、ガウスの消去法、ＳＯＲ法、ＣＧ法、ＳＤ法等の種々の方法が知られている。
【０００８】
【発明が解決しようとする課題】
上記（１）の方法は、全推定誤差が小さくなるように電流素片ｋのパラメータを微小量ずつ変化させるのであるから、最終的に正しい解析結果が得られるように思われる。しかし、１回の処理を行なっても全推定誤差が小さくなるという保証が全くなく、しかも上記処理の一部についてのみ並列処理が可能であり、他の部分（ｍ個の電流素片が測定点につくる磁場の演算処理、および全推定誤差を求める処理）については並列処理が不可能であるから、スーパーコンピュータを用いても著しく長時間がかかり、しかも処理回数を増加させても最終的な解を得ることができないという不都合がある。
【０００９】
上記適応ノイズキャンセラは、ＬＭＳアルゴリズムを採用しているのであるから推定ゲインの設定が必要であり、適切な推定ゲインの設定が困難であるという不都合がある。さらに詳細に説明すると、ノイズの伝播経路は著しく多いのであり、これらの全ての伝播経路を通ってノイズキャンセル対象位置に到達する実際のノイズを高精度に推定するためには、ノイズ源７２のノイズを微小時間ずつ遅延させ、それぞれに対して適切な推定ゲインを設定しなければならなくなる。この結果、設定すべき推定ゲインの数が著しく多くなるだけでなく、全ての推定ゲインを適切に設定しなければならないことになる。もし、全ての推定ゲインを同一の値に設定するのであれば、推定処理が発散しないように推定ゲインを小さな値に設定せざるを得ず、これに伴なって、解の収束が遅くなってしまうという不都合もある。
【００１０】
上記（３）の方法は、ＦＦＴを採用しているのであるから次の不都合が生じる。ＦＦＴはサンプリング定理に基づく処理を行なう方法であるから、通常測定信号に含まれている不要な高調波を除去するためにアンチエリアシングフィルタと呼ばれるローパスフィルタを設けることが必須であり、構成が複雑化するという不都合がある。また、サンプリング区間のデータが周期的に連続するという補償がなければ周波数解析結果の精度が著しく低下してしまうのであるから、適用可能な信号の種類が制限されてしまうという不都合もある。このような不都合を解消するために、ハミング、ハニング等の窓関数を用いることが提案されているが、逆フィルタの演算時に窓関数を用いると演算後にサンプリング区間全体にわたって波形が歪み、解析精度が低下するという新たな不都合が生じる。さらに、サンプリング間隔に基づいて定まる周波数間隔の出力しか得られないのであるから、広帯域の解析が必要な場合には必然的にサンプル数を増加させなければならないという不都合もある。さらにまた、周波数軸が対数目盛であっても等間隔にサンプリングを行なわなければならないので著しく多量のメモリが必要になり、また、サンプル数も２Nでなければならないという制約があるという不都合もある。
【００１１】
上記（４）の方法は、物理公式演算結果の偏微分値と誤差との積の総和に推定ゲインを乗算して得た値で補正を行なうのであるから、推定ゲインの設定が必要である。そして、設定すべき推定ゲインが多くなるとともに、適切な推定ゲインを設定しなければならないのであるから、推定ゲインの設定が著しく繁雑になってしまうという不都合がある。簡単化のために、もし、全ての推定ゲインを同一の値に設定するのであれば、推定処理が発散しないように推定ゲインを小さな値に設定せざるを得ず、これに伴なって解の収束が遅くなってしまうという不都合もある。
【００１２】
上記（５）の方法は、メッシュ分割数を増加させることにより物理量の解析精度を高めることができる。
しかし、連立一次方程式を解くためにガウスの消去法を採用すると、未知数の数がｎである場合に、ｎ3／３のオーダーの演算回数が必要であり、著しく演算負荷が大きくなるので、未知数の数が多い場合には、スーパーコンピュータを用いても実用上十分な速度での物理源解析を達成することができない。これに対してコレスキー法を採用すれば、演算回数をｎ3／６のオーダーにすることができるが、サイズが大きい連立一次方程式を解くための方法としては演算負荷が大きすぎる。そして、何れの方法においても誤差が累積されることになるので、物理源解析精度を余り高めることができない。
【００１３】
上記ＳＯＲ法は、ガウス＝ザイデル法の収束を加速するために、加速パラメータを導入して過大修正を行なう方法であるから、所定の条件を満足する係数マトリクスに関しては解の収束が保証されているが、所定の条件を満足しない係数マトリクスに関しては解の収束が保証されないという不都合がある。そして、解の収束が保証される係数マトリクスであっても、加速パラメータの最適値を簡単にかつ確実に求める方法がなく、加速パラメータの設定が不適当な場合には、収束を十分に加速することができない。
【００１４】
上記ＳＤ法は、最急降下法と呼ばれる方法であり、ＳＯＲ法のような係数マトリクスの制約はないが、反復回数が著しく多く、しかも最急降下方向の算出、未知数を補正するためのゲインの算出を毎回行なう必要があるとともに、これらの演算はベクトル演算および除算を含むのであるから、全体として演算負荷が著しく大きくなってしまう。
【００１５】
上記ＣＧ法は、共役傾斜法と呼ばれる方法であり、既に修正した方向の全てと直交する方向を修正方向として順次未知数を補正するのであるから、理論的には反復回数を著しく少なくすることができ（ｎ回以下にすることができ）、しかも物理源解析精度を著しく高くすることができる。しかし、実際には、殆ど全ての補正処理において誤差が発生するのであるから、係数マトリクスによって反復回数が大きくばらついてしまう。また、各回の処理を行なうに当ってベクトル演算、マトリクス演算が必須になるのであるから、全体として演算負荷が著しく大きくなってしまう。
【００１６】
【発明の目的】
この発明は上記の問題点に鑑みてなされたものであり、推定ゲインの設定が不要であり、しかも係数マトリクス、サンプル数の制約を排除して高速かつ高精度に物理量の解析を行なうことができる物理量解析方法およびその装置を提供することを目的としている。
【００１７】
【課題を解決するための手段】
上記の目的を達成するための、請求項１の物理量解析装置は、測定対象物が発生する磁場を複数の磁場センサで計測し、前記測定対象物内部の磁場源における物理量を解析する装置であって、
（１） 個々の磁場源ｉに対応する単位量の磁場源が各所定箇所に発生する磁場を示す、数１で与えられる定数αijを予め定めておく手段と、
（２） 物の内部に設定された格子点ｋ（ｋ＝１，２・・・，ｐ）に位置する各磁場源の、数２で示すx方向電流成分およびｙ方向電流成分からなる物理量ｕiを仮に定めておく手段と、
（３） 数３の演算を行なって観測点ｊに発生すると予測される磁場Ｏjを算出する手段と、
（４） 観測点ｊにおける磁場の実測値Ｓjを測定する手段と、
（５） 観測点ｊにおける磁場の実測値Ｓjと上記磁場Ｏjとの差Ｓj−Ｏjを算出する手段と、
（６） 数４の演算を行なって仮の解としての物理量ｕiを得る手段と、
（７） （３）から（６）までの各手段による処理を全ての磁場源ｉについて行なわせる手段と、
（８） 上記差Ｓj−Ｏjの二乗の総和が所定の閾値よりも小さくなるまで、（３）から（６）までの各手段による処理の反復処理を反復させる手段と、
（９） 最終的に得られた各磁場源の物理量ｕiを磁場源の物理量ｕiとして採用する手段とを含むものである。
【００４６】
請求項１の物理量解析装置であれば、測定対象物が発生する磁場を複数の磁場センサで計測し、前記測定対象物内部の磁場源における物理量を解析する装置であって、
（１） 個々の磁場源ｉに対応する単位量の磁場源が各所定箇所に発生する磁場を示す、数１で与えられる定数αijを予め定めておく手段と、
（２） 物の内部に設定された格子点ｋ（ｋ＝１，２・・・，ｐ）に位置する各磁場源の、数２で示すx方向電流成分およびｙ方向電流成分からなる物理量ｕiを仮に定めておく手段と、
（３） 数３の演算を行なって観測点ｊに発生すると予測される磁場Ｏjを算出する手段と、
（４） 観測点ｊにおける磁場の実測値Ｓjを測定する手段と、
（５） 観測点ｊにおける磁場の実測値Ｓjと上記磁場Ｏjとの差Ｓj−Ｏjを算出する手段と、
（６） 数４の演算を行なって仮の解としての物理量ｕiを得る手段と、
（７） （３）から（６）までの各手段による処理を全ての磁場源ｉについて行なわせる手段と、
（８） 上記差Ｓj−Ｏjの二乗の総和が所定の閾値よりも小さくなるまで、（３）から（６）までの各手段による処理の反復処理を反復させる手段と、
（９） 最終的に得られた各磁場源の物理量ｕiを磁場源の物理量ｕiとして採用する手段とを含むものであるから、推定ゲインの設定が不要になり、推定ゲインがないことに伴って解の収束を高速化できる。即ち、磁場源の解析を高速に達成することができる。さらに、磁場のｚ方向成分のみを計測することにより、各格子点におけるx方向電流成分、ｙ方向電流成分を解析することができる。
【００７３】
【実施例】
以下、実施例を示す添付図面によって詳細に説明する。
図１はこの発明の物理量解析方法の一実施例を説明するフローチャートであり、複数（ｍ個）の微小物理源の位置および複数（ｎ個）の観測点の位置がそれぞれ設定され、かつ系に線形加算性が成立する場合に、ステップＳＰ１において各微小物理源ｉ（ｉ＝１，２，・・・，ｍ）と各観測点ｊ（ｊ＝１，２，・・・，ｎ）の位置的関係から物理法則を用いて定数αｉｊを算出し、ステップＳＰ２において各微小物理源ｉの物理量ｕｉを仮に設定する（例えば、全ての物理量を０に設定する）。そして、ステップＳＰ３において数１の演算を行なうことにより観測点ｊに発生すると予測される物理場Ｏｊを算出し、ステップＳＰ４において観測点ｊにおける物理場の実測値Ｓｊを得る。
【００７４】
次いで、ステップＳＰ５において、実測値Ｓｊと物理場Ｏｊとの差および定数αｉ jを用いて数４の演算を行なうことにより仮の解としての物理量ｕｉを得る。
尚、ステップＳＰ３，ＳＰ４，ＳＰ５の処理は全ての微小物理源ｉについて順次行なわれる。
その後、ステップＳＰ６において実測値Ｓｊと物理場Ｏｊとの差の二乗の総和が十分に小さくなったか否かを判別し、十分には小さくなっていないと判別された場合には、再びステップＳＰ３の処理を行なう。逆に、十分に小さくなったと判別された場合には、ステップＳＰ７において、最終的に得られている物理量ｕｉを解析結果として取り出し、そのまま一連の処理を終了する。
【００７５】
さらに詳細に説明すると、誤差関数Ｅとして数５に示すように、実測値Ｓｊと物理場Ｏｊとの差の二乗の総和を採用する。
【００７６】
【数５】
ここで、
Ｓｊ−Ｏｊ＝（Ｓｊ−Ｏｊ＋αｉｊｕｉ）−αｉｊｕｉ
であるから、未知数の物理量ｕｉを含まない項と未知の物理量ｕｉを含む項とに分離できる。したがって、数５は数６と等価である。
【００７８】
【数６】
但し、Ａｉ ２，ＡｉＢｉ，Ｃiはそれぞれ数７，数８，数９で与えられる。
【００８０】
【数７】
【数８】
【数９】
数６から明らかなように、誤差関数Ｅは，図２に示すように、ｕｉに関して軸をＢｉ/Ａｉ、最小値をＣｉ−Ｂｉ ２とする放物線になる。即ち、誤差関数Ｅは，全体としてｍ次元放物面体であり、一意収束性がある。したがって、ｕｉの軸値への修正（ｕｉ←Ｂｉ/Ａｉ）は解への漸近を意味することになる。また、上記修正式に数７、数８を代入して整理することにより数４が得られる。
【００８４】
図３は分布電流を格子近似した磁場源モデルと磁場センサ配列面との関係を示す概略図であり、磁場センサにより対象物から発生する磁場を検出し、検出された磁場に基づいて対象物の内部の電流分布を測定する場合における磁場源モデルと磁場センサ配列面との関係を示している。図３においては、ｎ個の磁場センサが配列されているとともに、ｍ個の磁場源が想定されている。そして、磁場のｚ方向成分のみの計測を行なうと仮定すれば、磁場源はｘ方向とｙ方向の電流成分を有することになる。また、ｐ個の３次元格子で関心領域を近似すれば、各格子座標にｘ方向とｙ方向の２個の独立した電流双極子を割り当てることになる。即ち、ｍ＝２ｐとなる。
【００８５】
格子点ｋ（ｋ＝１，２，・・・，Ｐ）のｘ方向電流成分をＰｘｋ、ｙ方向電流成分をＰｙｋとし、その座標を（Ｘｋ，Ｙｋ，Ｚｋ）とし、観測点ｊの座標を（ｘｊ，ｙｊ，ｚｊ）とする。そして、磁場センサのピックアップコイル径を無視し、コイル中心の磁束密度を計測していると仮定すれば、数２の下で、ビオサバールの法則により、αｉｊは数１となる。そして、このように設定することにより、数３と同様に線形加算が成立する系になるので、前記と同様の処理を行なうことができ、磁場源の内部の電流分布を算出することができる。
【００８６】
そして、数４による推定を安定化するために、事前に正規化（αｉｊ’＝αｉｊ/Ａｉ）を行なって各放物線の開口度を１に統一する。また、正規化後の定数αｉｊ’による推定値を補正するためにｕｉ＝ｕｉ’/Ａｉの演算を行なう。
以上の一連の解析処理を行なうことにより、各格子点における電流成分を高速かつ正確に推定できた。
【００８７】
この実施例は上記具体例に限定されるものではなく、例えば、集光特性が既知の光量センサを採用することにより光量解析方法に適用することが可能であり、その他、指向性を持ち、かつ指向性が既知のセンサを採用することにより種々の解析方法に適用することが可能である。
【００８８】
【実施例２】
図４はこの発明の物理量解析装置の一実施例を示すブロック図であり、解析対象となる物理源の数と等しい数の格子ユニット１１，１２，・・・，１ｍと、各格子ユニット１１，１２，・・・，１ｍに対応する修正ユニット１１ａ，１２ａ，・・・，１ｍａとからなる格子ユニット層１と、観測点の数と等しい数の観測ユニット２１，２２，・・・，２ｎと、各観測ユニット２１，２２，・・・，２ｎに対応する物理場演算ユニット２１ａ，２２ａ，・・・，２ｎａとからなる観測ユニット層２と、格子ユニット層１、観測ユニット層２を制御する制御回路３とを有している。ここで、制御回路３は、全ての観測ユニット２１，２２，・・・，２ｎに対して互に同じタイミングで格納指示信号を供給し、各格子ユニット１１，１２，・・・，１ｍに対して順次修正指示信号を供給するものである。
【００８９】
図５は格子ユニット１ｉおよび修正ユニット１ｉａの構成を詳細に示すブロック図であり、修正量基本値演算セル３１ｉと、修正値演算セル３２ｉと、推定値保持セル３３ｉとを有している。
上記修正量基本値演算セル３１ｉは、ｎ個の入力端子を有しているとともに、各入力端子にそれぞれ定数乗算機能を持たせており、各観測ユニット２１，２２、・・・、２ｎからの出力値（Ｓｊ−Ｏｊ）を対応する入力端子に供給している。そして、定数αｉｊおよび出力値（Ｓｊ−Ｏｊ）に基づいて数１０で示す修正量基本値を出力する。
【００９０】
【数１０】
上記修正値演算セル３２ｉは、２個の入力端子を有しているとともに、各入力端子にそれぞれ１、数１１の定数を乗算する定数乗算機能を持たせており、１の定数を乗算する入力端子に推定値保持セル３３ｉから出力される今回の推定前の推定値ｕｉを供給しているとともに、他方の入力端子に修正量基本値を供給している。そして、定数および供給された値に基づいて数４で示す修正値を得、推定値保持セル３３ｉに供給する。
【００９２】
【数１１】
上記推定値保持セル３３ｉは、制御回路３から修正指示信号ｉが供給されたことに応答して、修正値演算セル３２ｉから供給される修正値を新たに保持するものである。
尚、例えば、上記修正量基本値演算セル３１ｉとしては、ニューロンデバイスを採用することが可能であり、しかもニューロンデバイスの飽和関数のうちリニアな部分のみを使用すれば足りるのであるから、ニューロンデバイスに代えて、乗算器、加算器のみを用いて修正量基本値演算セル３１ｉを構成することが可能である。
【００９４】
図６は観測ユニット２ｉおよび物理場演算ユニット２ｉａの構成を詳細に示すブロック図であり、物理場演算セル４１ｉと、誤差演算セル４２ｉと、観測値保持セル４３ｉとを有している。
上記物理場演算セル４１ｉは、ｍ個の入力端子を有しているとともに、各入力端子にそれぞれ定数乗算機能を持たせており、各格子ユニット１１，１２，・・・１ｍからの出力値ｕｉを対応する入力端子にを供給している。そして、定数αｉｊおよび出力値ｕｉに基づいて数３で示す物理場演算値を出力する。
【００９５】
上記観測値保持セル４３ｉは、制御回路３から格納指示信号ｊが供給されたことに応答して、観測値を新たに保持するものである。
上記誤差演算セル４２ｉは、所定の物理量を検出するセンサが検出した観測値から物理場演算値を減算して得た値を誤差として出力するものである。
したがって、制御回路３により修正指示信号ｉおよび格納指示信号ｊを出力することにより図１のフローチャートと同様に高速かつ高精度に物理量の解析を行なうことができる。
【００９６】
図７は、格子ユニット１ｉおよび修正ユニット１ｉａの構成を詳細に示すブロック図であり、図５の構成例と異なる点は、修正値演算セル３２ｉを省略し、修正量基本値演算セル３１ｉに、推定値を入力とし、かつ１の定数に基づく乗算機能を有する入力端子を追加するとともに、他の各入力端子に、図５の各入力端子が定数αｉｊに対して数１１を乗算して得た定数を乗算機能として持たせて新たな修正値演算セルを構成した点のみである。
【００９７】
したがって、この場合には、構成を簡素化でき、しかも図５と同様の作用を達成できる。
図８は観測ユニット２ｉおよび物理場演算ユニット２ｉａの構成を詳細に示すブロック図であり、図６の構成例と異なる点は、誤差演算セル４２ｉを省略し、物理場演算セル４１ｉの各入力端子に、図６の各入力端子が定数αijに対して−１を乗算して得た定数を乗算機能として持たせるとともに、観測値を入力とし、かつ１の定数に基づく乗算機能を有する入力端子を追加して物理場誤差演算セルを構成した点のみである。
【００９８】
したがって、この場合にも、構成を簡素化でき、しかも図６と同様の作用を達成できる。
図９は格子ユニット１ｉおよび修正ユニット１ｉａの構成を詳細に示すブロック図であり、図７の構成例と異なる点は、順次供給される誤差を保持する誤差格納セル（例えば、サンプルホールド回路）３４ｉ１，３４ｉ２，・・・，３４ｉｎをさらに有している点および推定値ｕiに対して定数αijを乗算する乗算セル３５ｉをさらに有している点のみである。
【００９９】
したがって、この場合には、観測ユニットを１つのみにすることにより構成を簡素化でき、しかも図７と同様の作用を達成できる。
【０１００】
【実施例３】
図１０はこの発明の他の実施例としての適応ノイズキャンセル方法を説明するフローチャートであり、ステップＳＰ１において、現時点を基準とする過去のノイズのリファレンスαj-kを得、ステップＳＰ２において各ノイズの実測値に対する相関係数（混ざり具合を示す値）ｕｉを仮に設定する。そしてステップＳＰ３において数１２の演算を行なうことにより観測点ｊに発生すると予測されるノイズＯｊを算出し、ステップＳＰ４において観測点ｊにおける実測値Ｓｊを得る。
【０１０１】
【数１２】
次いで、ステップＳＰ５において、実測値Ｓｊと予測ノイズＯｊとの差およびリファレンスαｋ−ｊを用いて数１３の演算を行なうことにより仮の解としての相関係数ｕｉを得る。
【０１０３】
【数１３】
尚、ステップＳＰ３，ＳＰ４，ＳＰ５の処理は全ての時点のノイズｉについて順次行なわれる。
その後、ステップＳＰ６において実測値Ｓjと予測ノイズＯjとの差の二乗の総和が十分に小さくなったか否かを判別し、十分には小さくなっていないと判別された場合には、再びステップＳＰ３の処理を行なう。逆に、十分に小さくなったと判別された場合には、ステップＳＰ７において、最終的に得られている相関係数ｕiを解析結果として採用し、ステップＳＰ８において、情報とノイズが混在した信号からリファレンスαijおよび相関係数ｕiに基づいて得られる予測ノイズＯj を減算して情報のみを抽出し、そのまま一連の処理を終了する。
【０１０５】
図１１はこの適応ノイズキャンセル方法の具体的適用例を示す波形図であり、５０Ｈｚの交流成分が重畳された心電図｛図１１（Ａ）参照｝から交流成分を除去した場合｛図１１（Ｂ）参照｝を示している。
図１１から明らかなように、従来の適応ノイズキャンセラのように推定ゲイン（ステップサイズとも呼ばれる）を設定する必要がなく、しかも高精度にノイズを除去して正確な心電図が得られる。
【０１０６】
【実施例４】
図１２はこの発明のさらに他の実施例としての音響探査方法を説明するフローチャートであり、ステップＳＰ１において、各サンプリング時点における送信音波αｋを得、ステップＳＰ２において探査対象となる物理源のインパルス応答ｕｋを仮に設定する。そして、ステップＳＰ３において数１４の演算を行なうことにより観測点ｊに発生すると予測される予測受信音波Ｏｊを算出し、ステップＳＰ４において観測点ｊにおける実測値Ｓｊを得る。但し、ｍは送信音波の１周期におけるサンプル数、ｎは受信音波の総サンプル数である。
【０１０７】
【数１４】
次いで、ステップＳＰ５において、実測値Ｓｊと予測受信音波Ｏｊとの差および送信音波αｋを用いて数１５の演算を行なうことにより仮の解としての相関係数ｕｋを得る。
【０１０９】
【数１５】
尚、ステップＳＰ３，ＳＰ４，ＳＰ５の処理は該当する全ての時点の送信音波αkについて順次行なわれる。
その後、ステップＳＰ６において実測値Ｓjと予測ノイズＯjとの差の二乗の総和が十分に小さくなったか否かを判別し、十分には小さくなっていないと判別された場合には、再びステップＳＰ３の処理を行なう。逆に、十分に小さくなったと判別された場合には、ステップＳＰ７において、最終的に得られているインパルス応答ｕk解析結果として採用し、そのまま一連の処理を終了する。
【０１１１】
尚、数１４、数１５の上段においてｍ個の総和を算出しているのは、送信音波が１周期分のみであると仮定しているからであり、数１５の下段において（ｎ-ｋ）個の総和を算出しているのは、受信音波の観測時間が有限であり、送信音波の一部のみに対応する受信音波しか得られていないからである。
図１３は音響探査の具体例を示す波形図であり、図１３（Ａ）に示す送信音波に基づいて図１３（Ｂ）に示す受信音波が得られた場合に、ノイズが存在していなければ図１３（Ｃ）に示すインパルス応答推定結果を得ることができ、ノイズ（１％）が存在している場合に図１３（Ｄ）に示すインパルス応答推定結果を得ることができた。
【０１１２】
図１３から明らかなように、ノイズが存在していなければ著しく高精度にインパルス応答を推定でき、ある程度のノイズが存在していてもかなり高い精度でインパルス応答を推定できることが分る。
尚、この実施例は音波のみについて説明したが、電磁波等にも適用できることはもちろんである。
【０１１３】
実施例３、実施例４の各方法は、実施例１の定数に代えて、解析目的に適合する値を採用し、推定値も実施例１の推定値と異なるのであるが、採用された値、推定値に基づく処理は実施例１とほぼ同様であるから、実施例２の装置に適用することにより簡単に実施例３、実施例４に対応する装置を得ることができる。但し、適応ノイズキャンセラに適用する場合には、各観測ユニット２１，２２，・・・，２ｎにそれぞれ別個の観測値を供給する代わりに、観測値が１つの観測ユニットに供給され、各観測ユニットの値が順次隣合う観測ユニットに供給されるよう構成すればよい。
【０１１４】
【実施例５】
図１４はこの発明の物理量解析方法のさらに他の実施例を説明するフローチャートであり、個々の物理源を含む、線形連立一次方程式が成立する領域内において、個々の物理源の物理量ｕiおよび上記領域に基づいて定まる比例定数αijに基づいて任意箇所で得られる物理量Ｏjを解析する方法を説明している。
【０１１５】
ステップＳＰ１において、領域の性質等に基づいて比例定数αijを得、ステップＳＰ２において各物理源の物理量ｕｉを仮に設定する（例えば、全ての物理量を０に設定する）。そして、ステップＳＰ３において数３の演算を行なうことにより任意箇所で得られる物理量Ｏｊを算出し、ステップＳＰ４において成立条件（領域内に与えた既知の物理量）Ｓｊと物理量Ｏｊとの差を算出する。
【０１１６】
次いで、ステップＳＰ５において上記差および比例定数αijを用いて数４の演算を行なうことにより仮の解としての物理量ｕｉを得る。
尚、ステップＳＰ３，ＳＰ４，ＳＰ５の処理は全ての物理源について順次行なわれる。
その後、ステップＳＰ６において成立条件Ｓｊと物理量Ｏｊとの差の二乗の総和が十分に小さくなったか否かを判別し、十分には小さくなっていないと判別された場合には、再びステップＳＰ３の処理を行なう。逆に、十分に小さくなったと判別された場合には、ステップＳＰ７において、最終的に得られている物理量ｕｉを解析結果として取り出し、そのまま一連の処理を終了する。
【０１１７】
さらに詳細に説明すると、誤差関数Ｅとして数５に示すように成立条件Ｓｊと物理量Ｏｊとの差の二乗の総和を採用すれば、数６から数９の関係が得られる。したがって、ｕｉの軸値への修正により解への漸近を達成できるとともに、収束を高速化できる。また、ｕｉの軸値への修正に当って、数７に基づく除算が必要であるが、数７の値は事前に算出可能な定数であるから、予め数７の逆数を算出しておくことにより、上記一連の処理における除算回数を著しく低減することができ、ひいては解析所要時間を一層短縮することができる。
【０１１８】
しかし、一連の処理を行なう毎に数３の演算を行なう必要があり、演算負荷を十分には低減できていない。この点に関連して、本件発明者は、物理量ｕｉを軸値に修正するための修正量△ｕｉ、対応する比例定数αij、および先行する物理量Ｏｊに基づいて数１６の演算を行なうことにより、物理量Ｏｊを修正することができることを見出した。
【０１１９】
【数１６】
したがって、一連の処理を行なう毎に数３の演算を行なう代わりに数１６の演算を行なうことにより、演算負荷を大幅に低減し、ひいては解析所要時間をさらに短縮することができる。
さらに、予めＤｊ＝Ｓｊ-Ｏｊを算出しておき、数１６の演算を行なう代わりに数１７の演算を行なうようにすれば、Ｓｊ、Ｏｊをそれぞれ保持しておく必要がなくなり、メモリ容量を低減することができる。
【０１２１】
【数１７】
この実施例を有限要素法に適用すれば、一連の処理を行なう毎に必要な演算負荷についても上述のように著しく低減することができる。したがって、係数マトリクスの次数が著しく大きい解析対象領域についても簡単に、かつ高速に物理量の解析を行なうことができる。
図１５は多数の電極により収集した体表面電位分布データに基づいて体内電気活動状態を可視化するための装置を概略的に示す図であり、図１６は図１５の装置を実現する処理の一部として電位パターンを得るための物理量解析方法の具体例を説明するフローチャートである。
【０１２３】
図１５の装置は、人体の頭部、胸部等に被着される基部部材２０１の所定位置に所定の配置パターンで電位計２０２が装着されてあり、電位計２０２からの出力信号が生体計測用絶縁アンプ２０３を介してＡ／Ｄコンバータ２０４に導かれてディジタル信号に変換され、ディジタル信号が処理部２０５に供給される。この処理部２０５は電位パターン保持部２０６を有しており、処理部２０５において図１６のフローチャートの処理を行なって得た電位パターンを電位パターン保持部２０６に保持させる。また、電位パターンが得られた後に、電位パターンに基づいて人体内部の電流の解析を行なう電流解析部２０７と、電流解析部２０７による解析結果を表示するための表示部２０８とを有している。
【０１２４】
上記電位計２０２は、例えば人体の表面に圧接できるように付勢機構を有しており、しかも、人体に対する圧接状態において、基部部材２０１に対する相対位置を検出できるように構成してあることが好ましい。
上記電流解析部２０７は、得られた電位パターンを比例定数αijとして採用し、例えば、実施例１と同様の処理を行なって人体内部の電流の解析を行なうものである。
【０１２５】
また、図１６のフローチャートに基づく処理を行なう前に、予めＭＲＩ、Ｘ線ＣＴ等により体内臓器の３次元形状を抽出し、別な計測法により得られ、または医学上の既知情報に基づいて得られた電気抵抗率を用いて、解析対象となる３次元領域を所定ピッチ（例えば５ｍｍピッチ）の３次元格子状抵抗ネットワークでモデル化しておく。また、上記装置は体表面電位分布データに基づいて体内電気活動状態を可視化するものであるから、その前提として、上記３次元領域内の任意の抵抗の一方の端部に単位電流を与えた状態における表面電位分布を、上記３次元領域内に含まれる全ての抵抗について算出しておくことが必要になる。ここで、３次元格子状抵抗ネットモデルのx、y、z方向のノード数をそれぞれＮｘ、Ｎｙ、Ｎｚで表し、図１７に示すように、ノード[ｋx，ｋy，ｋz]の電位をｕ[ｋx，ｋy，ｋz]で表し、ノード[ｋx，ｋy，ｋz]のｘ方向に接続されている抵抗の値をＲx[ｋx−１，ｋy，ｋz]、Ｒx[ｋx，ｋy，ｋz]で表し、ノード[ｋx，ｋy，ｋz]のｙ方向に接続されている抵抗の値をＲy[ｋx，ｋy-１，ｋz]、Ｒy[ｋx，ｋy，ｋz]で表し、ノード[ｋx，ｋy，ｋz]のｚ方向に接続されている抵抗の値をＲz[ｋx，ｋy，ｋz−１]、Ｒz[ｋx，ｋy，ｋz]で表し、外部からノード[ｋx，ｋy，ｋz]に強制的に供給される電流をＳ[ｋx，ｋy，ｋz]で表し、ノード[ｋx，ｋy，ｋz]から接続されている抵抗を通じて流出する電流をＯ[ｋx，ｋy，ｋz]で表す。そして、各抵抗値Ｒx[ｋx−１，ｋy，ｋz]、Ｒx[ｋx，ｋy，ｋz]、Ｒy[ｋx，ｋy−１，ｋz]、Ｒy[ｋx，ｋy，ｋz]、Ｒz[ｋx，ｋy，ｋz−１]、Ｒz[ｋx，ｋy，ｋz]の逆数をそれぞれＧx[ｋx−１，ｋy，ｋz]、Ｇx[ｋx，ｋy，ｋz]、Ｇy[ｋx，ｋy−１，ｋz]、Ｇy[ｋx，ｋy，ｋz]、Ｇz[ｋx，ｋy，ｋz−１]、Ｇz[ｋx，ｋy，ｋz]で表し、実際に接続されている抵抗に対応する総和ＧN[ｋx，ｋy，ｋz]を数１８で表す。
【数１８】
但し、数１８の右辺においては実際に接続されている抵抗に対応する項のみを加算の対象とする。
したがって、オームの法則に基づいて、流出電流Ｏ[ｋx，ｋy，ｋz]は数１９により算出することができる。
【０１２７】
【数１９】
また、キルヒホッフの第１法則を成立させるためには、全てのノードにおいて供給電流Ｓ[ｋx，ｋy，ｋz]と流出電流Ｏ[ｋx，ｋy，ｋz]とが等しくなるノード電位ｕ[ｋx，ｋy，ｋz]を求めればよい。尚、キルヒホッフの第２法則に関しては、未知の物理量をノード電位ｕ[ｋx，ｋy，ｋz]とすることにより自動的に満足される。即ち、数２０で示される評価関数Ｅの値を最小にするノード電位ｕ[ｋx，ｋy，ｋz]を算出すればよいことになる。
【０１２９】
【数２０】
上記の条件下において、ステップＳＰ１において（Ｎx−１）×Ｎy×Ｎz個のＧx[ｋx，ｋy，ｋz]、Ｎx×（Ｎy−１）×Ｎz個のＧy[ｋx，ｋy，ｋz]、Ｎx×Ｎy×（Ｎz−１）個のＧｚ[ｋx，ｋy，ｋz]、Ｎx×Ｎy×Ｎz個のＧN[ｋx，ｋy，ｋz]、およびＮx×Ｎy×Ｎz個の、数２１の関係を有するβ[ｋx，ｋy，ｋz]をそれぞれ算出して格納する。但し、数２１の右辺の各項は実際の接続を反映するように求められている。したがって、解析対象となるノードの数が著しく多くなっても、マトリクスの列方向の非０の比例定数の数は最大でも７にしかならない。したがって、係数マトリクスを保持しておくためのメモリ容量を、マトリクスの次数に比して十分に少なくすることができる。
【数２１】
そして、ステップＳＰ２においてＮx×Ｎy×Ｎz個の推定誤差Ｄ[ｋx，ｋy，ｋz]＝Ｓ[ｋx，ｋy，ｋz]−Ｏ[ｋx，ｋy，ｋz]を算出して格納する。
次いで、ステップＳＰ３において数２２に基づいて未知のノード電位ｕ[ｋx，ｋy，ｋz]を軸値に修正するための修正量Δｕを算出してノード電位ｕ[ｋx，ｋy，ｋz]を修正する。
【０１３２】
【数２２】
そして、ステップＳＰ４において数２３の演算を行なって推定誤差を修正する。
その後、ステップＳＰ５においてｋzがＮzに達したか否かを判別し、達していなければステップＳＰ６においてｋzをインクリメントし、再び、ステップＳＰ３の処理を行なう。逆に、ステップＳＰ５においてｋz＝Ｎzであると判別された場合には、ステップＳＰ７においてｋyがＮyに達したか否かを判別し、達していなければステップＳＰ８においてｋyをインクリメントし、再びステップＳＰ３の処理を行なう。逆に、ステップＳＰ７においてｋy＝Ｎyであると判別された場合には、ステップＳＰ９においてｋxがＮxに達したか否かを判別し、達していなければステップＳＰ１０においてｋxをインクリメントし、再びステップＳＰ３の処理を行なう。逆に、ステップＳＰ９においてｋx＝Ｎxであると判別された場合には、ステップＳＰ１１において評価関数Ｅの値が所定の閾値よりも小さいか否かを判別し、評価関数Ｅの値が所定の閾値よりも小さければステップＳＰ１２において、最終的に得られたノード電位ｕ[ｋx，ｋy，ｋz]を解析結果として採用し、そのまま一連の処理を終了する。逆に、評価関数Ｅの値が所定の閾値以上であれば、ｋx，ｋy，ｋzを初期値にリセットして再びステップＳＰ３の処理を行なう。
【数２３】
以上の一連の処理を行なうことにより、人体の頭部、胴体を解析対象領域とする未知の電位の解析を達成することができた。具体的には、人体の頭部を解析対象とする場合には未知の電位の数が約４００００であり、胴体を解析対象とする場合には未知の電位の数が約４８００００であるが、主記憶が６４Ｍバイトのワークステーションを用いて十分な精度の解析結果を得ることができた。また、所要時間に関しても、前者の場合に３０００回の反復処理を行なって１４４秒、後者の場合に３０００回の反復処理を行なって１９８６秒であった。
【０１３５】
以上のようにして得られた電位パターンを比例定数αijとして採用し、体表面電位分布データを既知の物理量Ｓjとして採用し、実施例１と同様の処理を行なって生体内部に設定された格子点の電流の流入流出状況を解析することができる。
さらに詳細に説明する。
【０１３６】
格子点［１，１，１］に単位電流が流入し、格子点［２，１，１］から電流が流出すると仮定して得た体表面電位分布データを、（格子点数−１）×３列の比例定数αijのうち、第１列の比例定数として採用し、格子点［１，１，１］に単位電流が流入し、格子点［１，２，１］から電流が流出すると仮定して得た体表面電位分布データを第２列の比例定数として採用し、第１列の比例定数として採用し、格子点［１，１，１］に単位電流が流入し、格子点［１，１，２］から電流が流出すると仮定して得た体表面電位分布データを第３列の比例定数として採用し、以下、同様にして他の全ての列の比例定数を採用する。
【０１３７】
このようにして全ての列の比例定数αijを採用すれば、数３の関係を満足することになるので、実施例１と同様の解析処理を行なうことにより、生体内部に設定された格子点の電流の流入流出状況を解析することができる。
具体的には、図１８中の（Ａ）の２次元等電位線図に示される電流源を仮定して図１６の処理を行なって電位分布データを得、この電位分布データを比例定数として採用して実施例１と同様の解析処理を行なった結果、図１８の（Ｂ）に示す電流源解析結果を得ることができた。
【０１３８】
尚、生体内部の電流分布を解析するためには、生体内部に設定された格子点の電流の流入流出状況の解析結果に基づいて図１６のフローチャートと同様の解析処理を行なう必要がある。但し、当初に既知の物理量として生体内部の供給電流Ｓ［ｋx，ｋy，ｋz］および流出電流Ｏ［ｋx，ｋy，ｋz］を採用する代わりに、抵抗を流れる電流を採用しておけば、図１６のフローチャートと同様の解析処理を行なって得た電位分布データを比例定数として採用して実施例１と同様の解析処理を行なうことにより、生体内部の電流分布を得ることができる。
【０１３９】
【実施例６】
図１９はこの発明の物理量解析方法のさらに他の実施例を説明するフローチャートであり、ステップＳＰ１において、解析対象となる多数の物理量の中から任意の２つの物理量を選択し、ステップＳＰ２において、選択した２つの物理量の一方（以下、物理量ｕkと称する）を対象として、実施例５と同様の方法により物理量ｕkを１回だけ修正し、ステップＳＰ３において他方の物理量（以下、物理量ｕiと称する）を対象として、実施例５と同様の方法により物理量ｕiを１回だけ修正する。次いで、ステップＳＰ４，ＳＰ５において、ステップＳＰ２，ＳＰ３と同様にして再び物理量ｕk，ｕiを修正する。その後、ステップＳＰ６において、ステップＳＰ３，ＳＰ４，ＳＰ５の処理により得られた修正量に基づいて、両物理量ｕk，ｕiに基づいて定まる評価関数Ｅの値を最小にするために必要な両物理量ｕk，ｕiの修正量Δｕk，Δｕiを算出し、ステップＳＰ７において、算出された修正量Δｕk，Δｕiに基づいて各物理量ｕk，ｕiを修正する。そして、ステップＳＰ８において推定誤差Ｄを上記修正量および対応する比例定数に基づいて修正し、ステップＳＰ９において全ての物理量に対する処理が行なわれたか否かを判別し、処理が行なわれていない物理量が存在すれば、ステップＳＰ１０において、他の２つの物理量を選択し、再びステップＳＰ２の処理を行なう。逆に、ステップＳＰ９において全ての物理量に対する処理が行なわれたと判別された場合には、ステップＳＰ１１において評価関数Ｅの値が所定の閾値よりも小さいか否かを判別し、評価関数Ｅの値が所定の閾値よりも小さければステップＳＰ１２において、最終的に得られた物理量を解析結果として採用し、そのまま一連の処理を終了する。逆に、評価関数Ｅの値が所定の閾値以上であれば、再びステップＳＰ１の処理を行なう。
【０１４０】
さらに詳細に説明する。
実施例５においては、修正すべき物理量が１つであるから、評価関数Ｅは放物線である。しかし、この実施例においては、修正すべき物理量が２つであるから、評価関数Ｅは図２０に示すように楕円になる。そして、ステップＳＰ２の処理を行なうことにより、図２０中にｐ1で示すように楕円状の点が得られ、次いで、ステップＳＰ３，ＳＰ４，ＳＰ５の処理を順次行なうことにより、楕円の中心により近い楕円状の点、ｐ2，ｐ3，ｐ4が得られる。これらの点ｐ1，ｐ2、ｐ3，ｐ4の関係は、図２１に示すとおりであり、ｐ1とｐ2との間隔をΔｕi´´、ｐ3とｐ4との間隔をΔｕi´、ｐ2とｐ3との間隔をΔｕk´、ｐ1とｐ3との間隔をｒ1、ｐ3と楕円中心との間隔をｒ2、ｐ2と楕円中心とのｕｉ軸方向の間隔をΔｕi、ｐ4と楕円中心とのｕk軸方向の間隔をΔｕkとすれば、数２４から数２６の比例関係式が得られる。
【０１４１】
【数２４】
【数２５】
【数２６】
したがって、数２４から数２７が得られる。
【０１４５】
【数２７】
そして、数２７を数２５，数２６に代入することにより、数２８のとおり、Δｕk，Δｕiを得ることができる。
【０１４７】
【数２８】
また、上記ステップＳＰ８における推定誤差Ｄの修正は、Ｄj←Ｄj−αijΔｕk−αkjΔｕiの演算により達成される。
以上から明らかなように、任意の２つの物理量を選択し、各物理量について交互に２回ずつ修正処理を行なった後は、両物理量を楕円中心に遷移させるために必要な修正量を簡単に算出することができる。
【０１４９】
【実施例７】
図２２はこの発明の物理量解析装置の他の実施例を示すブロック図であり、解析対象となる物理源の数と等しい数の格子ユニット１１１，１１２，・・・，１１ｍと、成立条件の数と等しい数の誤差ユニット１２１，１２２，・・・，１２ｎと、格子ユニット、誤差ユニットを制御する制御回路１３０とを有している。
【０１５０】
図２３は格子ユニット１１ｉの構成を詳細に示すブロック図であり、誤差ユニット１２ｊから供給される誤差値ＤJを入力として該当する物理量ｕiの修正量Δｕiを算出する修正量算出部１１ｉａと、物理量推定値ｕiを保持する推定値保持部１１ｉｂと、修正量算出部１１ａからの修正量Δｕiの出力を許容する状態と禁止する状態とを選択する選択部１１ｉｃと、選択部１１ｉｃにより出力が許容された修正量Δｕiと推定値保持部１１ｉｂに保持されている推定値ｕiとを加算して再び推定値保持部１１ｉｂにの供給する推定値修正部１１ｉｄとを有している。尚、選択部１１ｉｃにより出力が許容された修正量Δｕiはそのまま誤差ユニットに対する出力として使用される。また、修正量算出部１１ｉａは、誤差信号ＤJに対応して定数δij（＝βiαij）が設定されてあり、数２９の演算を行なって修正量Δｕiを算出する。
【０１５１】
【数２９】
図２４は誤差ユニット１２ｊの構成を詳細に示すブロック図であり、格子ユニット１１ｉから出力される修正量Δｕiを入力として該当する誤差Ｄjの修正量αijΔｕiを出力する誤差修正量演算部１２ｊａと、誤差値Ｄjを保持する誤差値保持部１２ｊｂと、誤差修正量演算部１２ｊａから出力される誤差修正量αijΔｕiと誤差値保持部１２ｊｂに保持されている誤差値Ｄjとを減算して再び誤差値保持部１２ｊｂに供給する誤差値修正部１２ｊｃとを有している。尚、誤差値保持部１２ｊｂに保持されている誤差値Ｄjはそのまま格子ユニットに対する出力として使用される。上記誤差修正量演算部１２ｊａは、修正量Δｕiに対応して定数αijが設定されてあり、αijΔｕiの演算を行なって誤差修正量を算出する。
【０１５３】
上記制御回路１３０は、選択部１１ｉｃを選択的に動作させる選択制御信号を出力するとともに、推定値保持部１１ｉｂに対して推定値修正部１１ｉｄからの出力値を保持すべきことを示す指示信号を出力し、また、誤差値保持部１２ｊｂに対して誤差値修正部１２ｊｃからの出力値を保持すべきことを示す指示信号を出力するものである。
【０１５４】
したがって、この実施例の場合には、実施例５と同様の作用を達成することができる。そして、上記の説明から明らかなように、メモリとして推定値保持部１１ｉｂと誤差値保持部１２ｊｂとが必要なだけであるから、メモリ容量を著しく低減することができる。即ち、メモリ容量の大幅な低減および物理源解析所要時間の大幅な短縮を両立させることができる。
【０１５５】
【実施例８】
図２５は格子ユニット２１ｉの他の構成例を示すブロック図であり、図２６は楕円中心に対する修正量を算出するための修正量算出ユニット２２の構成を示すブロック図である。尚、誤差修正量演算ユニットの構成は図２３と同一であるから図示および説明を省略する。
【０１５６】
上記格子ユニット２１ｉは、図２５に示すように、図２３の構成とほぼ同一であり、異なる点は、選択部１１ｉｃに代えて、修正量算出ユニット２２から出力される２つの修正量Δｕi，Δｕkおよび修正量算出部１１ｉａから出力される修正量Δｕiの何れかを選択するか、または何れも選択しないことを選択する選択部２１ｉｃを採用した点、および選択部２１ｉｃにより選択された修正量をそのまま修正量算出ユニット２２に対する出力としても採用した点のみである。
【０１５７】
修正量算出ユニット２２は、図２６に示すように、全ての格子ユニット２１ｉから出力される修正量Δｕiを入力とするマルチプレクサ２２ａと、マルチプレクサ２２ａにより選択された２つの修正量Δｕi，Δｕkをそれぞれ保持する１対の修正量保持部２２ｂ，２２ｃと、両修正量保持部２２ｂ，２２ｃにより保持されている修正量が修正方向を規定する基準値として供給され、しかも何れか一方の修正量を入力として数３０の演算を行なって２つの修正量Δｕi，Δｕkを出力する修正量算出部２２ｄとを有している。
【０１５８】
尚、この実施例においては、特には図示していないが、制御回路１３０によりマルチプレクサ２２ａに対して選択制御信号を供給するとともに、両修正量保持部２２ｂ，２２ｃに対して選択的に、マルチプレクサ２２ａから出力される修正値を保持すべきことを示す指示信号を供給するよう構成されている。したがって、マルチプレクサ２２ａにより１つの修正量を選択して一方の修正量保持部に保持させた後、マルチプレクサ２２ａにより他の修正量を選択して他方の修正量保持部に保持させ、次いで、マルチプレクサ２２ａにより交互に上記２つの修正量Δｕi，Δｕkを選択して、実施例６と同様の処理を行なって楕円中心に対応する２つの修正量Δｕi，Δｕkを算出することができる。尚、他の構成部分の作用は実施例７と同様であるから、詳細な説明を省略する。
【０１５９】
【発明の効果】
以上のように請求項１の発明は、磁場源の解析を高速に達成することができ、しかも磁場のｚ方向成分のみを計測することにより、各格子点におけるx方向電流成分、ｙ方向電流成分を解析することができるという特有の効果を奏する。
【図面の簡単な説明】
【図１】この発明の物理量解析方法の一実施例を説明するフローチャートである。
【図２】誤差関数の一例を示す図である。
【図３】分布電流を格子近似した磁場源モデルと磁場センサ配列面との関係を示す概略図である。
【図４】この発明の物理量解析装置の一実施例を示すブロック図である。
【図５】格子ユニットおよび修正ユニットの構成を詳細に示すブロック図である。
【図６】観測ユニットおよび物理場演算ユニットの構成を詳細に示すブロック図である。
【図７】格子ユニットおよび修正ユニットの構成を詳細に示すブロック図である。
【図８】観測ユニットおよび物理場演算ユニットの構成を詳細に示すブロック図である。
【図９】格子ユニットおよび修正ユニットの構成を詳細に示すブロック図である。
【図１０】この発明の他の実施例としての適応ノイズキャンセル方法を説明するフローチャートである。
【図１１】適応ノイズキャンセル方法の具体的適用例を示す波形図である。
【図１２】この発明のさらに他の実施例としての音響探査方法を説明するフローチャートである。
【図１３】音響探査の具体例を示す波形図である。
【図１４】この発明の物理量解析方法のさらに他の実施例を説明するフローチャートである。
【図１５】多数の電極により収集した体表面電位分布データに基づいて体内電気活動状態を可視化するための装置を示す概略図である。
【図１６】図１５の装置を具体化するための処理の一部として電位パターンを得るための物理量解析方法の具体例を説明するフローチャートである。
【図１７】３次元格子状抵抗ネットモデルの一部を拡大して示す図である。
【図１８】仮定された電流源と解析結果とを示す図である。
【図１９】この発明の物理量解析方法のさらに他の実施例を説明するフローチャートである。
【図２０】修正すべき物理量を２つ選択した場合における誤差評価関数の一例を示す図である。
【図２１】図１９のフローチャートに基づく修正量算出の原理を説明する概略図である。
【図２２】この発明の物理量解析装置の他の実施例を示すブロック図である。
【図２３】格子ユニットの構成を詳細に示すブロック図である。
【図２４】誤差修正量演算ユニットの構成を詳細に示すブロック図である。
【図２５】格子ユニットの他の構成例を詳細に示すブロック図である。
【図２６】楕円中心に対する修正量を算出するための修正量算出ユニットの構成を示すブロック図である。
【図２７】従来の適応ノイズキャンセラを示すブロック図である。
【図２８】従来の物理量解析装置を示すブロック図である。
【符号の説明】
１ｉａ 修正ユニット ２ｉ 観測ユニット
３１ｉ 修正量基本値演算セル ３２ｉ 修正値演算セル
３３ｉ 推定値保持セル ４２ｉ 誤差演算セル
１１ｉａ 修正量算出部 １１ｉｂ 推定値保持部
１１ｉｃ 選択部 １１ｉｄ 推定値修正部
１２ｊａ 誤差修正量演算部 １２ｊｂ 誤差値保持部
１２ｊｃ 誤差値修正部
Claims (1)
- 測定対象物が発生する磁場を複数の磁場センサで計測し、前記測定対象物内部の磁場源における物理量を解析する装置であって、
（１） 個々の磁場源ｉに対応する単位量の磁場源が各所定箇所に発生する磁場を示す、数１で与えられる定数αijを予め定めておく手段と、
（２） 物の内部に設定された格子点ｋ（ｋ＝１，２・・・，ｐ）に位置する各磁場源の、数２で示すx方向電流成分およびｙ方向電流成分からなる物理量ｕiを仮に定めておく手段と、
（３） 数３の演算を行なって観測点ｊに発生すると予測される磁場Ｏjを算出する手段と、
（４） 観測点ｊにおける磁場の実測値Ｓjを測定する手段と、
（５） 観測点ｊにおける磁場の実測値Ｓjと上記磁場Ｏjとの差Ｓj−Ｏjを算出する手段と、
（６） 数４の演算を行なって仮の解としての物理量ｕiを得る手段と、
（７） （３）から（６）までの各手段による処理を全ての磁場源ｉについて行なわせる手段と、
（８） 上記差Ｓj−Ｏjの二乗の総和が所定の閾値よりも小さくなるまで、（３）から（６）までの各手段による処理の反復処理を反復させる手段と、
（９） 最終的に得られた各磁場源の物理量ｕiを磁場源の物理量ｕiとして採用する手段とを含むことを特徴とする物理量解析装置。
Priority Applications (7)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP30557193A JP3718856B2 (ja) | 1993-06-03 | 1993-12-06 | 物理量解析装置 |
| US08/381,823 US5646868A (en) | 1993-06-03 | 1994-06-03 | Physical quantity analyzing method and apparatus thereof |
| PCT/JP1994/000901 WO1994029806A1 (fr) | 1993-06-03 | 1994-06-03 | Procede et appareil d'analyse d'une quantite physique |
| AU68565/94A AU683063B2 (en) | 1993-06-03 | 1994-06-03 | Method and apparatus for analysis of physical quantity |
| CA002141876A CA2141876A1 (en) | 1993-06-03 | 1994-06-03 | Method and apparatus for analysis of physical quantity |
| EP94917152A EP0654744A4 (en) | 1993-06-03 | 1994-06-03 | METHOD AND APPARATUS FOR ANALYSIS OF A PHYSICAL QUANTITY. |
| KR1019950700395A KR950702724A (ko) | 1993-06-03 | 1995-02-03 | 물리량 해석방법 및 그 장치 |
Applications Claiming Priority (3)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP13370793 | 1993-06-03 | ||
| JP5-133707 | 1993-06-03 | ||
| JP30557193A JP3718856B2 (ja) | 1993-06-03 | 1993-12-06 | 物理量解析装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JPH0749857A JPH0749857A (ja) | 1995-02-21 |
| JP3718856B2 true JP3718856B2 (ja) | 2005-11-24 |
Family
ID=26467984
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP30557193A Expired - Fee Related JP3718856B2 (ja) | 1993-06-03 | 1993-12-06 | 物理量解析装置 |
Country Status (7)
| Country | Link |
|---|---|
| US (1) | US5646868A (ja) |
| EP (1) | EP0654744A4 (ja) |
| JP (1) | JP3718856B2 (ja) |
| KR (1) | KR950702724A (ja) |
| AU (1) | AU683063B2 (ja) |
| CA (1) | CA2141876A1 (ja) |
| WO (1) | WO1994029806A1 (ja) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2007236737A (ja) * | 2006-03-10 | 2007-09-20 | Kanazawa Inst Of Technology | 電流素片推定方法および装置 |
Families Citing this family (14)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6064810A (en) * | 1996-09-27 | 2000-05-16 | Southern Methodist University | System and method for predicting the behavior of a component |
| JP5083194B2 (ja) * | 2008-12-18 | 2012-11-28 | 株式会社デンソーウェーブ | ロボットのキャリブレーション方法及びロボットの制御装置 |
| US9063882B1 (en) * | 2010-09-09 | 2015-06-23 | Sas Ip, Inc. | Matrix preconditioners for simulations of physical fields |
| US8738554B2 (en) | 2011-09-16 | 2014-05-27 | International Business Machines Corporation | Event-driven universal neural network circuit |
| US8874498B2 (en) | 2011-09-16 | 2014-10-28 | International Business Machines Corporation | Unsupervised, supervised, and reinforced learning via spiking computation |
| US8799199B2 (en) | 2011-12-14 | 2014-08-05 | International Business Machines Corporation | Universal, online learning in multi-modal perception-action semilattices |
| US8626684B2 (en) | 2011-12-14 | 2014-01-07 | International Business Machines Corporation | Multi-modal neural network for universal, online learning |
| US8826216B2 (en) | 2012-06-18 | 2014-09-02 | International Business Machines Corporation | Token-based current control to mitigate current delivery limitations in integrated circuits |
| US8863068B2 (en) | 2012-06-18 | 2014-10-14 | International Business Machines Corporation | Current-aware floorplanning to overcome current delivery limitations in integrated circuits |
| US8914764B2 (en) | 2012-06-18 | 2014-12-16 | International Business Machines Corporation | Adaptive workload based optimizations coupled with a heterogeneous current-aware baseline design to mitigate current delivery limitations in integrated circuits |
| US8826203B2 (en) * | 2012-06-18 | 2014-09-02 | International Business Machines Corporation | Automating current-aware integrated circuit and package design and optimization |
| JP2017191039A (ja) * | 2016-04-14 | 2017-10-19 | セイコーエプソン株式会社 | 磁場計測装置及び磁場計測装置の校正方法 |
| JP2019124595A (ja) * | 2018-01-17 | 2019-07-25 | 公立大学法人首都大学東京 | 物体内への物流入量算出方法 |
| JP6947194B2 (ja) * | 2019-02-13 | 2021-10-13 | Tdk株式会社 | 信号処理回路および磁気センサシステム |
-
1993
- 1993-12-06 JP JP30557193A patent/JP3718856B2/ja not_active Expired - Fee Related
-
1994
- 1994-06-03 WO PCT/JP1994/000901 patent/WO1994029806A1/ja not_active Application Discontinuation
- 1994-06-03 CA CA002141876A patent/CA2141876A1/en not_active Abandoned
- 1994-06-03 AU AU68565/94A patent/AU683063B2/en not_active Expired
- 1994-06-03 EP EP94917152A patent/EP0654744A4/en not_active Withdrawn
- 1994-06-03 US US08/381,823 patent/US5646868A/en not_active Expired - Lifetime
-
1995
- 1995-02-03 KR KR1019950700395A patent/KR950702724A/ko active IP Right Grant
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2007236737A (ja) * | 2006-03-10 | 2007-09-20 | Kanazawa Inst Of Technology | 電流素片推定方法および装置 |
Also Published As
| Publication number | Publication date |
|---|---|
| EP0654744A1 (en) | 1995-05-24 |
| EP0654744A4 (en) | 1997-07-23 |
| AU6856594A (en) | 1995-01-03 |
| CA2141876A1 (en) | 1994-12-22 |
| KR950702724A (ko) | 1995-07-29 |
| WO1994029806A1 (fr) | 1994-12-22 |
| JPH0749857A (ja) | 1995-02-21 |
| US5646868A (en) | 1997-07-08 |
| AU683063B2 (en) | 1997-10-30 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP3718856B2 (ja) | 物理量解析装置 | |
| US8369941B2 (en) | High definition impedance imaging | |
| CN108027409A (zh) | 时域mri | |
| US20050251062A1 (en) | Iterative approach for applying multiple currents to a body using voltage sources in electrical impedance tomography | |
| US9207197B2 (en) | Coil for magnetic induction to tomography imaging | |
| CN110990757A (zh) | 利用无相位数据解决高度非线性电磁逆散射问题的方法 | |
| JPH0542116A (ja) | 生体表面電位分布から生体内電流分布を推定する信号処理方法 | |
| Fan et al. | Maximum entropy regularization method for electrical impedance tomography combined with a normalized sensitivity map | |
| Lehtikangas et al. | Reconstruction of velocity fields in electromagnetic flow tomography | |
| Fischer et al. | Acoustic tomography of the atmosphere: Comparison of different reconstruction algorithms | |
| EP0477434A1 (en) | Analysis of biological signals using data from arrays of sensors | |
| Hakula et al. | On the hp-adaptive solution of complete electrode model forward problems of electrical impedance tomography | |
| ES2400025T3 (es) | Sensor de monedas | |
| US6697661B2 (en) | Apparatus, system and method for calibrating magnetic resonance receiver coils | |
| Abdallh et al. | Optimal needle placement for the accurate magnetic material quantification based on uncertainty analysis in the inverse approach | |
| US20160198975A1 (en) | System and method for analysing data from a microwave inverse scattering apparatus | |
| Zarinabad et al. | Modelling parameter role on accuracy of cardiac perfusion quantification | |
| Cutrupi et al. | An approach to the electrical resistance tomography based on meshless methods | |
| JPH06348870A (ja) | 物理量解析方法およびその装置 | |
| Marshall | Detecting and measuring sources at the noise limit | |
| JPH05220124A (ja) | 生体磁気計測装置 | |
| Imperiale et al. | Some fast calculations simulating measurements from single-photon emission computed tomography (SPECT) imaging | |
| JP2949983B2 (ja) | 物理源解析方法およびその装置 | |
| Janagal | Modeling Micro-gradient Magnetic Field Distribution with FEM | |
| Bondariev et al. | Imitational modeling of synchronous measuring detector with increased noise immunity and precision |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20040406 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20040607 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20050222 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20050816 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20050829 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090916 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100916 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100916 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110916 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110916 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120916 Year of fee payment: 7 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130916 Year of fee payment: 8 |
|
| LAPS | Cancellation because of no payment of annual fees |