第 26 回数値流体力学シンポジウム
講演番号 E08-3

逆拡散問題の逆ランダム．ウォーク法による計算の試み
A numerical method of reverse diffusion
日野幹雄 藤沢市片瀬山 4-20-6
Mikio Hino, Prof. Emeritus, Tokyo Institute of Technology, Katase-Yanma 4-20-6,Fujisawa-Shi

A method of the inverse source estimation is proposed, using inverse (random) walk. A numerical simulation of normal diffusion has been performed by usual random walk method. The flow
field is assumed to be given or known. A problem of reverse diffusion from the particle distribution obtained at a down stream section has been examined which proceeds upstream withhigher biased probability towards the local peak of concentration. Two examples are shown fora uniform and an undulating flow fields. The proposed method is very simple, but farther improvements are necessary to determine the reverse transition probability.
1.はじめに最近，八木．加藤．半場(2011)1) は，偏微分方程式により拡散の逆問題を解く試みを行っている．著者はこの発表を聞いたとき以前に似たような問題を解いたことをお思い出した．東工大定年後に移った私大文系学部(中央大総合政策学部)で文系に相応しい仕事をやろうと思い「数値文明学」とでも呼ぶべき論文を 8 編ほど書いた．
その中の一つに「草原の砂漠化と放牧」2) の研究をした．
このなかで，羊群が豊かな草を探し移動しながらお互いが付きすぎず，また群れ集団から離れすぎずに纏まろうとする羊の本能による「羊の群れの離合集散」の問題(日野 2) である．多数の羊が群れをなして，草をはみながら，自分の周りの草を食べ尽くすと移動する．草をはんでいる間は隣どうしが近づき過ぎないように散らばろうとするが，しかし，仲間の羊から離れ過ぎると，仲間の羊群へ集まろうとする．この関係を，羊の密度がある一定値より高いときは正の拡散係数の拡散とみなし，羊の密度がその一定値より低くなると，負の拡散係数を採る拡散とみなして，2 つの方法(偏微分方程式とランダム．ウォーク+個別要素法)で数値的に解いたことがある．偏微分方程式の場合には，拡散係数が負の値を採る逆拡散過程
では数値解は発散し，計算を収束させるのに大変苦労した．しかし，ランダム．ウォーク+個別要素法的に，個々の羊どうしの間に距離に応じた斥力や引力が働くとしたシミュレーションでは，いとも簡単に計算ができた．この論文は中央大学総合政策学部の紀要である「総合政策研究」に発表しただけであるので，ここに，その結果の一部を紹介する．
まず，偏微分方程式による解で，Fig.1 は敷地内に一様に広がった羊の群れが，凝集本能により集まる(逆拡散)が，ある密度以上になると互いの間隔を保とうとする本能により結局はある大きさの塊となる様子のシミュレーション結果である．また，Fig.2 は，草をはみながら周りの草が無くなると集団を保ちながら移動する様子のシミュレーション結果である．この偏微分方程式による解法では，逆拡散過程で収斂せず苦労した．10 年以上前の研究なので，その時どんなアルゴリズムを考えたのか思い出せない．
一方，羊の行動を一匹づつ追う「個別要素法」では，そんな苦労はなく楽に計算出来た．Fig.3 は最初バラバラに散らばっていた羊群が徐々にある大きさの集団を形成する様子を，Fig.4 は草をはみながら周りの草が少なくなるとリーダー的羊に導かれるように移動する様子である．
この論文で発表する方法は，以前著者が「羊群の離合集散」に用いたものであるが，あまりにも単純でシンポジウムでの発表は止めようかと思ったが，文献調査では第 26 回数値流体力学シンポジウム
講演番号 E08-3

この種の方法を見つけられなかったので，一応報告をする．その他の逆拡散の最近の論文を文献 3) 以下に挙げる．
2. 前提など
1) 流れの場の時間空間的変化は，与えられているとする．すなわち，流れ場の順計算．逆計算は別の問題であるとする．流れ場を解くには，初期値と境界条件を与えれば順計算としては問題なく解けるが，逆計算についてはほとんど検討されていないと思う．実験的には層流に関してはあまり問題はないようである(Taylor の実験cinema)．
2) 1 点もしくは数点の放出源からの下流のある断面上の放出粒子の濃度が与えられたとき，この拡散場の最終断面の粒子分布から時間を遡って上流各点での粒子濃度分布を逆計算する．
3) 以下で扱う例題は 2 次元問題であるが，3 次元問題にもこの方法は容易に拡張しうる．
4) 以下の例題では，放出源が 1 つの場合を示したが，この方法は複数の放出源の場合にも適用可能である．
3. 順拡散のシミュレーション
放出源から放出された第 L グループの第 M 番粒子の時刻 N の( I,J)面上の位置をIp(L,M,N) および Jp(L,M,N)とし，各粒子の 1 時間ステップの J 方向への遷移確率
prob(Jp+1，prob(Jp),prob(Jp-1)
を与え，多数の粒子の拡散をランダム．ウォークによりシミュレートし，順拡散の濃度場を求めた．こうして得られた濃度分布図が Fig.5 で，放出粒子数は 4,000 個である．
4. 逆推定の方法
最下流断面の濃度分布を順拡散の結果で与え，流れ場は与えられているとして時間を遡って逆拡散計算をおこなう．実際問題として，最終断面の濃度分布が完全に得られている訳ではないから，その場合には内挿により濃度分布を与える．
時間的空間的に変動する流れ場が分かっているとして，
2 つの逆解法を試みた．
弟 1 の方法では，各格子点の粒子が各段面の最大濃度点の方向への逆遷移確率が，反対方向への逆遷移確率よりも流れの逆方向への逆遷移確率が一番高くする．
弟 2 の方法は，各粒子に放出点での放出時刻(L),同時に放出される粒子数(M)と現時刻(N)を指標として逆計第 26 回数値流体力学シンポジウム
講演番号 E08-3

Fig.5：順拡散(ランダムウオーク法)
算をする方法であるが，各粒子の上流側への遷移確率を順拡散と同じに採れば，拡散はそのまま継続するだけである．そこで，本方法では，濃度分布が与えられた断面や逆推定された濃度分布の極大点 Jc(I,J,N)を求める．逆拡散課程では，粒子の遷移確率をこの極大点の方向に向かって大きく採る．
Jp > Jc : prob(Jp+1)<prob(Jp-1)
Jp < Jc : prob(Jp+1)>prob(Jp-1)
この方法の問題点は逆遷移確率の最適値を推定することである．偏微分方程式による場合も拡散係数の推定が難しいのと同様である．しかし，種々遷移確率の値を変えて計算したが，之にはそれほど敏感ではない．
4.1 一様流の場合
逆流跡線法では，流れ場が一様流の場合濃度分布はそのままの分布が上流に移動するだけである．本方法では，濃度分布が与えられた断面や逆推定された濃度分布の極大値をまず求める．
すなわち，下流端の 1 断面で濃度分布を与え，1 ステップずつ遡って濃度場を上記の方法で逆推定計算し，これを合成したものが Fig.6 である．Fig.5 と fig.6 の濃度分布は，等濃度線は放出源へ近づくに連れて収束度が悪化しており，今後，アルゴリズムや遷移確率係数の検討が必要である．
4.2 クロスフローのある場合
流れ場が一様ではなく横流れがある場合には，各点の粒子が 1 ステップ前に，確率平均的にどの点から移動してきたかの逆算をする．流れ場が分かっていればこれは簡単にできる．それ以外は一様流場の計算と同じである．
横断方向にクロスフローV(x,y,t) =A(y,t)sin[(x-ct)/Lx](Lx:x 方向の波長)がある場合には，逆遷移位置に移流分を補正すればよい．
Jp(L,M,N-1)=Jp(L,M,N)+prob(Ip,Jp,N)―NINT(V(Ip(L,M,N),Jp(L,M.N)))＊Dt/Dy

4.2.1 クロスフローのある場合のアルゴリズムの検証
この方法による逆拡散のアルゴリズムが巧く適用出来るかを種々検討したが，Fig.7 はその一例で，原点で放出された 100 個の粒子が「拡散の無いクロスフロー」に乗って移動し最遠断面に到達し，そこから逆算して元の軌道上を逆に辿って放出源に戻るかの計算である．粒子群は拡散することなく蛇行する流れによりもとの点に戻っており，逆計算のアルゴリズムが巧く働くことを示し

Fig.6 ：逆拡散の計算結果

Fig.7：無拡散の場合の逆拡散(逆移流)―アルゴズムの検討ている．
4.2.2 計算例
クロスフローのある場合の順拡散結果を Fig.8 に，上述の方法による逆拡散推定の濃度分布図の結果を Fig.9に示す．本方法による逆拡散の計算結果は放出源近傍を除けばまずまずである．なお，Figs.8 & 9 とも濃度の等値線が塊状になっているのは，VDt に比べて Dx が同程度と粗いのと，流下方向の拡散を考慮していないためと,使用した CG ソフトの特性(欠陥)によるものである．(人間の目による作図では，滑らかな拡散等濃度曲線となる．)
5. 結びにかえて
5. 1 名称について
この方法の適切な名称が思いつかない．プログラムを作るときに順拡散を「ランダム．ウォーク法」で作ったので，逆計算を何となく「逆ランダム．ウォーク法」と呼んでいたが，現段階では，「逆流跡線法」＋「？法」＝「？？法」である．今後，スキームの改良に伴い適切な名称を考えたい．
5.2 雄蛾の雌蛾探知行動
最近知ったのだが，カイコガの雄が雌の発散するフェロモンを感知して(嗅覚器は鼻ではなく触角にある)雌蛾に近づく行動は，まず少しだけ雌のフェロモンの漂って来る方向に進み，次に前進」しながらそれと直角方向に往復行動をし，やがて円周運動をするということである．

第 26 回数値流体力学シンポジウム
講演番号 E08-3

Fig.8 ：Cross flow のある場合の順拡散

Fig.10：カイコガの雌蛾の発するフェロモンに対する雄蛾の探索行動 (NHK 高校講座生物より)(Fig.10) このジクザク行動はフェロモン濃度の最も高い位置を探しているのであろう．何やら，本論文の方法に似ている．しかし，蛾の場合はその場その場でのフェロモン濃度分布が感知できるが，本論文では出発断面での濃度分布が与えられるだけで，それより風上の濃度分布は自ら逆推定しなければならない点が難しい点である，
参考文献(1)八木 徹．北沢大輔．安部 論．半場藤弘．加藤信介 矩

Fig.9 ：Cross flow のある場合の逆拡散形湾における汚染源特定のためのリバースシミュレー シ ョ ン ． 第 25 回 数 値 流 体 シ ン ポ ジ ウ ム ，E12-5,1-4.(2011)(2) 日野幹雄 草原の都市化過放と砂漠化．牧総合政策
研究 7，33-54.中央大学総合政策学部, (2001)(3) Islam,M.A. Application of Gaussian plume model
to determine the location of an unknown
souce.Water, Air and Soil Pollution 112, 241-245.,(1999)(4) Islam,M.A. & Dubey,B.N. An improved method forthe identification of an undeclared transboundaryemission
souce.
Environment
Protection
Engineering. 33, 67-74. (2007)(5) Bady,M., Kato,S. & Hung,H. Identification ofpollution sources in urban areas using revese
simulation with reversed time marching nethod. J.
of Asian Architecture and Building Engineering,275-282 (2009)(6) Abe,S. & Kato,S. Study on improving thenumerical stability of an algorithm to identify
pollutantbsources by reverse simulation. The Fifth
International Symposium on Computational WindEngineering (CWE2010)