カテゴリー別検索
書籍詳細ページ
出版
| タイトル | 不均質地盤の地下水モデリングにおける逆解析手法とその利用法の現状と展望(<特集>地盤工学への逆解析／データ同化の利用) | ||||
| 著者 | 増本 清 | ||||
| 出版 | 地盤工学会誌 Vol.65 No.10 No.717 | ||||
| ページ | 14〜17 | 発行 | 2017/10/01 | 文書ID | jk201707170010 |
| 内容 | 表示 報告不均質地盤の地下水モデリングにおける逆解析手法とその利用法の現状と展望Current State and Future Vision of Utilization of Inverse Modeling for Heterogeneous Aquifer増本島根大学大学院. は じ め に清（ますもときよし）総合理工学研究科. 順解析と逆解析の方法こうした問題に対処するために，筆者はこれまでに，地下水挙動を予測する上で，不均質な地盤を適切にモ 計測形態の工夫（多点非定常データ）◯ 計算アルゴリ◯デル化することが求められる。そのために多点で非定常 適切化ズムの効率化（ adjoint 法，準ニュートン法)2) ◯水圧挙動を計測し，数値モデルを用いて逆解析すること法（ノルム最小法，スムージング法）などを試みてきた。が有効と思われる。しかし，不均質性に対応するためにこれにより，多様な未知パラメータ（透水性だけでなく，は，膨大な未知パラメータを推定する必要から，計測・間隙率，初期条件や境界条件など）の設定，多様な目的解析コストが膨大になり，それに加えて逆問題の不適切関数の設定に対応できる数値逆解析法が実現でき，一定性の問題が生じる。それでも，ハードウエア及びソフトの有効性を示した。また，より複雑な二相流動や移流分ウエア技術の進展により，こうした逆解析が実用的なレ散現象とそれに伴う新たな未知パラメータ（毛管圧力曲ベルに近づいていると思われる。本稿では，筆者が行っ線，分散係数など）にも対応できることが示される3),4)。てきた不均質地盤の地下水モデル数値逆解析法の概要をまた，筆者は，逆解析結果が適切であるかどうかを判示し，現状の可能性と課題を示し，将来展望を述べる。定する方法を提案し検討した（増本5) ）。これはラグラ.計測法と解析法の要点. 想定する計測法非定常水圧データは，不均質な水理特性の影響を受けンジュ乗数法として定式化された問題に対して，ラグランジアンのヘッセ行列の固有値を調べることにより適切性を評価するものである。以下にこれらの解析法の要点を示す。支配方程式の離散化るため，それを逆解析することにより，水理特性分布がある程度推定できることが理論的には期待される。そこ定 常 あ る い は 非 定 常 流 動 で あ れ ， FDM あ る い はで，ここでは，パッカーで区切られたボーリング孔内のFEM であれ，支配方程式及び初期・境界条件を離散化区間で注水又は揚水を行い，多孔の複数深度の各区間ですると，（非）線形連立方程式が得られる。例えば，単多数の水圧を計測するハイドロパルステスト1)を想定す相流動の場合，非定常地下水流動の支配方程式は以下のる（図―）。このような試験形態による計測情報の増ようになる。大をはかったとしてもなお，実際に使用できる計測データは未知パラメータ量に比べて不足する上に，計測誤差やモデル誤差などの多様な誤差を含むため，適当な逆解R (x , y , z , t )≡析結果を得ることは一般には容易でない。特に逆問題の不適切性と，計算量（容量及び速度）が問題になる。() {}& qK－:･(:P＋rg:z ) ＋Q＝0&t BmB………………………………………………………(1)ここで，F間隙率，B容積係数，K浸透係数，m粘性係数，P間隙水圧，r水の密度，g重力加速度，z高さ（鉛直上向き正），Q流量である。また，R(x,y, z, t ）はこの支配方程式が空間変数 x, y, z と時間変数 tの関数であることを示している。差分法を用いて，式(1)を時間及び空間で離散化することにより，時間ステップ n ，空間ステップ i, j, k に対する差分方程式が得られる。得られた差分方程式は，R ni, j, k＝0 ………………………………………………(2)と書ける。これは境界条件及び初期条件も考慮にいれた式である。これを満たす間隙水圧を求めることが単相流動の順解析に相当する。さらに，時空間節点番号を通し番号で表現することにより，次のような一般的な形で書図―14ハイドロパルステスト概念図ける。地盤工学会誌，―（) 報Ri(X, P )＝0 ……………………………………………(3)告て，増本2)で示されている準ニュートン法を用いた。準ここで， i 時空間節点の通し番号， X 未知パラメーニュートン法では，目的関数を未知パラメータで微分しタベクトル，P(圧力等の）状態量ベクトルである。た値（以下，勾配とよぶ。）を求める必要がある。既往二相流動や移流分散方程式のように，格子点ごとに複文献2),3)に示されているように不均質性を表現するため数の支配方程式が適用される場合においても，離散化後には，浸透係数など多数の未知パラメータを扱うことにの方程式に通し番号をつけることにより式(3)と同様のなるため， adjoint 法のような高速勾配計算法を適用す表現が可能である。ることが有効である。次節で各未知パラメータに関する逆問題の定式化と適切化法勾配（微分値）を adjoint 法で求める方法を示す。ここで解くべき逆問題は，順解析で得られた圧力値と目的関数の微分値計算法観測で得られた圧力値の重み付き残差 2 乗和にペナルadjoint 法は通常変分法により説明されるが，離散化ティ項 G を加えた目的関数 J （式( 4)）を最小とする未した問題に対してはラグランジュ乗数法の考え方で説明知パラメータを求めることになる。することができる4)。具体的には，目的関数と制約条件∑∑W nm(P ncal, m－P nobs, m)＋GJ＝n……………(4)mここで m観測点番号，W 重み，Pcal計算圧力値，式を用いてラグランジュ関数 L を定義する。式(2)に対するラグランジュ関数は式( 6 )，より一般的な式( 3 )に対しては式(7)のようになる。Pobs観測圧力値， G ペナルティ項（後述）である。∑∑∑∑lnijk･R nijkL＝J＋すなわち，式( 3 )を拘束条件として式( 4 )を最小化するi未知パラメータを求める等式拘束条件付き最小化問題として浸透係数 k 及び初期条件 P 0 及び流量 Q （境界セルk………………………(6)nL(X, P, l)＝J(X, P)＋∑li･Ri(X, P )………(7)して定式化される。ここでは，単相流動問題例における未知パラメータとjiここで，l随伴ベクトルである。随伴ベクトル（ラグランジュ乗数ベクトル）は，式では外部との流入出量を表す）を未知パラメータとする( 2 )（又は式( 3 )）より順解析で状態量を求めたあと，場合を例として示す。実際には式中の任意のパラメータ随伴方程式を連立して解くことにより算出することがでを未知パラメータとして扱うことができる。ここで，浸きる。式( 2 )に対する随伴方程式は式( 8 )である。これ透係数と初期条件及び流量はいずれも空間的な位置に応は線形の連立方程式となる。じて定められる量であり，離散化モデルでは空間ステッd&R abc&L&J＝＋＝0 ………………(8)ld ･&P nijk &P nijk a, b, c, d abc &P nijk∑プごとに定められるものであることに注意されたい。Gは未知パラメータのみからなる関数であり，各項の係数ここで a, b, c, d はそれぞれ i, j, k, n に対応する一時変数を変えることにより未知パラメータに関する事前情報をである。用いた適切化の程度を調整することができる。一例としこの連立方程式の係数行列は，順解析で使用される係て，浸透係数のスムージング，浸透係数・初期条件・境数行列の転置となるので，順解析ができていればコーデ界流量にノルム最小法を適用する場合，次の形になる。ィングの手間はそれほど増えない。多くの問題では対称G＝GKsm＋GK ln＋GIC ln＋GBC ln ………………………(5)行列になるので同じサブルーチンをそのまま使用するこここで， GKsm 浸透係数スムージング項， GK ln 浸透とになる。非定常問題において順解析で時間ステップご係 数ノ ルム 最小 項， GIC ln 初 期条 件ノル ム最 小 項，とに状態量の空間分布を求める形であれば，随伴方程式GBC ln 境界条件ノルム最小項である。スムージングはについては終末条件を与えて逆時間ステップごとに同様空間的に隣接する浸透係数パラメータからなる二次差分の連立方程式を解く形になる。その際，終末条件は＋1………………………………………………(9)lNijk ＝0（曲率）の重み付き二乗和，ノルム最小項は，パラメータごとに設定される基準値とパラメータとの差の重み付とすればよい4) 。このことは式( 8 )において n ＝ N の場き二乗和である。ノルム最小法は未知パラメータの大き合の式を作り，式(9)と組み合わせることにより，ほかさ，スムージング法は未知パラメータの滑らかさに関すの n の式と同形にできることから分かる。終末時間スる制限を与える。これらはいわゆる正則化法に属するもテップ番号は N ではなく N ＋ 1 とする点に注意されたので，他にも多様なペナルティ項を作成できる。この多い。様性が実用上重要な意味をもつ。実際ノルムの取り方は未知パラメータ Xm に対する目的関数の微分値は式多様な形を選ぶことができ，異なるパラメータ間の関係( 10 )で表される。これは随伴ベクトルを用いて計算で式を組み込める。また，各パラメータに関して微分可能きる。4)な多様な変数変換が適用できる。例えば，浸透係数は対数を取ることにより 0 以上という制約を解消できる。&L&J＝＋&Xm &Xm&R∑li･&Xmi …………………………(10)i未知パラメータの多様な設定法は，モデリングの目的や以下に，式( 10 )による微分値計算式の具体例として，観測情報量に応じてエンジニアの主観を入れつつニーズ初期条件及び境界流量に関する勾配計算式を示す。とすり合わせながら決めることになろう。上記の非線形最適化問題を解くための逆解析手法としOctober, 2017初期条件に関する勾配は，順解析すなわち式(2)を解いて得られる圧力値と式(8)を解いて得られる随伴ベク15 報告トルを用いて，次式により算出することができる。&R dabc&L&Jldabc･ 0 ＝00 ＝0 ＋&P ijk &P ijk a, b, c, d&P ijk∑……………(11)逆解析推定値における状態量と乗数ベクトルが求められているので，これらを用いることにより，ラグランジュ関数のヘッセ行列（式(14)）を作ることができる5)。式( 11 )を離散化した式に適用すると，次のような簡潔222& L & L & L 2&p&p&u &p&lな式となる。( )&L Vijk･qijk 1 d＝･lijk･&P 0ijkDt ldP1B:2L＝………………………(12)ここで，V差分化したセルの体積，Dt時間ステップ&2L&u&p&2L&u2&2L&u&l………………………(14)&2L&2L &2L &l&p &l&u &l2 である。以上の計算法を用いることにより，未知パラメータ数に関わらず，順解析の 2 倍程度の計算時間でMasumoto5) では，このラグランジュ関数のヘッセ行列目的関数の勾配を算出することができる。流量に関する（ 2 階微分）の固有値を用いた信頼性指標を提案している。これにより，逆解析結果の適切性を評価できると考微分値は次のようになる。&L＝lnijk ……………………………………………(13)&Q nijkえられる。このヘッセ行列の固有値の計算には，ハウスホルダーこの式は，ラグランジュ乗数の実体が流量に関する微分法及び二分岐法を用いた。なお，負の固有値は不要なの値であることをも示している。これを逆用すれば，流量（未知パラメータの個数)で求めない。式(14)の次元は，に関する微分値を数値微分により計算すれば，随伴方程＋(等式制約条件式の個数）× 2 となるので，計算量は式が正しく解けたかの検証に使用できる。このようにな膨大になる。このヘッセ行列は正方行列になり，行列のるのは，制約条件式の流量 Q に関する微分値が 1 とな正の固有値数が制約条件式数と未知パラメータの総数のっ て い る た め で あ る 。 境 界 流 量に 関 す る 微 分値 は 式和と等しくなるとき，得られた解が適切であるといえる。(13)の節点番号が境界節点の場合に相当する。なお，式( 14 )の行列は，制約条件式の個数と同じ個数の負の固有値をもつことが示される。ヘッセ行列の正の非線形最小化以上に示すように，式(2)を解いて得られた状態量を使用して，式(8)を逆時間に解いて得られる乗数ベクト固有値の欠落数を NL とすると，欠落数は式( 15 )のように表せる。ルを用いると，式( 10 )により勾配が計算できる。式NL＝Nc＋Nm－Ne …………………………………(15)( 10 )は式( 2 )及び式( 8 )を解くのに比べて微小な計算量で求まる。したがって，未知パラメータの種類や個数にここで， Nc制約条件式（式( 2)）の個数， Nm 未知関わりなく，勾配計算時間は順解析 1 回と随伴方程式 1パラメータの個数， Ne L のヘッセ行列の正の固有値回を解く計算量となる。準ニュートン法などの勾配法に数である。より，この勾配が十分に 0 に近づくまで反復を繰り返逆解析結果の信頼性評価指標 E ( Masumoto )5) を， Lすことにより，逆解析が実行される。未知パラメータ数のヘッセ行列の正の固有値の最大値 emax 及び最小値 eminが大きくなる場合は準ニュートン法においても記憶容量及び欠落数 NL を用いて，式(16)としている。を節約する方法の使用が望ましい2)。adjoint法を適用すE＝る場合，実際には固定したくないパラメータを未知として扱うことにより，モデル誤差を防ぐことが期待できる。emin－NL …………………………………………(16)emaxただし，正の固有値が欠落する場合は emin は 0 であ未知パラメータの候補としては，物性値や流量などのほる。欠落数は解の一意性に関する情報を示し，固有値のかに，空間節点の座標や水分特性曲線のパラメータなど，比は解の安定性（感度）に関する情報を示すものと考え式(3)に含まれる状態量以外のもの全てが対象となり得られる。E が大きいほど，逆解析結果はより適切であるる。 adjoint 法を使用する限り，未知パラメータの種類ことが期待される。なお，欠落数が生じた場合は，本来をいくら増やしても勾配計算時間はほぼ不変である。は不適切問題といえる。しかし，実際の数値解析におい. 逆解析の信頼性評価指標ては，計算機誤差などの影響で正値と 0 が判別できな前節に示したように， adjoint 法では未知パラメータい場合も考えられる。したがって，欠落数がある場合にを増やしても勾配計算量が変わらないため，多数のパラついても相対的な評価指標として使用できる可能性はあメータを未知として逆解析を実施することができそうでると考えられる。ある。ただし，未知パラメータ数が多すぎると観測情報の量に限界があるため，逆問題が不適切になる可能性がある。その場合，ペナルティ項を適当に設定することに.数値逆解析の例. 多様な未知パラメータの逆解析例より適切化を図ることができる。ただし，逆解析結果が多様な未知パラメータ推定の例として，浸透特性と初適切な条件のもとで解かれたかどうか，その信頼性を評期・境界条件を同時に逆解析した例を紹介する。末光と価する必要性が出てくる。すなわち，適切化による未知増本6)では，浸透係数分布に加えて初期圧力分布と境界パラメータに対する縛りをどの程度に設定すればよいの流量を未知量とする逆解析を二次元（セル数20×13）かが問題となる。一方， adjoint 法を適用した場合は，数値実験で試みている。適切化法としてノルム最小法を16地盤工学会誌，―（) 報告用いており，各未知パラメータに対するノルム最小項のと考 え ら れる 。 adjoint 法 を 使 用す る 場 合， 未 知 パラ重み係数を調節することにより，3 種類の未知パラメーメータの種類や個数が増えても目的関数の微分値（勾配）タが適切に推定できることを示している。特に初期条件の計算量は変わらないので，あらゆるパラメータを未知として不均質性を考慮しないで作成した初期条件を既知として扱うことによるコストは大きくない。この場合はとして逆解析を実施すると浸透係数分布が正解値と大き未知パラメータの自由度をノルム最小法などにより抑えく異なることが指摘されており，初期条件を未知としてれば，適切性も解決できる可能性がある。しかし，未知扱うことの有効性が示されている。実データの解析におパラメータ数が増えすぎると，勾配法の計算速度や記憶いて，ノルム最小項の重みをどのように選べばよいかが容量の問題，信頼性評価計算量が増大する，といった問課題となる。題がある。そのほかに，ウェルモデルによる補正や，予. 逆解析結果の適切性評価例測誤差の評価なども実用上重要と考えられる。逆解析結果の信頼性評価の例として，スムージングパ今後の展開として，評価指標値が増大するように問題ラメータの違いが逆問題適切性に及ぼす影響を検討したを変換する方法の検討が考えられる。通常の行列計算に例を紹介する。評価指標の妥当性について，増本と仲おける CG 法のように条件数を 1 に近づけるような変換野7)に数値実験による検討例が示されている。これによを行い，高速で安定な計算を行う方法が威力を発揮するると，正解値の分かっている二次元数値モデルに対してことが知られているが，非線形逆問題においてもこのよスムージング項の係数を変化させたときの，浸透係数分うな変換法を開発することが期待される。これは“非線布推定結果と正解値との差に対する評価指標との関係が形問題の線形化”による効率化とみることもできる。そ良い相関性をもつことが示され，スムージング項の係数のためのツールとしてこのような評価指標が適用できるを決定する上で評価指標が適用できることが示されてい可能性がある。特に目的関数各項の重みづけや未知パラる。ただし，大規模問題に適用するには計算量が膨大にメータの変数変換は多様であり，実際に求めたいものなるという課題が残されている。.地下水流動数値逆解析法の現状の課題と展望（目的），前提とできること（事前情報），使用できるデータ（計測情報の質と量），実際の計測状況（単純化によるモデル誤差が許容できるか）等に応じて，総合的に適当な設定法を構築することが望まれる。順解析を行う場合は，パラメータの値を決めることがエンジニアの考えるべきことであった。しかし，逆解析の枠組みでは，パラメータ間の関係をどう設定するか，参1)目的関数を構成する各項の関係をどう設定するか，仮定したモデルや既知としたパラメータは適切かどうか，ということを検討することになる。特に，問題となるのは，モデル誤差の低減であろう。こうした問題を解決する上2)で，逆解析結果の信頼性を定量的に評価することは有効と考えられる。それに耐えうる計算効率の向上が求めら3)れる。また，予測誤差を評価したいというニーズもあるだろう。実用化に向けて，多様な解析を実行できるツールができつつあるものの，その発展をはかるだけでなく，4)使用法についての検討がより重要である。残された課題として以下が挙げられる。1)三次元大規模問題のための一層の計算効率化が5)必要。2)多相流動や移流分散現象等に対応するためには3)多様な未知パラメータの適切な設定法が不明。4)多様な目的関数の適切な設定法が不明。5)信頼性評価の計算量が大きい。計測及び未知パラメータの限定法の検討が必要。実用化のためには，逆解析結果の信頼性と計算速度を高めることが重要であろう。例えば，モデル誤差を解消するためには，極力パラメータを固定しないことが有効October, 20176)7)考文献Masumoto, K., Tosaka, H., Kojima, K., Itoh, K. andOtsuka, Y.: New Measuring System and High SpeedThree Dimensional Inversion Method for HydropulseTomography., Proc. Int. Congress on Rock Mechanics,ISRM, pp. 847850, 1995.増本 清地下水理逆解析における計算機記憶容量節約アルゴリズム，応用力学論文集，Vol. 7, No. 1, pp. 191～199, 2004.増本 清地下二相流動モデルにおける相対浸透率曲線の同定法，日本地下水学会， 2008 年秋季講演会講演要旨，pp. 28～31, 2008.増本 清逆解析による水理特性パラメータの評価，日本地下水学会原位置トレーサー試験に関するワーキンググループ編「地下水のトレーサー試験」第 8 章 2 節，技報堂出版，pp. 195～212, 2009.Kiyoshi Masumoto: Credibility evaluation of numericallyestimated heterogeneous hydraulic property by eigenvalues of Hessian of Lagrange function for constrainedgroundwater problem, Proc. ModelCARE2007, IAHSPublication 320, pp. 5863, 2008.末光明信・増本 清地下水理物性および初期・境界条件の同時逆解析における適切化手法の検討，日本地下水学 会 ， 2015 年 秋 季 講 演 会 講 演 予 稿 ， pp. 102 ～ 107,2015.増本 清・仲野允浩地下流動特性逆解析における適切な計測形態とスムージングパラメータの検討，地盤と建設，Vol. 28, No. 1, pp. 153～160, 2010.（原稿受理2017.7.3)17 | ||||
| ログイン | |||||
- タイトル
- 不均質地盤の地下水モデリングにおける逆解析手法とその利用法の現状と展望(<特集>地盤工学への逆解析／データ同化の利用)
- 著者
- 増本 清
- 出版
- 地盤工学会誌 Vol.65 No.10 No.717
- ページ
- 14〜17
- 発行
- 2017/10/01
- 文書ID
- jk201707170010
- 内容
- 報告不均質地盤の地下水モデリングにおける逆解析手法とその利用法の現状と展望Current State and Future Vision of Utilization of Inverse Modeling for Heterogeneous Aquifer増本島根大学大学院. は じ め に清（ますもときよし）総合理工学研究科. 順解析と逆解析の方法こうした問題に対処するために，筆者はこれまでに，地下水挙動を予測する上で，不均質な地盤を適切にモ 計測形態の工夫（多点非定常データ）◯ 計算アルゴリ◯デル化することが求められる。そのために多点で非定常 適切化ズムの効率化（ adjoint 法，準ニュートン法)2) ◯水圧挙動を計測し，数値モデルを用いて逆解析すること法（ノルム最小法，スムージング法）などを試みてきた。が有効と思われる。しかし，不均質性に対応するためにこれにより，多様な未知パラメータ（透水性だけでなく，は，膨大な未知パラメータを推定する必要から，計測・間隙率，初期条件や境界条件など）の設定，多様な目的解析コストが膨大になり，それに加えて逆問題の不適切関数の設定に対応できる数値逆解析法が実現でき，一定性の問題が生じる。それでも，ハードウエア及びソフトの有効性を示した。また，より複雑な二相流動や移流分ウエア技術の進展により，こうした逆解析が実用的なレ散現象とそれに伴う新たな未知パラメータ（毛管圧力曲ベルに近づいていると思われる。本稿では，筆者が行っ線，分散係数など）にも対応できることが示される3),4)。てきた不均質地盤の地下水モデル数値逆解析法の概要をまた，筆者は，逆解析結果が適切であるかどうかを判示し，現状の可能性と課題を示し，将来展望を述べる。定する方法を提案し検討した（増本5) ）。これはラグラ.計測法と解析法の要点. 想定する計測法非定常水圧データは，不均質な水理特性の影響を受けンジュ乗数法として定式化された問題に対して，ラグランジアンのヘッセ行列の固有値を調べることにより適切性を評価するものである。以下にこれらの解析法の要点を示す。支配方程式の離散化るため，それを逆解析することにより，水理特性分布がある程度推定できることが理論的には期待される。そこ定 常 あ る い は 非 定 常 流 動 で あ れ ， FDM あ る い はで，ここでは，パッカーで区切られたボーリング孔内のFEM であれ，支配方程式及び初期・境界条件を離散化区間で注水又は揚水を行い，多孔の複数深度の各区間ですると，（非）線形連立方程式が得られる。例えば，単多数の水圧を計測するハイドロパルステスト1)を想定す相流動の場合，非定常地下水流動の支配方程式は以下のる（図―）。このような試験形態による計測情報の増ようになる。大をはかったとしてもなお，実際に使用できる計測データは未知パラメータ量に比べて不足する上に，計測誤差やモデル誤差などの多様な誤差を含むため，適当な逆解R (x , y , z , t )≡析結果を得ることは一般には容易でない。特に逆問題の不適切性と，計算量（容量及び速度）が問題になる。() {}& qK－:･(:P＋rg:z ) ＋Q＝0&t BmB………………………………………………………(1)ここで，F間隙率，B容積係数，K浸透係数，m粘性係数，P間隙水圧，r水の密度，g重力加速度，z高さ（鉛直上向き正），Q流量である。また，R(x,y, z, t ）はこの支配方程式が空間変数 x, y, z と時間変数 tの関数であることを示している。差分法を用いて，式(1)を時間及び空間で離散化することにより，時間ステップ n ，空間ステップ i, j, k に対する差分方程式が得られる。得られた差分方程式は，R ni, j, k＝0 ………………………………………………(2)と書ける。これは境界条件及び初期条件も考慮にいれた式である。これを満たす間隙水圧を求めることが単相流動の順解析に相当する。さらに，時空間節点番号を通し番号で表現することにより，次のような一般的な形で書図―14ハイドロパルステスト概念図ける。地盤工学会誌，―（) 報Ri(X, P )＝0 ……………………………………………(3)告て，増本2)で示されている準ニュートン法を用いた。準ここで， i 時空間節点の通し番号， X 未知パラメーニュートン法では，目的関数を未知パラメータで微分しタベクトル，P(圧力等の）状態量ベクトルである。た値（以下，勾配とよぶ。）を求める必要がある。既往二相流動や移流分散方程式のように，格子点ごとに複文献2),3)に示されているように不均質性を表現するため数の支配方程式が適用される場合においても，離散化後には，浸透係数など多数の未知パラメータを扱うことにの方程式に通し番号をつけることにより式(3)と同様のなるため， adjoint 法のような高速勾配計算法を適用す表現が可能である。ることが有効である。次節で各未知パラメータに関する逆問題の定式化と適切化法勾配（微分値）を adjoint 法で求める方法を示す。ここで解くべき逆問題は，順解析で得られた圧力値と目的関数の微分値計算法観測で得られた圧力値の重み付き残差 2 乗和にペナルadjoint 法は通常変分法により説明されるが，離散化ティ項 G を加えた目的関数 J （式( 4)）を最小とする未した問題に対してはラグランジュ乗数法の考え方で説明知パラメータを求めることになる。することができる4)。具体的には，目的関数と制約条件∑∑W nm(P ncal, m－P nobs, m)＋GJ＝n……………(4)mここで m観測点番号，W 重み，Pcal計算圧力値，式を用いてラグランジュ関数 L を定義する。式(2)に対するラグランジュ関数は式( 6 )，より一般的な式( 3 )に対しては式(7)のようになる。Pobs観測圧力値， G ペナルティ項（後述）である。∑∑∑∑lnijk･R nijkL＝J＋すなわち，式( 3 )を拘束条件として式( 4 )を最小化するi未知パラメータを求める等式拘束条件付き最小化問題として浸透係数 k 及び初期条件 P 0 及び流量 Q （境界セルk………………………(6)nL(X, P, l)＝J(X, P)＋∑li･Ri(X, P )………(7)して定式化される。ここでは，単相流動問題例における未知パラメータとjiここで，l随伴ベクトルである。随伴ベクトル（ラグランジュ乗数ベクトル）は，式では外部との流入出量を表す）を未知パラメータとする( 2 )（又は式( 3 )）より順解析で状態量を求めたあと，場合を例として示す。実際には式中の任意のパラメータ随伴方程式を連立して解くことにより算出することがでを未知パラメータとして扱うことができる。ここで，浸きる。式( 2 )に対する随伴方程式は式( 8 )である。これ透係数と初期条件及び流量はいずれも空間的な位置に応は線形の連立方程式となる。じて定められる量であり，離散化モデルでは空間ステッd&R abc&L&J＝＋＝0 ………………(8)ld ･&P nijk &P nijk a, b, c, d abc &P nijk∑プごとに定められるものであることに注意されたい。Gは未知パラメータのみからなる関数であり，各項の係数ここで a, b, c, d はそれぞれ i, j, k, n に対応する一時変数を変えることにより未知パラメータに関する事前情報をである。用いた適切化の程度を調整することができる。一例としこの連立方程式の係数行列は，順解析で使用される係て，浸透係数のスムージング，浸透係数・初期条件・境数行列の転置となるので，順解析ができていればコーデ界流量にノルム最小法を適用する場合，次の形になる。ィングの手間はそれほど増えない。多くの問題では対称G＝GKsm＋GK ln＋GIC ln＋GBC ln ………………………(5)行列になるので同じサブルーチンをそのまま使用するこここで， GKsm 浸透係数スムージング項， GK ln 浸透とになる。非定常問題において順解析で時間ステップご係 数ノ ルム 最小 項， GIC ln 初 期条 件ノル ム最 小 項，とに状態量の空間分布を求める形であれば，随伴方程式GBC ln 境界条件ノルム最小項である。スムージングはについては終末条件を与えて逆時間ステップごとに同様空間的に隣接する浸透係数パラメータからなる二次差分の連立方程式を解く形になる。その際，終末条件は＋1………………………………………………(9)lNijk ＝0（曲率）の重み付き二乗和，ノルム最小項は，パラメータごとに設定される基準値とパラメータとの差の重み付とすればよい4) 。このことは式( 8 )において n ＝ N の場き二乗和である。ノルム最小法は未知パラメータの大き合の式を作り，式(9)と組み合わせることにより，ほかさ，スムージング法は未知パラメータの滑らかさに関すの n の式と同形にできることから分かる。終末時間スる制限を与える。これらはいわゆる正則化法に属するもテップ番号は N ではなく N ＋ 1 とする点に注意されたので，他にも多様なペナルティ項を作成できる。この多い。様性が実用上重要な意味をもつ。実際ノルムの取り方は未知パラメータ Xm に対する目的関数の微分値は式多様な形を選ぶことができ，異なるパラメータ間の関係( 10 )で表される。これは随伴ベクトルを用いて計算で式を組み込める。また，各パラメータに関して微分可能きる。4)な多様な変数変換が適用できる。例えば，浸透係数は対数を取ることにより 0 以上という制約を解消できる。&L&J＝＋&Xm &Xm&R∑li･&Xmi …………………………(10)i未知パラメータの多様な設定法は，モデリングの目的や以下に，式( 10 )による微分値計算式の具体例として，観測情報量に応じてエンジニアの主観を入れつつニーズ初期条件及び境界流量に関する勾配計算式を示す。とすり合わせながら決めることになろう。上記の非線形最適化問題を解くための逆解析手法としOctober, 2017初期条件に関する勾配は，順解析すなわち式(2)を解いて得られる圧力値と式(8)を解いて得られる随伴ベク15 報告トルを用いて，次式により算出することができる。&R dabc&L&Jldabc･ 0 ＝00 ＝0 ＋&P ijk &P ijk a, b, c, d&P ijk∑……………(11)逆解析推定値における状態量と乗数ベクトルが求められているので，これらを用いることにより，ラグランジュ関数のヘッセ行列（式(14)）を作ることができる5)。式( 11 )を離散化した式に適用すると，次のような簡潔222& L & L & L 2&p&p&u &p&lな式となる。( )&L Vijk･qijk 1 d＝･lijk･&P 0ijkDt ldP1B:2L＝………………………(12)ここで，V差分化したセルの体積，Dt時間ステップ&2L&u&p&2L&u2&2L&u&l………………………(14)&2L&2L &2L &l&p &l&u &l2 である。以上の計算法を用いることにより，未知パラメータ数に関わらず，順解析の 2 倍程度の計算時間でMasumoto5) では，このラグランジュ関数のヘッセ行列目的関数の勾配を算出することができる。流量に関する（ 2 階微分）の固有値を用いた信頼性指標を提案している。これにより，逆解析結果の適切性を評価できると考微分値は次のようになる。&L＝lnijk ……………………………………………(13)&Q nijkえられる。このヘッセ行列の固有値の計算には，ハウスホルダーこの式は，ラグランジュ乗数の実体が流量に関する微分法及び二分岐法を用いた。なお，負の固有値は不要なの値であることをも示している。これを逆用すれば，流量（未知パラメータの個数)で求めない。式(14)の次元は，に関する微分値を数値微分により計算すれば，随伴方程＋(等式制約条件式の個数）× 2 となるので，計算量は式が正しく解けたかの検証に使用できる。このようにな膨大になる。このヘッセ行列は正方行列になり，行列のるのは，制約条件式の流量 Q に関する微分値が 1 とな正の固有値数が制約条件式数と未知パラメータの総数のっ て い る た め で あ る 。 境 界 流 量に 関 す る 微 分値 は 式和と等しくなるとき，得られた解が適切であるといえる。(13)の節点番号が境界節点の場合に相当する。なお，式( 14 )の行列は，制約条件式の個数と同じ個数の負の固有値をもつことが示される。ヘッセ行列の正の非線形最小化以上に示すように，式(2)を解いて得られた状態量を使用して，式(8)を逆時間に解いて得られる乗数ベクト固有値の欠落数を NL とすると，欠落数は式( 15 )のように表せる。ルを用いると，式( 10 )により勾配が計算できる。式NL＝Nc＋Nm－Ne …………………………………(15)( 10 )は式( 2 )及び式( 8 )を解くのに比べて微小な計算量で求まる。したがって，未知パラメータの種類や個数にここで， Nc制約条件式（式( 2)）の個数， Nm 未知関わりなく，勾配計算時間は順解析 1 回と随伴方程式 1パラメータの個数， Ne L のヘッセ行列の正の固有値回を解く計算量となる。準ニュートン法などの勾配法に数である。より，この勾配が十分に 0 に近づくまで反復を繰り返逆解析結果の信頼性評価指標 E ( Masumoto )5) を， Lすことにより，逆解析が実行される。未知パラメータ数のヘッセ行列の正の固有値の最大値 emax 及び最小値 eminが大きくなる場合は準ニュートン法においても記憶容量及び欠落数 NL を用いて，式(16)としている。を節約する方法の使用が望ましい2)。adjoint法を適用すE＝る場合，実際には固定したくないパラメータを未知として扱うことにより，モデル誤差を防ぐことが期待できる。emin－NL …………………………………………(16)emaxただし，正の固有値が欠落する場合は emin は 0 であ未知パラメータの候補としては，物性値や流量などのほる。欠落数は解の一意性に関する情報を示し，固有値のかに，空間節点の座標や水分特性曲線のパラメータなど，比は解の安定性（感度）に関する情報を示すものと考え式(3)に含まれる状態量以外のもの全てが対象となり得られる。E が大きいほど，逆解析結果はより適切であるる。 adjoint 法を使用する限り，未知パラメータの種類ことが期待される。なお，欠落数が生じた場合は，本来をいくら増やしても勾配計算時間はほぼ不変である。は不適切問題といえる。しかし，実際の数値解析におい. 逆解析の信頼性評価指標ては，計算機誤差などの影響で正値と 0 が判別できな前節に示したように， adjoint 法では未知パラメータい場合も考えられる。したがって，欠落数がある場合にを増やしても勾配計算量が変わらないため，多数のパラついても相対的な評価指標として使用できる可能性はあメータを未知として逆解析を実施することができそうでると考えられる。ある。ただし，未知パラメータ数が多すぎると観測情報の量に限界があるため，逆問題が不適切になる可能性がある。その場合，ペナルティ項を適当に設定することに.数値逆解析の例. 多様な未知パラメータの逆解析例より適切化を図ることができる。ただし，逆解析結果が多様な未知パラメータ推定の例として，浸透特性と初適切な条件のもとで解かれたかどうか，その信頼性を評期・境界条件を同時に逆解析した例を紹介する。末光と価する必要性が出てくる。すなわち，適切化による未知増本6)では，浸透係数分布に加えて初期圧力分布と境界パラメータに対する縛りをどの程度に設定すればよいの流量を未知量とする逆解析を二次元（セル数20×13）かが問題となる。一方， adjoint 法を適用した場合は，数値実験で試みている。適切化法としてノルム最小法を16地盤工学会誌，―（) 報告用いており，各未知パラメータに対するノルム最小項のと考 え ら れる 。 adjoint 法 を 使 用す る 場 合， 未 知 パラ重み係数を調節することにより，3 種類の未知パラメーメータの種類や個数が増えても目的関数の微分値（勾配）タが適切に推定できることを示している。特に初期条件の計算量は変わらないので，あらゆるパラメータを未知として不均質性を考慮しないで作成した初期条件を既知として扱うことによるコストは大きくない。この場合はとして逆解析を実施すると浸透係数分布が正解値と大き未知パラメータの自由度をノルム最小法などにより抑えく異なることが指摘されており，初期条件を未知としてれば，適切性も解決できる可能性がある。しかし，未知扱うことの有効性が示されている。実データの解析におパラメータ数が増えすぎると，勾配法の計算速度や記憶いて，ノルム最小項の重みをどのように選べばよいかが容量の問題，信頼性評価計算量が増大する，といった問課題となる。題がある。そのほかに，ウェルモデルによる補正や，予. 逆解析結果の適切性評価例測誤差の評価なども実用上重要と考えられる。逆解析結果の信頼性評価の例として，スムージングパ今後の展開として，評価指標値が増大するように問題ラメータの違いが逆問題適切性に及ぼす影響を検討したを変換する方法の検討が考えられる。通常の行列計算に例を紹介する。評価指標の妥当性について，増本と仲おける CG 法のように条件数を 1 に近づけるような変換野7)に数値実験による検討例が示されている。これによを行い，高速で安定な計算を行う方法が威力を発揮するると，正解値の分かっている二次元数値モデルに対してことが知られているが，非線形逆問題においてもこのよスムージング項の係数を変化させたときの，浸透係数分うな変換法を開発することが期待される。これは“非線布推定結果と正解値との差に対する評価指標との関係が形問題の線形化”による効率化とみることもできる。そ良い相関性をもつことが示され，スムージング項の係数のためのツールとしてこのような評価指標が適用できるを決定する上で評価指標が適用できることが示されてい可能性がある。特に目的関数各項の重みづけや未知パラる。ただし，大規模問題に適用するには計算量が膨大にメータの変数変換は多様であり，実際に求めたいものなるという課題が残されている。.地下水流動数値逆解析法の現状の課題と展望（目的），前提とできること（事前情報），使用できるデータ（計測情報の質と量），実際の計測状況（単純化によるモデル誤差が許容できるか）等に応じて，総合的に適当な設定法を構築することが望まれる。順解析を行う場合は，パラメータの値を決めることがエンジニアの考えるべきことであった。しかし，逆解析の枠組みでは，パラメータ間の関係をどう設定するか，参1)目的関数を構成する各項の関係をどう設定するか，仮定したモデルや既知としたパラメータは適切かどうか，ということを検討することになる。特に，問題となるのは，モデル誤差の低減であろう。こうした問題を解決する上2)で，逆解析結果の信頼性を定量的に評価することは有効と考えられる。それに耐えうる計算効率の向上が求めら3)れる。また，予測誤差を評価したいというニーズもあるだろう。実用化に向けて，多様な解析を実行できるツールができつつあるものの，その発展をはかるだけでなく，4)使用法についての検討がより重要である。残された課題として以下が挙げられる。1)三次元大規模問題のための一層の計算効率化が5)必要。2)多相流動や移流分散現象等に対応するためには3)多様な未知パラメータの適切な設定法が不明。4)多様な目的関数の適切な設定法が不明。5)信頼性評価の計算量が大きい。計測及び未知パラメータの限定法の検討が必要。実用化のためには，逆解析結果の信頼性と計算速度を高めることが重要であろう。例えば，モデル誤差を解消するためには，極力パラメータを固定しないことが有効October, 20176)7)考文献Masumoto, K., Tosaka, H., Kojima, K., Itoh, K. andOtsuka, Y.: New Measuring System and High SpeedThree Dimensional Inversion Method for HydropulseTomography., Proc. Int. Congress on Rock Mechanics,ISRM, pp. 847850, 1995.増本 清地下水理逆解析における計算機記憶容量節約アルゴリズム，応用力学論文集，Vol. 7, No. 1, pp. 191～199, 2004.増本 清地下二相流動モデルにおける相対浸透率曲線の同定法，日本地下水学会， 2008 年秋季講演会講演要旨，pp. 28～31, 2008.増本 清逆解析による水理特性パラメータの評価，日本地下水学会原位置トレーサー試験に関するワーキンググループ編「地下水のトレーサー試験」第 8 章 2 節，技報堂出版，pp. 195～212, 2009.Kiyoshi Masumoto: Credibility evaluation of numericallyestimated heterogeneous hydraulic property by eigenvalues of Hessian of Lagrange function for constrainedgroundwater problem, Proc. ModelCARE2007, IAHSPublication 320, pp. 5863, 2008.末光明信・増本 清地下水理物性および初期・境界条件の同時逆解析における適切化手法の検討，日本地下水学 会 ， 2015 年 秋 季 講 演 会 講 演 予 稿 ， pp. 102 ～ 107,2015.増本 清・仲野允浩地下流動特性逆解析における適切な計測形態とスムージングパラメータの検討，地盤と建設，Vol. 28, No. 1, pp. 153～160, 2010.（原稿受理2017.7.3)17
- ログイン