• 検索結果がありません。

On the Prediction of the Concentration Distribution of One-Dimensional Constant Coefficient Stochastic Convective-Dispersion Equation Based on the Second-Order Filter

N/A
N/A
Protected

Academic year: 2021

シェア "On the Prediction of the Concentration Distribution of One-Dimensional Constant Coefficient Stochastic Convective-Dispersion Equation Based on the Second-Order Filter"

Copied!
6
0
0

読み込み中.... (全文を見る)

全文

(1)

第33回 水理 講演 会論 文集 1989年2月

2次 フィル ターに よる定 係数一 次元確 率 移流 分散 方程式 の濃 度予 測 につい て

On the Prediction of the Concentration Distribution

of One-Dimensional Constant Coefficient Stochastic Convective-Dispersion Equation Based on the Second-Order Filter

九 州 大 学 大 学 院 吉 永 宙 司 Hiroshi YOSHINAGA 九 州 大 学 工 学 部 河 村 明 Akira KAWAMURA 九 州 大 学 工 学 部 神 野 健 二 Kenji JINNO 九 州 大 学 工 学 部 上 田 年 比 古 Toshihiko UEDA

Two methods, one uses the second-order filter and the other the extended Kalman filter are proposed, which are used to predict the concentration distribution and to identify the parameters of the one-dimensional stochastic convective-dispersion equation is discussed. It is shown that the proposed methods identify the parameters effectively and predict the concentration distribution accurately. It is clarified that the computational time for the second-order filter requires about two times longer than that

for the extended Kalman filter, whereas the accuracy of prediction of the second-order filter is almost the same as that of the extended Kalman filter.

Keyword: concentration distribution, convective-dispersion equation, second-order filter, parameter identification

1. は じめ に

移流分 散現 象 を考 える場合 、従 来、 確率 変動 項 は無 視 され確 定論 的 アプ ローチ がな され て、初 期条 件、境

界条件 お よび物理 パ ラ メー タは既 知情 報 として取扱 われ て いる。 また移 流分 散現 象の濃 度 の観測 デ ー タにつ

いて は、空 間的 に任意 に設 置 され た観 測 点 において経 時 観測 が行 われ てい る が、得 られた観 測情 報 は移流 分

散方程 式 の解析解 、数 値解 の 吟味 に使 わ れ る程 度 であ る1)。 著者 らは 、上述 の問題 の解決 への アプ ロ ーチ と

して、確率 移流 分散方 程式 を フー リエ級 数展 開 を用 い て フー リエ係数 に関す る連 立常 微分 方程式 に変 形 し、

次 いで、初 期条 件 、境界 条件 お よび物理 パ ラ メータす なわ ち物理 量(流 速 、分散係 数 、一 次反 応係数)は 不 明

で あ ると し、空 間 的 に任 意 に設置 され た観測 点 か らの定 時観 測 デ ータ を用 い て、各 時点毎 に物理 量 お よび フ

ー リエ係 数 を拡張 カルマ ンフィル ターに よ り逐次 同定 しなが ら濃 度予 測 を行 う手法 を提案 してい る2)∼4)。

本論 文 では さ らに、テ イ ラー級数 展開 の2次 項 まで を考慮 して 、 カル マ ン フィル ター の連 続 時間非 線形 状態

方 程式 をbias correctionと 呼 ばれ るベ ク トル外 力関 数で 補正す る2次 フ ィル ター5)∼7)によ り、定 係数 一 次

元確率移 流分 散方 程式 の濃度 予測 を行 う手 法 を定 式化 す る。 次い で、本手 法 お よび1次 フィル タ ー(拡張 カ

ルマ ン フィル タ-)を 模 擬発 生 デ ータ に適 用 し、物理 量 の観測 の有無 の効 果 、組み 込 むべ き フー リエ級 数 の

項数 お よび、観測 情報 の入 手時 間間隔 の 影響 と濃度 予測 精度 との 関係 な どにつ いて検討 を行 い 、両者 の フィ

ル ターの特 性 につ いて比較 したも ので あ る。

2. 濃 度予測 手法 の定式 化

2. 1. 2次 フィル タ ーの算定5)∼7)

2次 フ ィル ターは式(1)で 示 され る連続 時 間非 線形状 態方 程式 に した がって 遷移す るシステ ム状態 量Xを 推

定 す るの に適 用 され る理 論 で、本 研究 の濃度 予測 の問 題 では 、Xは線 形 の観測 方程式(2)を 通 して あ る離 散時

間間隔 ご とに観測 される 。

(2)

ここ に、 t:時 間 を表す 変数

k:観 測量 の得 られ る時点

X:シ ステム の状態 量ベ ク トル

f:ベ ク トルXの 非線 形 ベ ク トル関 数

y:観 測量 ベ ク トル

H:既 知 の観測 行列

v:シ ステ ム雑 音 ベ ク トルで 、平均 値0、 共 分散行 列Vの 正 規性 白色雑 音

w:観 測雑 音 ベ ク トル で、平 均値0、 共 分散行 列Wの 正規 性 白色雑 音

次 に あ る 時 刻t(k<t≦k+1)に お い て 、t=kま で の 観 測 量 を用 い て 得 られ る シ ス テ ム 状 態 量Xの 最 適 推 定 値 お よび そ の 推 定 誤 差 共 分 散 行 列 をそ れ ぞ れ 父(t│k)、P(t│k)と 表 示 す る と、 これ らは 次 の 微 分 方 程 式 を解 く こ と に よ っ て 求 め られ る 。

(3)

(4)

こ こ に 、 μ:bias correctionで 、

(5)

Φi:i番 目 の 要 素 が1、 そ の他 の要 素 が0で あ るベ ク トル tr(・):行 列 の ト レ ー ス Bi:ベ ク トル 関 数fのHessian行 列 A:ベ ク トル 関 数fのJacobian行 列T:転 置 記 号 観 測 量yがt=k+1で 得 ら れ る と す れ ば 、X(k+1│k+1)お よ びP(k+1│k+1)は カル マ ン フ ィ ル タ ー理 論 に よ り次 式 で 求 め られ る 。

(6)

(7)

こ こ に 、 I:単 位 行 列K:カ ル マ ン ・ゲ イ ン ・マ トリ ッ ク ス で 、

(8)

な お 、1次 フ ィル タ ー(拡 張 カル マ ン フ ィ ル タ ー)の 場 合 は μ(t)=0で あ る 。

2. 2. 2次 フ ィル ターに よる濃度 予測 計算 の定式 化

本論 文で は次式 に示 す よ うな定 係数 一 次元確 率移流 分 散方 程式 を考 え る ことにす る。

(9)

こ こ に 、 C:濃 度(g/m3)u:流 速(m/day) D:分 散 係 数(m2/day) γ:一 次 反 応 係 数(1/day)

x:距 離(m)t=時 刻(day) ε:平 均 値0,分 散 σε2の正 規 性 白色 雑 音(g/(m3day))

なお物理 量u,D,γ

は時空 間的 に一定 の場 合 を想定 して い る。

さて、式(9)に は右辺 第3項 に確率 項 が含 まれて お り、 また初期条 件 、境界 条件 お よび物理 量 が未知 の場

合 に は、式(9)を 解 析 的 に解 くこ とは不可 能で あ る。 また通常 の差分 法 や有限 要素 法 をで は、移流 項 が分散

項 よ り卓 越す ると きには移流 項 の離散 化誤 差 が大 き く影 響 し、不確 実性 が増 大す る こ とが予想 され る8)。 し

たが って ここで は、空 間的 に任意 に設 置 され た観 測点 よ り得 られ る情報 を基 に物 理量 を推定 しつ つ濃 度Cを

予測 す る こと を考 え る。式(9)を 常微 分方程式 に直す た めに濃 度C(x,t)お よび 雑音 項 ε(x,t)を次 式 のよ うに

変数Xに 関 して フー リエ級数 に展 開 し、任意 の波 数 の フー リエ係 数 につ いて解 析 を行 う。

(10)

(11)

こ こ に 、 B0:濃 度Cの 平 均 値 に相 当 す る フー リエ 係 数(g/m3)M:取 り込 む フ ー リ エ級 数 の 項 数 Am,Bm:濃 度Cの 波 数mに 関 す る フー リエ 係 数(g/m3)2:基 本 波 長 Em,Fm:ε の波 数mに 関 す る フ ー リエ 係 数(g/m3) 式(10),(11)を 式(9)に 代 入 す る と、B0(t)お よ び 波 数mの フ ー リエ 係 数 に 関 して 次式 の よ う な 連 立 常 微 分 方 程 式 が得 られ る 。

(3)

(12)

こ こ に 、

(13)

(14)

こ こで は 、 カ ル マ ン フ ィ ル タ ー に お け る シ ス テ ム 状 態 量 と して 式(9)の 物 理 量u,D,γ お よ び 式(10)の フ ー リエ 係 数B0,A m,Bm(m=1∼M)を ま と め てXと お く。 す な わ ち 、

(15)

こ こで 、 式(9)のu,D,γ の 値 は 時 空 間 的 に一 定 の 場 合 を考 え て い る の で 、 こ れ ら に 対 す る遷 移 方 程 式 は 次 式 と な る 。

(16)

よ っ てX(t)の シ ス テ ム方 程 式 は 式(12)と 式(16)を ま と め て 次 式 と な る 。

(17)

次 に式(2)に お け る観 測 量yと して 、x軸 上 に 任 意 に 配 置 さ れ た 観 測 点Xn(n=1∼N)か ら一 定 の 時 間 間 隔 で 得 ら れ る 濃 度Cを 考 え る 。 こ の 場 合 観 測 雑 音wが 混 入 す る も の とす る 。 した が っ て 式(2)に お け るyの あ る観 測 点xに 対 応 す る成 分yxは 次 式 で 表 され る 。

(18)

式(18)の 右 辺 第1項 の ベ ク トル をす べ て の 観 測 地 点 に つ い て ま と め て マ ト リ ッ ク ス 表 示 す る と 、式(2)の 観 測 行 列H(k)[{N×(2M+4)}行 列]は 次 式 の よ う に な る 。

(19)

な お 上 記 の観 測 点 以 外 の 点 でu,D,γ が観 測 さ れ る場 合 に は 、 式(19)の 行 数 を増 やせ ば よ い 。 以 上 の よ う に して カ ル マ ン フ ィ ル タ ー の シ ス テ ム方 程 式(1)と 観 測 方 程 式(2)が 定 式 化 さ れ た の で 、 あ と は 式(11)のl、M、 濃 度 観 測 点 数Nお よ び そ の 地 点xnを 定 め 、 状 態 量Xの 初 期 推 定 値X(0│0)お よ び そ の 推 定 誤 差 共 分 散 行 列 の初 期 値P(0│0)、 シ ス テ ム 雑 音,観 測 雑 音 の 分 散 共 分 散 行 列V(k),W(k)を 適 当 に 与 え れ ば 、2. 1.の2次 フ ィ ル タ ー に よ り、 式(3)、(4)の 離 散 時 間 間 隔 Δtご と にXの 最 適 推 定 値 が 求 め られ 、 これ を用 い て 濃 度Cの 予 測 が行 わ れ る 。 3. 適 用 例 こ こ で は 、 本 手 法 の有 用 性 お よ び 特 性 の 把 握 の た め に 、 以 下 の 方 法 で 発 生 さ せ た 真 値 の 分 か っ て い る デ ー タ に対 し て検 討 す る 。 確 率 項 の加 わ っ た 移 流 分 散 方 程 式 に 従 う濃 度 を模 擬 発 生 させ る場 合 、 ま ず 図-1に 示 し た初 期 濃 度 分 布 よ り フ ー リ エ 係 数 の 初 期 値 を求 め る 。 次 い で 、u=1.0(m/day),D=0.5(m2/day),γ=0. 01(1/day),Δt=0.1(day),基 本 波 長l=100(m),フ ー リエ 項 数M=20と し、 式(14)の 雑 音 項Em(t),Fm(t)に は平 均0、 分 散(0.005)2の 正 規 乱 数 を与 え 、 ル ンゲ ・ク ッ タ ー ・ギル(Runge-Kutta-Gill)法 に よ り式(12)を 数

(4)

擬 発 生 す る 。 さ ら に模 擬 発 生 し て 得 られ た フ ー リエ 級 数 か ら式 (10)に よ り濃 度 分 布 を計 算 し、 これ を 真 値 と し た4)。 こ こで 図 -1の 初 期 濃 度 分 布 は 式(9)の 方 程 式 で 雑 音 項 ε(x,t)を 除 い た 確 定 方 程 式 に お い て 、C(X,0) =1(0≦x≦5),C(x,0)=0(-∞ <x<0,5<x< ∞)の 初 期 条 件 で 解 析 解 にt=10(day)を 与 え た も の で あ る 。以 上 の よ う に して 式 (9)の 定 係 数 確 率 移 流 分 散 方 程 式 に 従 う濃 度 分 布 の 真 値 を模 擬 発 生 させ た結 果 を図-1の 実 線 に示 して い る 。 な お 、 実 際 に 各 地 点 で 観 測 さ れ る 濃 度 と して 上 述 の真 値 に 平 均0, 分 散(0.01)2の 観 測 雑 音 を加 え た 。 ま た 物 理 量 が 任 意 の1点 で 観 測 さ れ る 場 合 に は そ の観 測 物 理 量 に は 真 値 の10%の2乗 の 分 散 を も つ 雑 音 を 加 え て い る 。 次 に模 擬 発 生 さ せ た濃 度 分 布 に 対 して 2.で 述 べ た 手 法 に よ り濃 度 予 測 を行 う。 こ の 場 合 物 理 量u,D,γ お よ び フ ー リエ 係 数Bo,Am,Bmの 初 期 値 す な わ ち 状 態 量X(0 │0)に は 真 値 の50%の 値 を与 え た 。 こ こ で 濃 度 観 測 点 が 図-1の 矢 印 に 示 した よ う に ラ ン ダ ム に20地 点 、物 理 量u,D,γ の 観 測 は な い 場 合 で 、 観 測 時 間 間 隔 を10ス テ ッ プ す な わ ち10Δt=1(day)、 フ ー リ エ 級 数 の 項 数M=20と し た 場 合 の2次 フ ィ ル タ ー に よ る10ス テ ッ プ 先 の 予 測 を 図-1に 示 して い る。 ま た物 理 量u,D,γ フ ー リエ 係 数Boお よ びAm,Bmの 同 定 結 果 の 一 部 を図-2に 示 して い る 。 次 に、 図-1 の条 件 に お い て 物 理 量 が 観 測 され る場 合 、 観 測 時 間 間 隔 が 変 化 す る 場 合 、 お よび フ ー リエ 級 数 の項 数 が 変 化 す る 場 合 の1次 お よ び2次 フ ィ ル タ ー の 濃 度 予 測 精 度J の変 化 を そ れ ぞ れ図-3,図-4,図-5 に示 し て い る 。 こ の 場 合 、Jと して 次 の よ う な 指 標 を用 い た 。

図-

1 模 擬 発 生 濃 度(真 値)と そ の予 測 値

図-2 物 理 量 お よび フ ー リエ 係 数 の 同定 結 果

(5)

(20)

(21)

こ こ に 、C(x,k):濃 度 の 真 値 C(x,k│k-1):t=k-1ま で の 観 測 量 が得 られ た 場 合 の 濃 度C(x,k)の 予 測 値 N1:評 価 地 点 数 で 、 こ の場 合0mか ら100mま で1m お き に101地 点 を取 っ た 。 N2:評 価 時 点 数 で 、 こ の 場 合100時 点 か ら500時 点 まで の100時 点 お き の5時 点 を と っ た 。

4. 考察

図-1の

実 線 よ り時間 の経過 と共 に濃 度分布 の乱 れ が大

き くなっ てい るの が分 かる。 また2次 フィル ターによ る10

ステ ップ先の予 測値 をみ ると、雑音 に埋 もれ て い る高周 波

成分 の小 さな変動分 布 の予測 は困難 な もの の、時 点の進 行

と ともに濃度 分布 が乱 れて も精度 良 く濃 度分 布形 状 を予測

して いる といえ よ う。

次 に図-2のu,D,γ

の同定結 果 をみ る と、お よそ150

ステ ップ くらいで同定 値 が収束 して お り、いず れ も精 度 良

く同定 されて いる といえ る。 また フー リエ係 数B。,A1に つ

いて は同定 値 は100ス テ ップ まで には収束 し、 そ れ以後真

値 の変動 を精度 良 く同定 して い る。A10に つ いて は高周波

成分 の ため、 シス テム雑 音 の影響 を強 く受 け真値 お よび 同

定値 とも かな り変動 して お り、 同定精度 は低 くなって い る。

これ は高 波数 の フ ー リエ係数 が時 点 の進行 に伴 い雑音 レベ

ルの 中に埋 もれて しま うため と考 え られ る4)。

次 に図-3よ

り、1次 、2次 フィル ター とも濃度予 測 精

度 に差 は な く、濃 度観 測点 数Nが10以 下 の場合 には、物理

量 を観 測 した方 が予 測精度 が向上 して い る。 また物 理量 の

うちuを観 測 した場合 がもっ とも精度 が良 く、つ ついてD,

γの順 とな ってい る。 またNを15点 以上 に取 ると 、物理

量 の観測 の有 無 にか かわ らずJの 値は0.008程 度 に落 ち着

いて いる 。以 上 よ り観測 点数 を多 く取 る こ と、お よび 観

測 点数 が十分 で な い とき には流速uを 観測す る ことが予

測精度 向上 につ な がる と考 え られ る。

次 に図-4よ

り、1次 、2次 フィル ター ともNを15以

上 とる と予 測精 度 は安定 して くるが、 その精 度は観 測 時

間間 隔 が短 いほ ど高 くなって い る。 これ は観 測時 間間 隔

が短 い程 、観測 情報 が多 く入 手 され る ことお よび短 期 間

予測 とな るためで あ る と考 え られ る。 ま たNが10以 下 の

場 合 は観 測時 間間 隔 が小 さ くな って も必ず しも精 度 向上

して いな い。 これ は 、カル マ ン フィル ターは濃 度観測 地

図-3 物 理量 の観 測 の有 無 に 対 す る予 測 精 度 図-4 観 測 時 間 間 隔 の変 化 に 対す る予 測 精 度 図-5 フ ー リエ級 数 の項 数 変 化 に対 す る予 測 精 度

(6)

全体 的な分 布 は逆 に合わ な くな るため と考え られ る。事 実、 濃度観 測 地点 のみ で予測 精度 を評価 す る と観 測

時間 間隔 の短 い方 が高精 度 となっ てい る。

次 に図-5よ

り、1次 、2次 フィル タ ーともNを15以 上 とれ ば フ ー リエ級数 の項 数10の 濃 度予 測精度 は 項

数20の そ れ とほ とん ど同程度 にな って い る。 またNが10以 下 と少 な くな ると項 数10の 方 が項 数20よ りも濃 度

予測精 度 が良 くな って い る。 これ は層が少 ない場合 、項数 が多 くな る と同定 すべ きパ ラメ-タ が多 くな るの

で、 これ らのパ ラ メー タ をかえっ て精 度 良 く同定 で きず 濃度 予測精 度 が低 くなるも の と考 え られ る。

最後 に図 −1と 同 じ条 件で1次 フィル ター を用 いて濃 度予 測 を行っ た場合 、 お よそ半分 の演算 時間 で

図-1と ほぼ 同様 の濃 度 予測結 果 が得 られ た。 これは状 態方 程式 の非線 形 ベ ク トル関 数fに そ れほ ど強 い非 線形

性 がな く、 したが って本 手法 に関す る限 り、2次 フィル ターのbias correctionは 演算 時 間 を要す るだ けで

予測精 度 の向上 に有 効で はな かっ た と考 え られ る。

5. むす び

本報 で は、空間 的 に任意 に設置 され た観測点 か ら、 あ る時 間間 隔で 入手 され るデ ー タを用 い て、2次 フィ

ル ターに よ り定係 数 一次元 確率移 流分 散 方程式 で記述 され る濃度分 布 を予測 す る手法 を定式 化 し、 次いで 本

手 法 を模擬発 生濃 度 分布 に適用 した。 その結果 、本手 法 によ り精度 良 く濃度 分布 を予 測す る こと がで きる こ

とが示 され 、 また物理 量 の観測 の有無 、観 測 時間間 隔 お よび 読み込 む べ き フー リエ級数 の項 数 な どの濃度 分

布 予測 精度 に対す る影響 が評価 され た 。また 、2次 フィル タ ーの予 測精 度 は1次 フィル ターの それ と比較 し

てほ とん ど同様で あ る が、演算 時間 は およそ2倍 を必 要 としてお り、実 際上 は1次 フィル ター で十分 高い 予

測 が可 能で あ るこ とが示 され た。

参考 文 献

1) 神 野健二 ・上 田年 比古 ・安 田

裕: 実 時 間観測 デ ー タに よる定 係数2次 元移 流分 散方程 式 のパ ラ メー タ

同定 と濃 度予 測 、第30回 水理 講演 会論 文集 、pp. 313∼318, 昭 和61年2月.

2) 河 村

明 ・神 野健 二 ・上 田年 比古: 定 係数 一 次元確 率移流 分散 方程 式 の濃度 予測手 法 につ いて、 土木 学

会 第42回 年次 学術 講演 会概要 集第2部 、pp. 6∼7, 昭 和62年9月.

3) 吉 永宙 司 ・河 村

明 ・神野健 二 ・上 田年比古: 定 係 数一 次元確 率移 流分 散方 程式 の濃度 予測 手法 の特 性

につ いて(第2報)、 土木 学会 第43回 年 次学術 講演 会 概要集 第2部 、pp. 194∼195, 昭和63年10月.

4) 河村

明 ・神 野健 二 ・上 田年 比古 ・吉 永宙 司: カル マ ン フィル ター によ る定 係数一 次元確 率 移流 分散 方

程式 のオ ンライ ン濃度 分布予 測 につ いて、 九州大 学 工学集 報 、第62巻 、 第1号 、pp. 17∼24, 1989年1月.

5) Athans, M., Wisher, R. P. and Bertolini, A.: Suboptimal state estimation for

continuous-time nonlinear systems from discrete noisy measurements. IEEE Trans. on Automatic Control,

Vol. AC-13, No. 5, pp. 504•`514, October, 1968.

6) Geb, A. (ed.) : Applied Optimal Estimation, The MIT Press, pp.180•`228, 1974.

7) 神 田 徹 ・藤 田 睦 博: 新 体 系 土 木 工 学26水 文 学 、 技 報 堂 出 版 、pp. 231∼246, 1982年1月.

8) 神 野 健 二 ・上 田 年 比 古:粒 子 の 移 動 に よ る移 流 分 散 方 程 式 の 数 値 解 法 の 検 討 、 土 木 学 会 論 文 報 告 集 、 第 271号 、pp. 45∼46, 1978年3月.

参照

関連したドキュメント

For the multiparameter regular variation associated with the convergence of the Gaussian high risk scenarios we need the full symmetry group G , which includes the rotations around

The torsion free generalized connection is determined and its coefficients are obtained under condition that the metric structure is parallel or recurrent.. The Einstein-Yang

[56] , Block generalized locally Toeplitz sequences: topological construction, spectral distribution results, and star-algebra structure, in Structured Matrices in Numerical

By interpreting the Hilbert series with respect to a multipartition degree of certain (diagonal) invariant and coinvariant algebras in terms of (descents of) tableaux and

Kilbas; Conditions of the existence of a classical solution of a Cauchy type problem for the diffusion equation with the Riemann-Liouville partial derivative, Differential Equations,

Maria Cecilia Zanardi, São Paulo State University (UNESP), Guaratinguetá, 12516-410 São Paulo,

Debreu’s Theorem ([1]) says that every n-component additive conjoint structure can be embedded into (( R ) n i=1 ,. In the introdution, the differences between the analytical and

Then it follows immediately from a suitable version of “Hensel’s Lemma” [cf., e.g., the argument of [4], Lemma 2.1] that S may be obtained, as the notation suggests, as the m A