概要

　地すべりの発生において、地下水の影響を知ることはとても重要です。この事について簡単に説明すれば、地下水の影響により浮力が生じ、地盤の垂直応力が弱まることで、せん断応力が（相対的に）大きくなり、滑りやすくなります。（塑性変形が起こりやすくなると言っていいのだろうか？）

　地下水流動解析（浸透流解析）を行う際、透水係数 $K$ などのパラメーターを設定する必要があります。また、井戸揚水や注水、降雨による浸透により外部から地盤に水が流入または流出します。解析において、その流入出量 $q$ を設定する必要があります。

　一般的には、地下の構造は不明であり、透水係数 $K$ や流入出量 $q$ を上手く設定できないことがあり、浸透流解析の結果と実際の観測結果と一致しない場合があります。

　この記事では、アジョイント法を用いて、浸透係数 $K$ と流入出量 $q$ のパラメータ推定に関することをメモ程度に書き込んでいきます。

　アジョイント法については、4次元変分法の記事で説明しました。

最尤推定の復習

　今回の記事では、4次元変分法の技術であるアジョイント法が使われる。4次元変分法で求めたい値は、シミュレーションの最適な初期値 $x_0$ である。観測値が与えられたとき、確率が最も高くなる初期値を求めることを考える。確率分布は、ガウス分布に従うと仮定する。

　シミュレーションで解きたい方程式は、

\begin{align}
x_t = M(x_{t-1})
\end{align}

であるとする。$t$ は時間を表す。

　初期値 $x_0$ が従う確率分布は、

\begin{align}
P(x_0)\propto\exp\left\{- \frac{1}{2} \left(x_0-x^b \right)^T B^{-1} \left(x_0-x^b \right) \right\}
\end{align}

とする。$B$ は背景誤差共分散行列であり、$x_b$ は、背景場（事前情報）である。指数関数の前の係数は、今後の計算において不要だったので省いた。

　初期値 $x_0$ が与えられたとき、観測値 $y$ が従う確率分布は、

\begin{align}
P(y^{obs}_0,y^{obs}_1,\cdots,y^{obs}_T|x_0)	\propto \exp\left\{- \frac{1}{2} \sum_{t=0}^T \left(H(x_t)-y^{obs}_t \right)^T R_t^{-1} \left(H(x_t)-y^{obs}_t \right)\right\}
\end{align}

とする。$R_t$ は観測誤差共分散行列であり、$H_t$ は観測演算子である。（観測演算子は数値計算の結果を代入することと考えてよい。）

　ここで求めたい確率分布は、観測値 $y$ が与えられたとき、初期値 $x_0$ が従う確率分布である。ベイズの定理を用いれば

\begin{align}
P(x_0|y^{obs}_0,y^{obs}_1,\cdots,y^{obs}_T)　= \frac{P(y^{obs}_0,y^{obs}_1,\cdots,y^{obs}_T|x_0)P(x_0)	}{P(y^{obs}_0,y^{obs}_1,\cdots,y^{obs}_T)}
\end{align}

と書けるので、最終的に、

\begin{align}
\newcommand{\argmax}{\mathop{\rm arg~max}\limits}
x_0 &= \argmax_{x_0} \left\{P(x_0|y^{obs}_0,y^{obs}_1,\cdots,y^{obs}_T) \right\}
\\
&=  \argmax_{x_0}\left\{ P(y^{obs}_0,y^{obs}_1,\cdots,y^{obs}_T|x_0)P(x_0)\right\}
\end{align}

を解けば、求めたい初期値 $x_0$ が求まる。2行目は、ベイズの定理から得られた確率分布の分母の値（観測値が従う確率分布）は、初期値 $x_0$ に依存しないため省いた。

　上式をそのまま解くのは難しいため、以下の対数尤度関数を導入する。

\begin{align}
J(x_0)&\equiv -\log \left\{P(y^{obs}_0,y^{obs}_1,\cdots,y^{obs}_T|x_0)P(x_0)\right\}\\
&= \frac{1}{2} \left(x_0-x^b \right)^T B^{-1} \left(x_0-x^b \right) +
\frac{1}{2} \sum_{t=0}^T \left(H(x_t)-y^{obs}_t \right)^T R_t^{-1} \left(H(x_t)-y^{obs}_t \right)+\cdots
\end{align}

$\cdots$ は、$x_0$ に依存しない、対数尤度関数の最小化に寄与しない項である。対数尤度関数は、一般的には誤差関数（目的関数）と呼ばれる。求めたい初期値は、誤差関数が最小となる $x_0$ である。対数尤度関数の最小化は、準ニュートン法などを使い求める。

　この記事では、初期値 $x_0$ を求めることを目的としておらず、物性値であるパラメータをアジョイント法で求めることを目的とする。

飽和地下水流れの基礎方程式

　飽和地下水流れの基礎方程式は、質量保存則から導かれる。今回の記事では、以下導出される方程式を用いる。

　微小な領域（微小立方体）を考える。微小立方体に流体が流入・流出する。このとき、微小立方体における流体の質量保存則は、

\begin{align}
（時刻t+\Delta tにおける&流体質量）=（時刻tにおける流体質量）\\
&+（微小立方体の各面を通じた流入出質量）-（人為的に取り去られる質量）
\\
\therefore \ \ M_{t+\Delta t} &=M_t +\Delta m_{x}+ \Delta m_y+ \Delta m_z - m_q
\end{align}

となる。時刻$t$ における流体質量 $M_t$ は、間隙率を $\varphi$ （間隙体積➗全体体積）、密度を $\rho$ とすると

\begin{align}
M_{t} = \varphi(t,x,y,z)\rho(t,x,y,z) \Delta x \Delta y \Delta z
\end{align}

と書ける。

微小立方体の $x$ 方向の面を通じた流入質量は、

\begin{align}
（x 方向の面を通じた流入質量）&= (密度) \times （面積） \times （x 方向の流入速度）\times\Delta t \\

&=\rho(t,x,y,z)\times\Delta y \Delta z \times V_x(t,x,y,z)\times \Delta t
\end{align}

流出質量は、

\begin{align}
（x 方向の面を通じた流出質量）&= (密度) \times （面積） \times （x 方向の流出速度）\times \Delta t　\\

&=\rho(t,x+\Delta x,y,z)\times\Delta y \Delta z \times V_x(t,x+\Delta x,y,z)\times \Delta t
\end{align}

$y,z$ も同様である。

　微小立方体の各面を通じた流入出質量 $\Delta m_{x}, \Delta m_y, \Delta m_z$ は、各方向の流速を $V_x,V_y,V_z$ とすれば、

\begin{align}
\Delta m_{x} &= \left(\rho(t,x,y,z)V_x(t,x,y,z)-\rho(t,x+\Delta x ,y,z)V_x(t,x+\Delta x ,y,z) \right)\Delta t \Delta y \Delta z \\

\Delta m_{y} &= \left(\rho(t,x,y,z)V_y(t,x,y,z)-\rho(t,x,y+\Delta y,z)V_y(t,x ,y+\Delta y,z) \right)\Delta t \Delta x \Delta z \\

\Delta m_{z} &= \left(\rho(t,x,y,z)V_z(t,x,y,z)-\rho(t,x,y,z+\Delta z)V_z(t,x ,y,z+\Delta z) \right)\Delta t \Delta x \Delta y \\
\end{align}

と書ける。

　人為的に取り去られる質量 $m_q$ は、

\begin{align}
m_q = \rho_0 q \Delta x \Delta y \Delta z \Delta t
\end{align}

と書ける。$\rho_0$ は基準密度、$q\Delta x \Delta y \Delta z$ は単位時間あたりの揚水量である。

　まとめると、

\begin{align}
&\varphi(t+\Delta t,x,y,z)\rho(t+\Delta t,x,y,z) \Delta x \Delta y \Delta z
-\varphi(t,x,y,z)\rho(t,x,y,z) \Delta x \Delta y \Delta z \\
&=\left(\rho(t,x,y,z)V_x(t,x,y,z)-\rho(t,x+\Delta x ,y,z)V_x(t,x+\Delta x ,y,z) \right)\Delta t \Delta y \Delta z \\

&+\left(\rho(t,x,y,z)V_y(t,x,y,z)-\rho(t,x,y+\Delta y,z)V_y(t,x ,y+\Delta y,z) \right)\Delta t \Delta x \Delta z \\

&+\left(\rho(t,x,y,z)V_z(t,x,y,z)-\rho(t,x,y,z+\Delta z)V_z(t,x ,y,z+\Delta z) \right)\Delta t \Delta x \Delta y \\

&-\rho_0 q \Delta x \Delta y \Delta z \Delta t

\end{align}

となり、両辺を $\Delta t \Delta x \Delta y \Delta z$ で割り

\begin{align}
&\frac{\varphi(t+\Delta t,x,y,z)\rho(t+\Delta t,x,y,z)
-\varphi(t,x,y,z)\rho(t,x,y,z) }{\Delta t}\\

&=\frac{\rho(t,x,y,z)V_x(t,x,y,z)-\rho(t,x+\Delta x ,y,z)V_x(t,x+\Delta x ,y,z) }{\Delta x} \\

&+\frac{\rho(t,x,y,z)V_y(t,x,y,z)-\rho(t,x,y+\Delta y,z)V_y(t,x ,y+\Delta y,z) }{\Delta y} \\

&+\frac{\rho(t,x,y,z)V_z(t,x,y,z)-\rho(t,x,y,z+\Delta z)V_z(t,x ,y,z+\Delta z) }{\Delta z}\\

&-\rho_0 q 

\end{align}

極限 $\Delta t,\Delta x,\Delta y,\Delta z \to 0$ を実行すれば、

\begin{align}
\frac{\partial  }{\partial t}\left\{\varphi\rho \right\}

=&-\frac{\partial  }{\partial x}\left\{\rho V_x \right\}
-\frac{\partial  }{\partial y}\left\{\rho V_y \right\} 
-\frac{\partial  }{\partial z}\left\{\rho V_z \right\} 
-\rho_0 q 
\end{align}

が得られる。

　間隙率と密度の微小線形変化は、$h$ を水頭として、

\begin{align}
\rho = \rho_0\left\{1+c_{fh}(h-h_0) \right\} \ \ , \ \  \varphi = \varphi_0\left\{1+c_{rh}(h-h_0) \right\}
\end{align}

と書ける。$c_{fh}=\rho gc_f$ であり、$c_{fh}$ は水頭が変化した際の流体の密度変化率、 $c_f$ は水の圧縮率である。また、$c_{rh}$ は地層間隙の変化率である。これらの係数は、微小量とする。よって、時間微分の項は、

\begin{align}
\frac{\partial  }{\partial t}\left\{\varphi\rho \right\}
&= \frac{\partial  }{\partial h}\left\{\varphi\rho \right\}\frac{\partial h}{\partial t} \\

&= \left\{\frac{\partial \varphi }{\partial h}\rho +\varphi \frac{\partial \rho }{\partial h} \right\}\frac{\partial h}{\partial t} \\

&\approx \varphi_0\rho_0(c_{fh}+c_{rh})\frac{\partial h}{\partial t}\\

&＝ S_s \rho_0\frac{\partial h}{\partial t}
\end{align}

と表せる。$S_s =\varphi_0(c_{fh}+c_{rh}) $ は比貯留係数である。

　さらに、ダルシー流速

\begin{align}
V_x = -K_x \frac{\partial h }{\partial x}
 \ \ , \ \
V_y = -K_y \frac{\partial h }{\partial y}
 \ \ , \ \
V_z = -K_z \frac{\partial h }{\partial z}
\end{align}

を使い、水頭 $h$ の１次まで取り出し、$\rho_0$ で割れば、浸透流の支配方程式が得られる。

\begin{align}
S_s \frac{\partial h}{\partial t}
=&\frac{\partial  }{\partial x}\left\{K_x \frac{\partial h }{\partial x}\right\}
+\frac{\partial  }{\partial y}\left\{K_y \frac{\partial h }{\partial y} \right\} 
+\frac{\partial  }{\partial z}\left\{K_z \frac{\partial h }{\partial z} \right\} 
-q 
\end{align}

ラグランジュ関数

　ラグランジュ関数を構成する。上記で説明した支配方程式を離散化すると、

\begin{align}
S_s \frac{h^{n}_{i,j,k}-h^{n-1}_{i,j,k} }{\Delta t} &= \frac{1}{\Delta x}\left\{ 
K \frac{\partial h }{\partial x} \Biggr|_{i+1/2}-K \frac{\partial h }{\partial x} \Biggr|_{i-1/2}\right\}
+\frac{1}{\Delta y}\left\{ 
K \frac{\partial h }{\partial y} \Biggr|_{j+1/2}-K \frac{\partial h }{\partial y} \Biggr|_{j-1/2}\right\}
\\
&+\frac{1}{\Delta z}\left\{ 
K \frac{\partial h }{\partial z} \Biggr|_{k+1/2}-K \frac{\partial h }{\partial z} \Biggr|_{k-1/2}\right\}
-q
\\
&=\frac{1}{\Delta x}\left\{ 
K_{i+1/2,j,k}\frac{h^{n}_{i+1,j,k}-h^{n}_{i,j,k}}{\Delta x}-K_{i-1/2,j,k} \frac{h^{n}_{i,j,k}-h^{n}_{i-1,j,k}}{\Delta x}  \right\}
\\
&+\frac{1}{\Delta y}\left\{ 
K_{i,j+1/2,k}\frac{h^{n}_{i,j+1,k}-h^{n}_{i,j,k}}{\Delta y}-K_{i,j-1/2,k} \frac{h^{n}_{i,j,k}-h^{n}_{i,j-1,k}}{\Delta y}  \right\}
\\
&+\frac{1}{\Delta z}\left\{ 
K_{i,j,k+1/2}\frac{h^{n}_{i,j,k+1}-h^{n}_{i,j,k}}{\Delta z}-K_{i,j,k-1/2} \frac{h^{n}_{i,j,k}-h^{n}_{i,j,k-1}}{\Delta z}  \right\}
-q
\end{align}

となる。時間微分の離散化は、陰解法を用いる。陽解法だと時間ステップを小さくしなければ発散してしまう。

　まとめて

\begin{align}
R^n_{i,j,k} &=S_s \Delta x \Delta y \Delta z  \left( h^{n}_{i,j,k}-h^{n-1}_{i,j,k} \right) +q \Delta t \\
&-\Delta t \Delta y \Delta z \left\{ 
K_{i+1/2,j,k}\frac{h^{n}_{i+1,j,k}-h^{n}_{i,j,k}}{\Delta x}-K_{i-1/2,j,k} \frac{h^{n}_{i,j,k}-h^{n}_{i-1,j,k}}{\Delta x}  \right\}
\\
&-\Delta t \Delta x \Delta z \left\{ 
K_{i,j+1/2,k}\frac{h^{n}_{i,j+1,k}-h^{n}_{i,j,k}}{\Delta y}-K_{i,j-1/2,k} \frac{h^{n}_{i,j,k}-h^{n}_{i,j-1,k}}{\Delta y}  \right\}
\\
&-\Delta t \Delta x \Delta y \left\{ 
K_{i,j,k+1/2}\frac{h^{n}_{i,j,k+1}-h^{n}_{i,j,k}}{\Delta z}-K_{i,j,k-1/2} \frac{h^{n}_{i,j,k}-h^{n}_{i,j,k-1}}{\Delta z}  \right\}
=0

\end{align}

と書く。浸透流解析を行う際は、$h^{n}$ について上記の方程式を解く（連立方程式を解く必要がある）。

　ラグランジュ関数を

\begin{align}
\mathcal{L} = E(h,K,q)+ \sum_{n}\sum_{i,j,k}\lambda_{i,j,k}^n R^n_{i,j,k} 
\end{align}

のように構成する。$\lambda$ はラグランジュ未定乗数である。

　$E(h,K,q)$ は、誤差関数（目的関数）であり、

\begin{align}
E(h,K,q)=\frac{1}{2}\sum_{n}\sum_{i,j,k}\frac{\left(h^{n}_{i,j,k}-{h_{obs}}^{n}_{i,j,k} \right)^2}{{\sigma_{i,j,k}^n}^2}
\end{align}

と書ける。 $h_{obs}$ は観測量、$\sigma^2 $ は $h$ に対する誤差分散である。ただし、実際の浸透流解析において、観測値はすべてのグリッド点で与えられることはないだろう。

　誤差関数は、適切な $K,q$ が決まれば観測値と数値解は一致するので、 $h$ は陽に表さず

\begin{align}
E(K,q)=E(h,K,q)
\end{align}

と書ける。
（ただし、ラグランジュ関数 $\mathcal{L}$ の中では、$E(h,K,q)$ として計算を行う。）

　支配方程式が成立するならば、

\begin{align}
\mathcal{L} - E(K,q)=0 
\end{align}

であり、さらに、

\begin{align}
\delta \left(\mathcal{L} - E(K,q)\right)=0 
\end{align}

が成立する。

　もちろん、支配方程式は、ラグランジュ関数を $\lambda$ で微分すれば得られる。

\begin{align}
&\frac{\delta }{\delta \lambda^{n}_{i,j,k}}\left(\mathcal{L} - E(K,q)\right) =\frac{\partial \mathcal{L}}{\partial \lambda^{n}_{i,j,k}} = 0
\\
\therefore \ \ &R^n_{i,j,k} = 0
\end{align}

　浸透流解析は、通常のシミュレーションと同様に初期条件・境界条件を与え、その他 $K,q$ などのパラメータを設定しシミュレーションを行う。

アジョイント方程式

　水頭 $h$ に関するアジョイント方程式を導出する。

　まず、

\begin{align}
\frac{\delta }{\delta h^{n}_{i,j,k}}\left(\mathcal{L} - E(K,q)\right) = \frac{\partial \mathcal{L}}{\partial h^{n}_{i,j,k}} = 0
\end{align}

が成立する。よって、

\begin{align}
\frac{\partial \mathcal{L}}{\partial h^{n}_{i,j,k}} &= \frac{h^{n}_{i,j,k}-{h_{obs}}^{n}_{i,j,k} }{{\sigma_{i,j,k}^n}^2}
+S_s \Delta x \Delta y \Delta z  \left(\lambda^{n}_{i,j,k}-\lambda^{n+1}_{i,j,k} \right)
\\
&-\Delta t \Delta y \Delta z \left\{ 
K_{i+1/2,j,k}\frac{\lambda^{n}_{i+1,j,k}-\lambda^{n}_{i,j,k}}{\Delta x}-K_{i-1/2,j,k} \frac{\lambda^{n}_{i,j,k}-\lambda^{n}_{i-1,j,k}}{\Delta x}  \right\}
\\
&-\Delta t \Delta x \Delta z \left\{ 
K_{i,j+1/2,k}\frac{\lambda^{n}_{i,j+1,k}-\lambda^{n}_{i,j,k}}{\Delta y}-K_{i,j-1/2,k} \frac{\lambda^{n}_{i,j,k}-\lambda^{n}_{i,j-1,k}}{\Delta y}  \right\}
\\
&-\Delta t \Delta x \Delta y \left\{ 
K_{i,j,k+1/2}\frac{\lambda^{n}_{i,j,k+1}-\lambda^{n}_{i,j,k}}{\Delta z}-K_{i,j,k-1/2} \frac{\lambda^{n}_{i,j,k}-\lambda^{n}_{i,j,k-1}}{\Delta z}  \right\}
=0
\end{align}

となり, $\lambda^{n}$ について上記の方程式を解けば（連立方程式を解く必要がある）、アジョイント方程式が得られる。

　支配方程式 $R_{i,j,k}^n$ と比べると、右辺の一項目を除いて、アジョイントになっていることがわかる。つまり、 $h^n$ を $\lambda^n$ に入れ替え、時間ステップ $n$ を逆にした（行列表示したなら転置行列にした）方程式になっている。

　アジョイント方程式は、浸透流解析における最終時間ステップを $n=T$ とすれば、アジョイント方程式の初期値を

\begin{align}
\lambda^{T+1}_{i,j,k} = 0
\end{align}

として、時間方向が逆になるように

\begin{align*}
\lambda^{T+1}_{i,j,k} = 0 \Rightarrow \lambda^{T}_{i,j,k} \Rightarrow  \lambda^{T-1}_{i,j,k}\Rightarrow \cdots \Rightarrow \lambda^{0}_{i,j,k}
\end{align*}

と計算する。

　通常のデータ同化における4次元変分法では、初期値 $h^0$ を求めたい場合が多く、 $n=0$ における勾配 $g$ は

\begin{align}
g_{i,j,k} = \lambda^{0}_{i,j,k}
\end{align}

である。(誤差関数　$E(K,q)\to E(h^0,K,q)$　と変更する必要がある。）

　最急降下法を用いるならば

\begin{align}
h^{0}_{i,j,k} \leftarrow h^{0}_{i,j,k} - \alpha g_{i,j,k}
\end{align}

として初期値 $h^0$ を更新する。観測結果と一致する初期値 $h^0$ を反復的に求める。

　実際の計算では、準ニュートン法が使われる。

浸透係数・流入出量に関する勾配

　今回の記事で求めたい値は、初期値 $h^0$ ではなくパラメータ $K,q$ である。つまり、求めたい勾配 $g^K,g^q$ は、

\begin{align}
g^{K}_{i,j,k} &= \frac{\delta E}{\delta K_{i,j,k}} \\
g^{q} &= \frac{\delta E}{\delta q} \\
\end{align}

である。最急降下法を用いるならば

\begin{align}
K_{i,j,k} &\leftarrow K_{i,j,k} - \alpha g^{K}_{i,j,k} \\
q &\leftarrow q - \alpha g^{q} \\
\end{align}

としてパラメータ $K,q$ を更新する。

　パラメータ $K,q$ における勾配は、アジョイント方程式と同様に、ラグランジュ関数を微分すれば得られる。

浸透係数

　浸透係数 $K$ に関しては、

\begin{align}
\frac{\delta }{\delta K_{i,j,k}} \left(\mathcal{L} - E(K,q)\right)= 0 
\end{align}

が成立する。ここで、

\begin{align}
K_{i+1/2,j,k} &= \frac{1}{2}\left(K_{i+1,j,k}+K_{i,j,k} \right)
 \ \ , \ \
K_{i-1/2,j,k} = \frac{1}{2}\left(K_{i,j,k}+K_{i-1,j,k} \right)
\\
K_{i,j+1/2,k} &= \frac{1}{2}\left(K_{i,j+1,k}+K_{i,j,k} \right)
 \ \ , \ \
K_{i,j-1/2,k} = \frac{1}{2}\left(K_{i,j,k}+K_{i,j-1,k} \right)
\\
K_{i,j,k+1/2} &= \frac{1}{2}\left(K_{i,j,k+1}+K_{i,j,k} \right)
 \ \ , \ \
K_{i,j,k-1/2} = \frac{1}{2}\left(K_{i,j,k}+K_{i,j,k-1} \right)
\\
\end{align}

とすれば、

\begin{align}
\frac{\delta \mathcal{L}}{\delta K_{i,j,k}} &= \frac{\partial \mathcal{L}}{\partial K_{i+1/2,j,k}}\frac{\delta  K_{i+1/2,j,k}}{\delta K_{i,j,k}} 
+\frac{\partial \mathcal{L}}{\partial K_{i-1/2,j,k}}\frac{\delta  K_{i-1/2,j,k}}{\delta K_{i,j,k}} 
\\
&+\frac{\partial \mathcal{L}}{\partial K_{i,j+1/2,k}}\frac{\delta  K_{i,j+1/2,k}}{\delta K_{i,j,k}} 
+\frac{\partial \mathcal{L}}{\partial K_{i,j-1/2,k}}\frac{\delta  K_{i,j-1/2,k}}{\delta K_{i,j,k}} 
\\
&+\frac{\partial \mathcal{L}}{\partial K_{i,j,k+1/2}}\frac{\delta  K_{i,j,k+1/2}}{\delta K_{i,j,k}} 
+\frac{\partial \mathcal{L}}{\partial K_{i,j,k-1/2}}\frac{\delta  K_{i,j,k-1/2}}{\delta K_{i,j,k}} 
\\
\end{align}

\begin{align}
&=\frac{1}{2}\left(\frac{\partial \mathcal{L}}{\partial K_{i+1/2,j,k}}
+\frac{\partial \mathcal{L}}{\partial K_{i-1/2,j,k}}
+\frac{\partial \mathcal{L}}{\partial K_{i,j+1/2,k}}
+\frac{\partial \mathcal{L}}{\partial K_{i,j-1/2,k}} 
+\frac{\partial \mathcal{L}}{\partial K_{i,j,k+1/2}}
+\frac{\partial \mathcal{L}}{\partial K_{i,j,k-1/2}} \right)
\\
\end{align}

と書ける。

　ラグランジュ関数を $K_{i+1/2,j,k}$ で微分すると

\begin{align}
\frac{\partial \mathcal{L}}{\partial K_{i+1/2,j,k}}　&=\sum_n  \Delta t \Delta y \Delta z \left\{ 
\frac{h^{n}_{i+1,j,k}-h^{n}_{i,j,k}}{\Delta x}\lambda^{n}_{i,j,k}-\frac{h^{n}_{i+1,j,k}-h^{n}_{i,j,k}}{\Delta x}\lambda^{n}_{i+1,j,k} \right\}
\\
&=\sum_n \frac{\Delta t \Delta y \Delta z}{\Delta x} \left(h^{n}_{i+1,j,k}-h^{n}_{i,j,k} \right) \left(\lambda^{n}_{i+1,j,k} -\lambda^{n}_{i,j,k}  \right)
\end{align}

となる。他も同様に計算すれば、

\begin{align}
\frac{\partial \mathcal{L}}{\partial K_{i-1/2,j,k}}　&=
\sum_n \frac{\Delta t \Delta y \Delta z}{\Delta x} \left(h^{n}_{i,j,k}-h^{n}_{i-1,j,k} \right) \left(\lambda^{n}_{i,j,k} -\lambda^{n}_{i-1,j,k}  \right)
\\
\frac{\partial \mathcal{L}}{\partial K_{i,j+1/2,k}}　&=
\sum_n \frac{\Delta t \Delta x \Delta z}{\Delta y} \left(h^{n}_{i,j+1,k}-h^{n}_{i,j,k} \right) \left(\lambda^{n}_{i,j+1,k} -\lambda^{n}_{i,j,k}  \right)
\\
\frac{\partial \mathcal{L}}{\partial K_{i,j-1/2,k}}　&=\sum_n \frac{\Delta t \Delta x \Delta z}{\Delta y} \left(h^{n}_{i,j,k}-h^{n}_{i,j-1,k} \right) \left(\lambda^{n}_{i,j,k} -\lambda^{n}_{i,j-1,k}  \right)
\\
\frac{\partial \mathcal{L}}{\partial K_{i,j,k+1/2}}　&=\sum_n \frac{\Delta t \Delta x \Delta y}{\Delta z} \left(h^{n}_{i,j,k+1}-h^{n}_{i,j,k} \right) \left(\lambda^{n}_{i,j,k+1} -\lambda^{n}_{i,j,k}  \right)
\\
\frac{\partial \mathcal{L}}{\partial K_{i,j,k-1/2}}　&=\sum_n \frac{\Delta t \Delta x \Delta y}{\Delta z} \left(h^{n}_{i,j,k}-h^{n}_{i,j,k-1} \right) \left(\lambda^{n}_{i,j,k} -\lambda^{n}_{i,j,k-1}  \right)
\\
\end{align}

となる。最終的に $K$ の勾配は、アジョイント方程式から計算されたラグランジュ未定乗数を使い、

\begin{align}
\frac{\delta E}{\delta K_{i,j,k}}=
&\frac{1}{2}\sum_n \frac{\Delta t \Delta y \Delta z}{\Delta x} \left(h^{n}_{i+1,j,k}-h^{n}_{i,j,k} \right) \left(\lambda^{n}_{i+1,j,k} -\lambda^{n}_{i,j,k}  \right)
\\
&+\frac{1}{2}\sum_n \frac{\Delta t \Delta y \Delta z}{\Delta x} \left(h^{n}_{i,j,k}-h^{n}_{i-1,j,k} \right) \left(\lambda^{n}_{i,j,k} -\lambda^{n}_{i-1,j,k}  \right)
\\
&+\frac{1}{2}\sum_n \frac{\Delta t \Delta x \Delta z}{\Delta y} \left(h^{n}_{i,j+1,k}-h^{n}_{i,j,k} \right) \left(\lambda^{n}_{i,j+1,k} -\lambda^{n}_{i,j,k}  \right)
\\
&+\frac{1}{2}\sum_n \frac{\Delta t \Delta x \Delta z}{\Delta y} \left(h^{n}_{i,j,k}-h^{n}_{i,j-1,k} \right) \left(\lambda^{n}_{i,j,k} -\lambda^{n}_{i,j-1,k}  \right)

\\
&+\frac{1}{2}\sum_n \frac{\Delta t \Delta x \Delta y}{\Delta z} \left(h^{n}_{i,j,k+1}-h^{n}_{i,j,k} \right) \left(\lambda^{n}_{i,j,k+1} -\lambda^{n}_{i,j,k}  \right)

\\
&+\frac{1}{2}\sum_n \frac{\Delta t \Delta x \Delta y}{\Delta z} \left(h^{n}_{i,j,k}-h^{n}_{i,j,k-1} \right) \left(\lambda^{n}_{i,j,k} -\lambda^{n}_{i,j,k-1}  \right)

\end{align}

と求まる。

流入出量

　流入出量 $q$ も同様に、

\begin{align}
\frac{\delta }{\delta q} \left(\mathcal{L} - E(K,q)\right)= 0 
\end{align}

が成立する。ラグランジュ関数を $q$ で微分すると

\begin{align}
\frac{\delta \mathcal{L}}{\delta q}　&=\sum_n \sum_{i,j,k}\lambda^{n}_{i,j,k} \Delta t  
\end{align}

が成立するので、最終的に $q$ の勾配は、

\begin{align}
\frac{\delta E}{\delta q}=\sum_n\sum_{i,j,k} \lambda^{n}_{i,j,k} \Delta t 
\end{align}

と求まる。

　実際の $q$ の値は、時間や空間に依存するので、

\begin{align}
\frac{\delta E}{\delta q^n_{i,j,k}}= \lambda^{n}_{i,j,k} \Delta t 
\end{align}

を解くことになる。

その他

背景誤差項

　背景誤差項は、事前分布を仮定すれば得られる項で、機械学習における正則化の役割を果たす。浸透係数 $K$ において、事前分布における背景誤差分散の逆数を $\chi_K$ とすれば、背景誤差項は、

\begin{align}
E_b(K)=\frac{1}{2}\chi_K \sum_{i,j,k} \left(K_{i,j,k}-K^{b}_{i,j,k} \right)^2
\end{align}

と書ける。 $K^{b}$ は $K$ の事前情報であるが、設定することは難しい。しかし、背景誤差項を数値的な振動を取り除くために使用することができる。

　$K^{b}$ は平滑フィルタ（２次精度数値フィルタ）を各方向に用いて、

\begin{align}
\phi_{i,j,k} &= \frac{1}{4} \left(K_{i-1,j,k}+2K_{i,j,k}+K_{i+1,j,k}  \right)
\\
\varphi_{i,j,k} &= \frac{1}{4} \left(\phi_{i,j-1,k}+2\phi_{i,j,k}+\phi_{i,j+1,k}  \right)
\\
K^{b}_{i,j,k}&= \frac{1}{4} \left(\varphi_{i,j,k-1}+2\varphi_{i,j,k}+\varphi_{i,j,k+1}  \right)
\end{align}

とする。　流入出量 $q$ も同様である。

　ただし、浸透係数 $K$ は負になることはないので、単純な二乗誤差が使えるか疑問が残る。二乗誤差は、ガウス分布を対数尤度の形で表したものである。ガウス分布のとりうる値は、$-\infty$ から $\infty$ であり、浸透係数 $K$ が従う分布には適さないと考えられる。

浸透係数の上限・下限

　浸透係数 $K$ は負の値は取らないが、 $K$ の値を更新するときに負になる可能性がある。その問題を避けるため、勾配を計算する際に変数変換を行う。

　とあるパラメータ $\kappa$ を使い $K$ を

\begin{align}
K= \frac{K_{max}+K_{min}e^{-\kappa} }{1+e^{-\kappa}}
\end{align}

と表す。$K_{max}$ は $K$ の最大値であり、$K_{min}$ は $K$ の最小値である。$K_{max}=1,K_{min}=0$ ならば、深層学習で馴染みのあるシグモイド関数である。この関数を使えば、負の値をとりうることはない。

　勾配を計算するときは、 $K$ を $\kappa$ の関数に変換して

\begin{align}
\kappa = - \log\left(\frac{K_{max}-K}{K-K_{min}} \right)
\end{align}

勾配も

\begin{align}
\frac{\delta E}{\delta \kappa_{i,j,k}}=\frac{\partial K_{i,j,k}}{\partial  \kappa_{i,j,k}}\frac{\delta E}{\delta K_{i,j,k}}
\end{align}

\begin{align}
\frac{\partial K}{\partial \kappa}= \frac{e^{-\kappa} }{(1+e^{-\kappa})^2}(K_{max}-K_{min})
\end{align}

と変換して計算を行う。

参考文献

記事は、以下の書籍を参考にしています。

地下水のトレーサー試験 ―地下水の動きを知る―

地圏水循環の数理―流域水環境の解析法
登坂博行・著

データ同化流体科学
大林茂著・三坂孝志著・加藤博司著・菊地亮太著・照井伸彦編・小谷元子編・赤間陽二編・花輪公雄編