JP5386472B2 - 拡散物質の発生源推定装置および発生源推定方法 - Google Patents

Description

本発明は、拡散物質の発生源を推定する拡散物質の発生源推定装置および発生源推定方法に関するものである。

プラント施設（火力発電所、ゴミ焼却施設、化学プラント等）での事故等による汚染物質の放出や、テロ等による毒性ガス等の放出に対して、その発生源を同定して即座に対処する必要性から、事故または事件の現場情報（濃度計測値等）から発生源情報（放出地点の位置および発生量）を推定する拡散物質の発生源推定装置および発生源推定方法について各種技法が提案されている。

他方、多様な分野で応答や結果などの出力から原因や入力を推定する逆問題(inverse problem)の重要性が認識されてきており、盛んに適用されるようになってきている。ここで、逆問題の取扱いを逆問題解析あるいは逆解析(inverse analysis)という。例えば、非特許文献１の「時間逆解析での数値安定性確保についての研究」では、拡散物質の発生源を推定する手法にこの逆解析を用いており、拡散物質の輸送方程式を時間方向に逆に解析して、放出地点の位置、時刻および放出量を同定する手法が提案されている。

安部諭,加藤信介，"時間逆解析での数値安定性確保についての研究"，第２３回数値流体力学シンポジウム，Ｇ２−５（２００９年）

一般的に、偏微分方程式等に関する逆問題では、偏微分方程式等を差分化して問題を解くことが多く、安定性が得られず非適切な問題となることが多いことが指摘されている。非特許文献１においても、拡散物質の輸送方程式を差分化しているため、数値不安定性の問題が常に伴い、フィルタ操作など何かしらの工夫を要するという事情があった。また、非特許文献１では、観測地点の濃度をそのまま輸送方程式に代入していると考えられ、解が１つに定まらない可能性がある。

本発明は、上記従来の事情に鑑みてなされたものであって、逆解析を用いて発生源を推定する際の数値不安定性の問題を解消した拡散物質の発生源推定装置および発生源推定方法を提供することを目的としている。

上記課題を解決するため、本発明は以下の手段を採用する。
本発明に係る拡散物質の発生源推定装置は、複数個の観測器からの情報に基づき、ガス発生源情報を推定する拡散物質の発生源推定装置であって、前記観測器から位置情報、計測した濃度情報および計測時刻を入手する観測情報入手手段と、各観測器の位置をそれぞれ仮想放出点とし、各仮想放出点の放出強度をそれぞれ対応する観測器で得た濃度情報に逆比例した値とする逆問題設定手段と、各仮想放出点からの放出について、粒子法を用いた拡散の数値解析で時間を前記計測時刻から逆に進行させて逆解析を行う逆解析手段と、逆に進行させたそれぞれの時刻において、各仮想放出点からの粒子濃度分布を評価し、粒子濃度が略等しい領域の内、粒子濃度の分散が最小となる位置および時刻をそれぞれ放出地点および放出時刻と推定する推定手段と、を備えることを特徴とする。

本発明によれば、各観測器の位置をそれぞれ仮想放出点とし、各仮想放出点の放出強度をそれぞれ対応する観測器で得た濃度情報に逆比例した値として逆問題を設定し、粒子法を用いた拡散の数値解析で時間を計測時刻から逆に進行させる逆解析に基づく評価によって発生源を推定することとしたので、差分化による数値不安定性の問題を解消した拡散物質の発生源推定装置を提供することができる。

また、本発明は、上記記載の拡散物質の発生源推定装置において、記推定手段は、粒子濃度の分散が最小となる位置における粒子濃度の逆数値を前記放出地点の放出強度と推定することを特徴とする。

本発明によれば、放出地点、放出時刻および放出強度をそれぞれ推定することができる。

また、本発明は、上記記載の拡散物質の発生源推定装置において、前記逆解析手段は、最も遅い計測時刻との相対時刻で各仮想放出点から順次放出させていくのが望ましい。

また、本発明は、上記記載の拡散物質の発生源推定装置において、相互の離間距離が一定の評価格子上で各格子線が交差する位置を評価点として設定する評価点設定手段、を備え、前記推定手段は、逆に進行させたそれぞれの時刻において、設定された評価点について各仮想放出点からの粒子濃度分布を評価し、粒子濃度の分散が最小となる位置および時刻をそれぞれ放出地点および放出時刻と推定することを特徴とする。

本発明によれば、評価格子上での評価点の設定により、推定手段における粒子濃度分布評価の処理が簡便となり、装置設計をより容易なものとすることができる。

また、本発明は、上記記載の拡散物質の発生源推定装置において、前記評価点設定手段は、前記推定手段で推定された放出地点を含み、且つ相互の離間距離がより短い評価格子上で各格子線が交差する位置を評価点として再設定することを特徴とする。

本発明によれば、格子線の相互の離間距離が相対的に長い仮想格子から相対的に短い仮想格子へと絞り込みながら評価点を再設定して発生源を推定するようにすれば、離間距離を最短とした１面の仮想格子上で評価点を設定して発生源を推定する場合と比べて、１面の仮想格子上での評価点数を格段に小さくすることができ、全体的な処理の計算量を抑制してより速い発生源の推定を行うことができる。

また、本発明に係る拡散物質の発生源推定方法は、複数個の観測器からの情報に基づき、ガス発生源情報を推定する拡散物質の発生源推定方法であって、前記観測器から位置情報、計測した濃度情報および計測時刻を入手する観測情報入手ステップと、各観測器の位置をそれぞれ仮想放出点とし、各仮想放出点の放出強度をそれぞれ対応する観測器で得た濃度情報に逆比例した値とする逆問題設定ステップと、各仮想放出点からの放出について、粒子法を用いた拡散の数値解析で時間を前記計測時刻から逆に進行させて逆解析を行う逆解析ステップと、逆に進行させたそれぞれの時刻において、各仮想放出点からの粒子濃度分布を評価し、粒子濃度が略等しい領域の内、粒子濃度の分散が最小となる位置および時刻をそれぞれ放出地点および放出時刻と推定する推定ステップと、備えることを特徴とする。

本発明によれば、各観測器の位置をそれぞれ仮想放出点とし、各仮想放出点の放出強度をそれぞれ対応する観測器で得た濃度情報に逆比例した値として逆問題を設定し、粒子法を用いた拡散の数値解析で時間を計測時刻から逆に進行させる逆解析に基づく評価によって発生源を推定することとしたので、差分化による数値不安定性の問題を解消した拡散物質の発生源推定方法を提供することができる。

また、本発明は、上記記載の拡散物質の発生源推定方法において、記推定ステップは、粒子濃度の分散が最小となる位置における粒子濃度の逆数値を前記放出地点の放出強度と推定することを特徴とする。

本発明によれば、放出地点、放出時刻および放出強度をそれぞれ推定することができる。

また、本発明は、上記記載の拡散物質の発生源推定方法において、前記逆解析ステップは、最も遅い計測時刻との相対時刻で各仮想放出点から順次放出させていくのが望ましい。

また、本発明は、上記記載の拡散物質の発生源推定方法において、相互の離間距離が一定の評価格子上で各格子線が交差する位置を評価点として設定する評価点設定ステップ、を備え、前記推定ステップは、逆に進行させたそれぞれの時刻において、設定された評価点について各仮想放出点からの粒子濃度分布を評価し、粒子濃度の分散が最小となる位置および時刻をそれぞれ放出地点および放出時刻と推定することを特徴とする。

本発明によれば、評価格子上での評価点の設定により、推定手段における粒子濃度分布評価の処理が簡便となり、プログラム設計をより容易なものとすることができる。

また、本発明は、上記記載の拡散物質の発生源推定方法において、前記評価点設定ステップは、前記推定手段で推定された放出地点を含み、且つ相互の離間距離がより短い評価格子上で各格子線が交差する位置を評価点として再設定することを特徴とする。

本発明によれば、格子線の相互の離間距離が相対的に長い仮想格子から相対的に短い仮想格子へと絞り込みながら評価点を再設定して発生源を推定するようにすれば、離間距離を最短とした１面の仮想格子上で評価点を設定して発生源を推定する場合と比べて、１面の仮想格子上での評価点数を格段に小さくすることができ、全体的な処理の計算量を抑制してより速い発生源の推定を行うことができる。

本発明によれば、差分化による数値不安定性の問題を解消した拡散物質の発生源推定装置および発生源推定方法を提供することができるという効果を奏する。

本発明の第１実施形態に係る拡散物質の発生源推定装置の構成図である。 拡散現象における逆問題および逆解析を説明する説明図である。 第１実施形態の拡散物質の発生源推定方法を説明するフローチャートである。 仮想放出点からの拡散を例示する説明図である。 仮想放出点からの放出タイミングおよび評価タイミングを例示するタイムチャートである。 第２実施形態の拡散物質の発生源推定方法を説明するフローチャートである。 仮想格子に基づく評価点を用いた放出候補地点の推定を例示する説明図である。 仮想格子のネスティングを説明する説明図である。

以下、本発明の拡散物質の発生源推定装置および発生源推定方法の実施形態について、第１実施形態、第２実施形態の順に図面を参照して詳細に説明する。

〔第１実施形態〕
図１は本発明の第１実施形態に係る拡散物質の発生源推定装置の構成図である。
同図において、本実施形態の拡散物質の発生源推定装置３は、通信インタフェース１１，入力部１３、発生源推定処理部１６、記憶部１７および出力部１９を備えて構成されている。すなわち、発生源推定装置３はいわゆるコンピュータシステムの構成であり、発生源推定処理部１６はＭＰＵ（マイクロプロセッサ）やＤＳＰ（ディジタルシグナルプロセッサ）等のプロセッサで具現される。また、記憶部１７はＲＡＭ（Random Access Memory）、ＲＯＭ（Read Only Memory）、ＨＤＤ（Hard Disk Drive）等の記憶装置で、入力部１３はキーボード、マウス等の入力装置で、出力部１９はディスプレイ、プリンタ等の出力装置でそれぞれ具現される。

また、図１において、発生源推定装置３は、ｎ個の観測器５−１〜５−ｎ（ｎは正整数）からの情報、即ち該観測器の位置情報並びに該観測器が計測した濃度情報および計測時刻情報を、通信インタフェース１１を介して入手する構成である。観測器５−１〜５−ｎは、少なくともその設置場所における大気中の所望ガスの濃度を計測する機能を備え、定期的に該観測器の位置情報並びに該観測器が計測した濃度情報および計測時刻情報を発信する機能も備える。なお、観測器５−１〜５−ｎからリアルタイムに情報入手が可能な場合には、観測器の位置情報および濃度情報のみを受信して、発生源推定装置３側での受信タイミングで計測時刻を代用するようにしても良い。

また、観測器５−１〜５−ｎは固定設置または可動設置の何れでも良く、可動設置の場合には、ＧＰＳ機能を備えて、常時、該観測器の位置情報を取得可能な構成とすれば良い。またさらに、固定設置の場合には、観測器の位置情報を固有の識別コード（機器番号等）で代用することも可能である。この場合、発生源推定装置３側で、観測器の識別コードに対応したテーブル等に基づき観測器の位置情報を導出することになる。

また、観測地点から情報を入手する方法は、図１に例示した構成に限定されない。例えば、民間または公共機関等によるガス濃度観測システムが存在して、各観測地点における濃度データを逐次蓄積していくデータベース機能を備えたものである場合には、通信インタフェース１１を介しインターネット等のネットワーク上に設置されたデータベースにアクセスして、各観測地点の位置データ並びに該観測器が計測した濃度データおよび観測時刻データを入手するようにしても良い。なお、ガス放出に対してその発生源を即座に同定する観点からは、本実施形態の構成（図１）が望ましい。

また、発生源推定装置３の発生源推定処理部１６には、観測情報入手部２１、逆問題設定部２３、評価点設定部２５、逆解析部２６および推定部２７を備えている。これら各構成要素は、ＭＰＵやＤＳＰ等のプロセッサ上で実行されるプログラムの機能的まとまりとして具現されるものである。

ここで、発生源推定処理部１６の各構成要素の具体的な機能を説明する前に、本発明における拡散物質の発生源を推定する手法の基本的な考え方について、図２を参照して説明する。図２は拡散現象における逆問題および逆解析を説明する説明図である。

まず、図２（ａ）に示すように、放出地点Ｐｏからの拡散物質が放出され、観測地点Ｐｖで濃度を観測する場合について考える。なお、放出地点近傍でｘ方向に一様な風が吹いていると仮定し、風に直角な方向をｙ方向とする。このとき、放出位置を（ｘ_０，ｙ_０，ｔ_０）、放出強度をｑ_０とすると、任意の観測位置（ｘ_ｉ，ｙ_ｉ，ｔ_ｉ）における濃度は次式で表される。
［数１］
Ｄ_ｉ（ｘ_ｉ，ｙ_ｉ，ｔ_ｉ）＝ｑ_０Ｄ_ｉ０（ｘ_ｉ−ｘ_０，ｙ_ｉ−ｙ_０，ｔ_ｉ−ｔ_０） （１）
ここに、Ｄ_ｉ０（ｘ_ｉ−ｘ_０，ｙ_ｉ−ｙ_０，ｔ_ｉ−ｔ_０）は観測位置（ｘ_ｉ，ｙ_ｉ，ｔ_ｉ）と放出位置（ｘ_０，ｙ_０，ｔ_０）との相対位置および相対時刻によって決まる影響関数である。

また、複数の観測地点（即ち、ｎ箇所（ｎは正整数）の評価位置）における濃度（Ｄ_ｉ（ｘ_ｉ，ｙ_ｉ，ｔ_ｉ）；ｉ＝１，２，…，ｎ）が計測されていれば、次式が成り立つ。

一方、観測地点Ｐｖで観測した濃度から放出地点Ｐｏを推定することを逆問題として設定し、図２（ｂ）に示すように、観測地点Ｐｖを仮想放出点として逆拡散させる場合を考える。観測位置（ｘ_ｉ，ｙ_ｉ，ｔ_ｉ）から強度ｑ_ｉを仮想的に放出して時間を逆進させると、放出位置（ｘ_０，ｙ_０，ｔ_０）における濃度は式（１）と同様に次式で表される。
［数３］
Ｄ_０（ｘ_０，ｙ_０，ｔ_０）＝ｑ_ｉＤ_０ｉ（ｘ_０−ｘ_ｉ，ｙ_０−ｙ_ｉ，ｔ_０−ｔ_ｉ） （３）

ここで、拡散現象の性質から、Ｄ_ｉ０＝Ｄ_０ｉが成立し、次式を得る。

したがって、

とすれば、即ち、観測濃度に逆比例する強度を観測地点（仮想放出点）からの放出量にすれば、式（４）は次式に変形できる。
［数６］
Ｄ_０（ｘ_０，ｙ_０，ｔ_０）＝α／ｑ_０ （６）

式（６）から、「どの観測地点からの仮想的な放出に対しても、該仮想放出に対する放出地点Ｐｏでの濃度は同じであり、且つ放出強度はその濃度の逆数である」ことが分かる。本発明はこの性質を利用したもので、各観測地点をそれぞれ仮想放出点とし、各仮想放出点の放出強度をそれぞれ観測濃度に逆比例した値として逆問題を設定し、各仮想放出点からの放出について、粒子法を用いた拡散の数値解析で時間を計測時刻から逆に進行させて逆解析を行い、逆に進行させたそれぞれの時刻において、各仮想放出点からの粒子濃度分布を評価し、粒子濃度が略等しい領域の内、最も収束した位置および時刻をそれぞれ放出地点および放出時刻と推定する拡散物質の発生源推定手法を提案する。

例えば、図２（ｃ）に示すように、放出地点Ｐｏからの放出を放出地点Ｐｏによる拡散領域内に設置された２箇所の観測地点Ｐｖ１およびＰｖ２で観測するとした場合、各観測地点Ｐｖ１およびＰｖ２をそれぞれ仮想放出点とし、それら仮想放出点の放出強度をそれぞれの観測濃度に逆比例した値として仮想的な逆拡散を行えば、仮想放出点Ｐｖ１およびＰｖ２によるそれぞれの逆拡散領域が重複する領域が大まかな推定範囲とすることができ、また、仮想放出点Ｐｖ１およびＰｖ２からの粒子濃度分布について、粒子濃度が略等しい領域の内、最も収束した位置を放出地点Ｐｏと推定できる。

次に、以上説明した基本的な理論を踏まえて、発生源推定処理部１６の各構成要素（即ち、観測情報入手部２１、逆問題設定部２３、評価点設定部２５、逆解析部２６および推定部２７）の具体的な機能について説明する。

まず、観測情報入手部２１は、観測器５−１〜５−ｎからの情報、即ち該観測器の位置情報、該観測器が計測した濃度情報および計測時刻情報を、通信インタフェース１１を介して入手する。入手した各情報は関連づけて記憶部１７の所定領域に保存しておくのが望ましい。

また、逆問題設定部２３は、各観測器５−ｉ（ｉ＝１〜ｎ）の位置をそれぞれ仮想放出点とし、各仮想放出点の放出強度をそれぞれ対応する観測器５−ｉで得た濃度情報（観測濃度）に逆比例した値として逆問題を設定する。ここで、式（５）における係数αは任意に設定する。なお、後述のように粒子法による逆解析を行うので、具体的には、各仮想放出点の放出粒子数を全て同数とし、粒子強度を観測濃度に逆比例させる。なお、粒子濃度＝粒子数×粒子強度である。ここで、放出粒子数を同数とするのは、各仮想放出点からの粒子数依存度を減らすためである。

また、評価点設定部２５は、相互の離間距離が一定の仮想格子を想定して、該仮想格子の領域上で各格子線が交差する位置を評価点として設定する。例えば、ｒ×ｒの仮想格子上でｒ^２箇所の評価点Ｐｏｅｊ（ｊ＝１〜ｒ^２）が設定されるが、この評価点数が多ければ多いほど放出地点の推定精度は向上し、そのためには仮想格子の離間距離をより小さくして交差点数を増やせば良いが、その分だけ計算量は増大することになる。したがって、要求される処理時間（観測地点から情報を入手してから放出地点を推定するまでの時間）に応じて、予め発生源推定処理部１６を具現するプロセッサの処理性能等を勘案して評価点数を概算しておくのが望ましい。

また、逆解析部２６は、各仮想放出点からの放出について、粒子法を用いた拡散の数値解析で時間を計測時刻から逆に進行させて逆解析を行う。ここで、逆解析の計算コードは公知のものを使用すれば良く、例えば、コロラド州立大学と米国ATMET社で開発されたＨＹＰＡＣＴ（Hybrid Particle and Concentration Transport Model）を用いて良い。なお、ＨＹＰＡＣＴは、粒子法により拡散物質の大気中濃度を評価するものを使用する。ＨＹＰＡＣＴによれば、粒子群の重力による沈降（乾式沈着）、降雨による沈降（湿式沈着）、放射性核種の減衰、熱・運動エネルギーによる排気上昇をリアルタイムに計算し、空間濃度分布を求めることができる。

なお、粒子法における基本式は次式で表される。
［数７］
Ｘnew＝Ｘold＋ｕΔｔ
Ｙnew＝Ｙold＋ｖΔｔ （７）
Ｚnew＝Ｚold＋ｗΔｔ
ここで、Δｔは現時刻と次時刻間の時間刻み（正値）を、Ｘold，Ｙold，Ｚoldは現時刻の粒子座標を、Ｘnew，Ｙnew，Ｚnewは次時刻の粒子座標を、ｕ，ｖ，ｗは風速成分をそれぞれ表す。

また、時間を逆に進行させて逆解析を行う場合には、式（７）は次式のように変形される。
［数８］
Ｘold＝Ｘnew＋ｕΔｔ
Ｙold＝Ｙnew＋ｖΔｔ （８）
Ｚold＝Ｚnew＋ｗΔｔ
ここで、現時刻と次時刻間の時間刻みΔｔは負値をとるものとする。したがって、既存のＨＹＰＡＣＴの計算コードについて、時間刻みΔｔに負値を与えるように修正すれば良く、新たに計算コードを設計する必要はない。

また、逆解析部２６による逆解析では、観測器５−ｉ（ｉ＝１〜ｎ）それぞれの観測時刻が異なる場合には、最も遅い計測時刻の観測器に対応する仮想放出点から放出を開始し、その後、最も遅い計測時刻との相対時刻の短い順に、最も遅い計測時刻の観測器に対応する仮想放出点からの放出タイミングから該相対時刻遅れたタイミングで、他の仮想放出点から順次放出させていく。

また、推定部２７は、時間刻みΔｔずつ逆に進行させたそれぞれの時刻について、各仮想放出点からの粒子濃度分布を評価し、粒子濃度が略等しい領域の内、粒子濃度の分散が最小となる位置および時刻をそれぞれ放出地点および放出時刻と推定する。なお、粒子濃度分布の評価は、最も速い計測時刻の観測器に対応する仮想放出点からの放出タイミングより前の時刻まで逆に進行させた後に行うものとする。つまり、全観測器５−１〜５−ｎにそれぞれ対応する仮想放出点から放出されている状態に至ってから評価を開始する。

具体的に、粒子濃度分布の評価は以下のようにして行われる。すなわち、評価点設定部２５によりｒ×ｒの仮想格子上でｒ^２個の評価点Ｐｏｅｊ（ｊ＝１〜ｒ^２）が設定されているとき、各評価点Ｐｏｅｊに対する各仮想放出点からの粒子濃度分布の分散σ_ｊに基づき評価する。つまり、逆解析部２６により、ある時刻における各観測器５−ｉに対応する仮想放出点からの放出による粒子濃度Ｃ_ｉｊが略等しく、最も収束した位置を見いだすために、粒子濃度分布の分散σ_ｊが最小となる位置を求め、粒子濃度の分散が最小となる位置および時刻をそれぞれ放出地点および放出時刻と推定し、また、粒子濃度の分散が最小となる位置における粒子濃度の逆数値を放出地点の放出強度と推定する。

ここで、評価点Ｐｏｅｊに対する粒子濃度分布の分散σ_ｊおよび放出強度ｑ_ｊは次式で求められる。

次に、以上のような構成要素を備えた発生源推定装置３における拡散物質の発生源推定方法について、図３〜図５を参照して説明する。ここで、図３は本実施形態に係る拡散物質の発生源推定方法を説明するフローチャートであり、図４は仮想放出点からの拡散を例示する説明図であり、図５は仮想放出点からの放出タイミングおよび評価タイミングを例示するタイムチャートである。

まず、ステップＳ１０１では、観測情報入手部２１において、観測器５−ｉ（ｉ＝１〜ｎ）の位置情報、該観測器５−ｉが計測した濃度情報および計測時刻情報を、通信インタフェース１１を介して入手する。

次に、ステップＳ１０２では、逆問題設定部２３において、各観測器５−ｉ（ｉ＝１〜ｎ）の位置をそれぞれ仮想放出点Ｐｖｉ（ｉ＝１〜ｎ）とし、各仮想放出点の放出強度（粒子強度）をそれぞれ対応する観測器５−ｉで得た濃度情報に逆比例した値として逆問題を設定する。

次に、ステップＳ１０３では、評価点設定部２５により、相互の離間距離が一定の仮想格子を想定して、該仮想格子の領域上で各格子線が交差する位置を評価点Ｐｏｅｊ（ｊ＝１〜ｒ^２）として設定する。

次に、ステップＳ１０４では、逆解析部２６により、各仮想放出点からの放出について、粒子法を用いた拡散の数値解析で時間を計測時刻から逆に進行させて逆解析を行う。なお、観測器５−ｉ（ｉ＝１〜ｎ）それぞれの観測時刻が異なる場合には、最も遅い計測時刻の観測器に対応する仮想放出点から放出を開始し、その後、最も遅い計測時刻との相対時刻の短い順に、最も遅い計測時刻の観測器に対応する仮想放出点からの放出タイミングから該相対時刻遅れたタイミングで、他の仮想放出点から順次放出させていく。

例えば、図４に示すように、観測器５−１〜５−３で濃度が計測（計測時刻はそれぞれＴｖ１〜Ｔｖ３とする）されたとき、各観測器の位置をそれぞれ仮想放出点Ｐｖ１〜Ｐｖ３として、各仮想放出点から逆拡散させるとすれば、各仮想放出点からの放出タイミングは、図５に示すようになる。つまり、時刻Ｔｖ３に最も遅い計測時刻の観測器５−３に対応する仮想放出点Ｐｖ３から放出を開始し、その後、時刻Ｔｖ１に観測器５−１に対応する仮想放出点Ｐｖ１から、時刻Ｔｖ２に観測器５−２に対応する仮想放出点Ｐｖ２からそれぞれ放出させていく。

また、図５中の鉛直方向に引かれた破線は、それぞれ逆解析部２６による数値解析のタイミングを刻むが、逆解析部２６によって求められる仮想放出点Ｐｖ１〜Ｐｖ３からの放出による粒子濃度Ｃ_ｉｊ等のデータを推定部２７に渡すのは、最も速い計測時刻の観測器５−２に対応する仮想放出点Ｐｖ２からの放出タイミングより前まで逆に進行させた時点、即ち時刻Ｔｅ１のときからである。

次に、ステップＳ１０５では、推定部２７により、時間刻みΔｔずつ逆に進行させたそれぞれの時刻について、各仮想放出点Ｐｖｉ（ｉ＝１〜ｎ）からの粒子濃度分布を評価し、粒子濃度が略等しい領域の内、粒子濃度の分散が最小となる評価点Ｐｏｅｊの位置および時刻をそれぞれ放出地点および放出時刻と推定する。すなわち、各評価点Ｐｏｅｊに対する各仮想放出点からの粒子濃度分布の分散σ_ｊが最小となる評価点Ｐｏｅｊの位置を求め、粒子濃度の分散が最小となる評価点Ｐｏｅｊの位置および時刻をそれぞれ放出地点および放出時刻と推定し、また、粒子濃度の分散が最小となる評価点Ｐｏｅｊの位置における粒子濃度の逆数値を放出地点の放出強度と推定する。

以上説明したように、本実施形態の拡散物質の発生源推定装置および発生源推定方法では、観測器５−ｉ（ｉ＝１〜ｎ）の位置情報および計測した濃度情報を観測情報入手部２１によって入手し、逆問題設定部２３により、各観測器の位置をそれぞれ仮想放出点Ｐｖｉ（ｉ＝１〜ｎ）とし、各仮想放出点の放出強度をそれぞれ対応する観測器で得た濃度情報に逆比例した値として逆問題を設定し、逆解析部２６により、各仮想放出点からの放出について、粒子法を用いた拡散の数値解析で時間を計測時刻から逆に進行させて逆解析を行い、推定部２７により、逆に進行させたそれぞれの時刻において、各仮想放出点からの粒子濃度分布を評価し、粒子濃度が略等しい領域の内、粒子濃度の分散が最小となる位置および時刻をそれぞれ放出地点および放出時刻と推定する。

このように、各観測器５−ｉの位置をそれぞれ仮想放出点とし、各仮想放出点Ｐｖｉの放出強度をそれぞれ対応する観測器で得た濃度情報に逆比例した値として逆問題を設定し、粒子法を用いた拡散の数値解析で時間を計測時刻から逆に進行させる逆解析に基づく評価によって発生源（放出地点および放出時刻）を推定することとし、差分手法を用いないで粒子法を使用したので、数値不安定性の問題を解消した拡散物質の発生源推定装置および発生源推定方法を提供することができる。

また、本実施形態の拡散物質の発生源推定装置および発生源推定方法では、推定部２７において、粒子濃度の分散が最小となる位置における粒子濃度の逆数値を放出地点の放出強度と推定するので、放出地点および放出時刻のみならず、放出強度も推定することができる。

さらに、本実施形態の拡散物質の発生源推定装置および発生源推定方法では、評価点設定部２５により相互の離間距離が一定の仮想格子上で各格子線が交差する位置を評価点Ｐｏｅｊ（ｊ＝１〜ｒ^２）として設定して、推定部２７において、評価点について各仮想放出点からの粒子濃度分布を評価し、粒子濃度の分散が最小となる位置および時刻をそれぞれ放出地点および放出時刻と推定するので、粒子濃度分布評価の処理が簡便となり、装置設計およびプログラム設計をより容易なものとすることができる。

〔第２実施形態〕
次に、本発明の第２実施形態に係る拡散物質の発生源推定装置および発生源推定方法について説明する。本実施形態における拡散物質の発生源推定装置の構成は第１実施形態の構成（図１参照）と同等であり、各構成要素についての詳細な説明は省略する。

但し、第１実施形態では１面の仮想格子上で評価点Ｐｏｅｊ（ｊ＝１〜ｒ^２）を設定したのに対し、本実施形態では、格子線の相互の離間距離が段階的に異なるｓ面（ｓは正整数）の仮想格子（一面の仮想格子における格子線の相互の離間距離は一定）上で、それぞれ評価点ＰｏｅＮｊ（Ｎ＝１〜ｓ，ｊ＝１〜ｐ；ｐはｐ＜ｒ^２で正整数）を設定する点が異なる。なお、仮想格子の面の適用は、格子線の相互の離間距離が相対的に長い仮想格子から相対的に短い仮想格子へ段階的に使用していくものとし、ｓ段目の仮想格子の面が（離間距離が最短の）最小グリッドであるものとする。

次に、本実施形態における拡散物質の発生源推定方法について、図６〜図８を参照して説明する。ここで、図６は本実施形態に係る拡散物質の発生源推定方法を説明するフローチャートであり、図７は仮想格子に基づく評価点を用いた放出候補地点の推定を例示する説明図であり、図８は仮想格子のネスティングを説明する説明図である。

まず、ステップＳ２０１では、観測情報入手部２１において、観測器５−ｉ（ｉ＝１〜ｎ）の位置情報、該観測器５−ｉが計測した濃度情報および計測時刻情報を、通信インタフェース１１を介して入手する。

次に、ステップＳ２０２では、逆問題設定部２３において、各観測器５−ｉ（ｉ＝１〜ｎ）の位置をそれぞれ仮想放出点Ｐｖｉ（ｉ＝１〜ｎ）とし、各仮想放出点の放出強度（粒子強度）をそれぞれ対応する観測器５−ｉで得た濃度情報に逆比例した値として逆問題を設定する。

次に、ステップＳ２０３では、評価点設定部２５により、Ｎ段目の仮想格子を想定して、該仮想格子の領域上で各格子線が交差する位置を評価点ＰｏｅＮｊ（Ｎ＝１〜ｓ，ｊ＝１〜ｐ）として設定する。

次に、ステップＳ２０４では、逆解析部２６により、各仮想放出点からの放出について、粒子法を用いた拡散の数値解析で時間を計測時刻から逆に進行させて逆解析を行う。なお、観測器５−ｉ（ｉ＝１〜ｎ）それぞれの観測時刻が異なる場合には、最も遅い計測時刻の観測器に対応する仮想放出点から放出を開始し、その後、最も遅い計測時刻との相対時刻の短い順に、最も遅い計測時刻の観測器に対応する仮想放出点からの放出タイミングから該相対時刻遅れたタイミングで、他の仮想放出点から順次放出させていく。

次に、ステップＳ２０５では、推定部２７により、時間刻みΔｔずつ逆に進行させたそれぞれの時刻について、各仮想放出点Ｐｖｉ（ｉ＝１〜ｎ）からの粒子濃度分布を評価し、粒子濃度が略等しい領域の内、粒子濃度の分散が最小となる評価点Ｐｏｅｊの位置および時刻をそれぞれ放出地点および放出時刻と推定する。すなわち、各評価点Ｐｏｅｊに対する各仮想放出点からの粒子濃度分布の分散σ_ｊが最小となる評価点Ｐｏｅｊの位置を求め、粒子濃度の分散が最小となる評価点Ｐｏｅｊの位置および時刻をそれぞれ放出候補地点および放出時刻と推定し、また、粒子濃度の分散が最小となる評価点Ｐｏｅｊの位置における粒子濃度の逆数値を該放出候補地点の放出強度と推定する。

次に、ステップＳ２０６では、現在適用している仮想格子の面が最小グリッド（ｓ段目で）であるか否かを判断し、最小グリッドである場合には終了し、また最小グリッドでなければ、ステップＳ２０７に進んでＮをインクリメントした後、ステップＳ２０３に戻る。

ここで、図７を参照して、仮想格子に基づく評価点を用いた放出候補地点の推定について説明する。図７では、図４と同様に、観測器５−１〜５−３で濃度が計測され、各観測器の位置をそれぞれ仮想放出点Ｐｖ１〜Ｐｖ３として各仮想放出点から逆拡散させて放出候補地点を推定している。なお、図７では、本来の放出地点Ｐｏに対して放出候補地点Ｐｏｅが推定されている。同図から推察できるように、１面の仮想格子上で評価点Ｐｏｅｊを設定する場合に、放出地点の推定精度を向上させるためには、格子線の離間距離を短くして評価点数を多く設定する必要があるが、その分だけ計算量は増大することになる。そこで、本実施形態では、仮想格子のネスティング手法を用いて計算量をより低減することとした。

仮想格子のネスティングは、図８に示すように、格子線の相互の離間距離が相対的に長い仮想格子から相対的に短い仮想格子へと、段階的に行われる。図８では、離間距離が１０００［ｍ］の仮想格子Ｇ１と、離間距離が５００［ｍ］の仮想格子Ｇ２を例示している。より具体的に、格子線の離間距離は、例えば、１０００［ｍ］→５００［ｍ］→２５０［ｍ］→１２５［ｍ］→６７．５［ｍ］→３１．３［ｍ］→１５．６［ｍ］→７．８［ｍ］→３．９［ｍ］→１．９［ｍ］→１．０［ｍ］というふうに段階的に変えていき、放出候補地点Ｐｏｅが本来の放出地点Ｐｏにより近づくよう絞り込んでいく。

なお、Ｎがインクリメントされた（ステップＳ２０６）後、Ｎ＋１段目の仮想格子を用いて評価点を設定する場合に、Ｎ＋１段目の仮想格子は、前段（Ｎ段目）の仮想格子上で放出候補地点とされた点を、その仮想格子面内に含むように設定される。

このように、格子線の相互の離間距離が相対的に長い仮想格子から相対的に短い仮想格子へと、段階的に、（離間距離が最短で）最小グリッドとなるｓ段目の仮想格子の面まで適用していき、ｓ段目の仮想格子の面の適用において、推定部２７により、粒子濃度の分散が最小となる評価点Ｐｏｅｊの位置および時刻がそれぞれ放出候補地点および放出時刻と推定され、また、粒子濃度の分散が最小となる評価点Ｐｏｅｊの位置における粒子濃度の逆数値が該放出候補地点の放出強度と推定されたとき、該放出候補地点が放出地点として推定されることとなる。

以上説明したように、本実施形態の拡散物質の発生源推定装置および発生源推定方法では、評価点設定部２５において、推定部２７で推定された放出地点を含み、且つ相互の離間距離がより短い評価格子上で各格子線が交差する位置を評価点として再設定することとし、格子線の離間距離が相対的に長い仮想格子から相対的に短い仮想格子へと絞り込みながら評価点を再設定して発生源を推定するようにしている。これにより、第１実施形態（離間距離を最短とした１面の仮想格子上で評価点を設定して発生源を推定する場合）と比べて、１面の仮想格子上での評価点数を格段に小さくすることができ、全体的な処理の計算量を抑制してより速い発生源の推定を行うことができる。

以上、本発明の実施形態について図面を参照して詳述したが、本発明はこれら実施形態に限定されるものではなく、本発明の要旨を逸脱しない範囲の設計変更等があっても本発明に含まれる。

３ 拡散物質の発生源推定装置
５−１〜５−ｎ 観測器
１１ 通信インタフェース
１３ 入力部
１６ 発生源推定処理部
１７ 記憶部
１９ 出力部
２１ 観測情報入手部
２３ 逆問題設定部
２５ 評価点設定部
２６ 逆解析部
２７ 推定部

Claims (10)

1. 複数個の観測器からの情報に基づき、ガス発生源情報を推定する拡散物質の発生源推定装置であって、
   前記観測器から位置情報、計測した濃度情報および計測時刻を入手する観測情報入手手段と、
   各観測器の位置をそれぞれ仮想放出点とし、各仮想放出点の放出強度をそれぞれ対応する観測器で得た濃度情報に逆比例した値とする逆問題設定手段と、
   各仮想放出点からの放出について、粒子法を用いた拡散の数値解析で時間を前記計測時刻から逆に進行させて逆解析を行う逆解析手段と、
   逆に進行させたそれぞれの時刻において、各仮想放出点からの粒子濃度分布を評価し、粒子濃度が略等しい領域の内、粒子濃度の分散が最小となる位置および時刻をそれぞれ放出地点および放出時刻と推定する推定手段と、
   を有することを特徴とする拡散物質の発生源推定装置。
2. 前記推定手段は、粒子濃度の分散が最小となる位置における粒子濃度の逆数値を前記放出地点の放出強度と推定することを特徴とする請求項１に記載の拡散物質の発生源推定装置。
3. 前記逆解析手段は、最も遅い計測時刻との相対時刻で各仮想放出点から順次放出させていくことを特徴とする請求項１または請求項２に記載の拡散物質の発生源推定装置。
4. 相互の離間距離が一定の評価格子上で各格子線が交差する位置を評価点として設定する評価点設定手段、を有し、
   前記推定手段は、逆に進行させたそれぞれの時刻において、設定された評価点について各仮想放出点からの粒子濃度分布を評価し、粒子濃度の分散が最小となる位置および時刻をそれぞれ放出地点および放出時刻と推定することを特徴とする請求項１〜請求項３の何れか１項に記載の拡散物質の発生源推定装置。
5. 前記評価点設定手段は、前記推定手段で推定された放出地点を含み、且つ相互の離間距離がより短い評価格子上で各格子線が交差する位置を評価点として再設定することを特徴とする請求項４に記載の拡散物質の発生源推定装置。
6. 複数個の観測器からの情報に基づき、ガス発生源情報を推定する拡散物質の発生源推定方法であって、
   前記観測器から位置情報、計測した濃度情報および計測時刻を入手する観測情報入手ステップと、
   各観測器の位置をそれぞれ仮想放出点とし、各仮想放出点の放出強度をそれぞれ対応する観測器で得た濃度情報に逆比例した値とする逆問題設定ステップと、
   各仮想放出点からの放出について、粒子法を用いた拡散の数値解析で時間を前記計測時刻から逆に進行させて逆解析を行う逆解析ステップと、
   逆に進行させたそれぞれの時刻において、各仮想放出点からの粒子濃度分布を評価し、粒子濃度が略等しい領域の内、粒子濃度の分散が最小となる位置および時刻をそれぞれ放出地点および放出時刻と推定する推定ステップと、
   を有することを特徴とする拡散物質の発生源推定方法。
7. 前記推定ステップは、粒子濃度の分散が最小となる位置における粒子濃度の逆数値を前記放出地点の放出強度と推定することを特徴とする請求項６に記載の拡散物質の発生源推定方法。
8. 前記逆解析ステップは、最も遅い計測時刻との相対時刻で各仮想放出点から順次放出させていくことを特徴とする請求項６または請求項７に記載の拡散物質の発生源推定方法。
9. 相互の離間距離が一定の評価格子上で各格子線が交差する位置を評価点として設定する評価点設定ステップ、を有し、
   前記推定ステップは、逆に進行させたそれぞれの時刻において、設定された評価点について各仮想放出点からの粒子濃度分布を評価し、粒子濃度の分散が最小となる位置および時刻をそれぞれ放出地点および放出時刻と推定することを特徴とする請求項６〜請求項８の何れか１項に記載の拡散物質の発生源推定方法。
10. 前記評価点設定ステップは、前記推定手段で推定された放出地点を含み、且つ相互の離間距離がより短い評価格子上で各格子線が交差する位置を評価点として再設定することを特徴とする請求項９に記載の拡散物質の発生源推定方法。
    

Images