離散形河川水質モデルと不完全データ問題

青山 智夫, 神部 順子, 長嶋 雲兵


Return

1 はじめに

昔,川は上流から下流にかけて連続的に変化していく水の流れであったが,人口増加により用水を必要とし取水が行われるようになり,近年では下水道の整備により処理水が河川に流入するようになった.その量は支流量を超えることもある[1].下水処理水の水質は炭素に関する生物学的酸素要求量(C-BOD, Carbon-Biochemical Oxygen Demand)については次第に改善され,近年では支流とほぼ同じレベルである.しかし化学的酸素要求量(COD, Chemical Oxygen Demand),全窒素(T-N, Total-Nitrogen),全燐(T-P, Total-Phosphorus),N-BOD(Nitrogen-Biochemical Oxygen Demand)については1オーダ大きいこともある[2].このような環境に対する小さくない負荷が断続的に生じているのが現代の都市圏を流れる河川である.
河川の水質問題についてさまざまな研究が成されてきた.それらの研究の多くは膨大な観測を必要としたり,現象記述方程式(支配方程式)の知識を必要とする[3, 4].環境に対する大負荷であっても詳細に測定でき,支配方程式を決定できれば詳細な検討が可能である.しかし一般的に公開されている数km間隔の観測点間距離の河川データを用いて河川水質を論じるには既存の方法では困難である.都市部を流域とする河川では環境負荷が大きく水質状況を定量的に把握するには従来とは別の,観測数をそれ程要求せず,明示的な支配方程式を要求しない,新しい処理が求められている.我々は点と点の間の水質指標の離散的変化として河川水質を定式化し,河川水質の断続的変化を的確に表したい.
現象を離散的に定式化する方法の一つに階層型神経回路網(以下神経回路網という)がある.神経回路網が河川水量に適用された例[5]はある.水質変化について適用された例はないように思われる.
神経回路網は非線形関数適合(fittingという),補完(補間と補外の総称),分類機能を有し,目的現象の支配方程式が不明であっても,説明変数と現象の観測値から同方程式の模擬関数(emulation functionという)を自動的に生成する.この技術は薬理学では確立されている.自動生成された関数は離散データが基になっていても,連続かつ微分可能な関数である.その性質を利用して予測が可能である.これらの機能は河川水質問題を研究する上で有用ではないかと考えられる.ただし,その諸機能が発現するためには,説明変数と観測値に不足があってはならない.
不足とは観測データが本来あるべきことが自明であるのに何らかの理由で不明な場合と,現象を説明する観測自体がされていないか,そういう観測自体の存在不確実の三種類がある.前者を欠測という.欠測を含んだデータの情報処理は近年かなり研究されている[6].一方,後二者は前者以上に重要であると思われるが,ほとんど検討されていない.当然のことで,少ない観測から非線形現象を再現すること自体容易でないのに,さらにその再現の背景に不足データの存在を指摘するのは極めて難しい.しかし,我々は環境問題の真の理解のために後者は回避できない問題と考える.同時にそれは神経回路網出力の検定に関わる問題である.我々は神経回路網の検定を出力値の信頼性という観点から調査したが,そこには後者の概念が検討されていない.本論文の目的はこの観測データ不足問題を考察することである.
ほとんど全ての環境データに欠測が見出されるが,水質問題に関する諸数値はデータ不足の状況である.それにも関わらず,環境保護団体や行政府の議論にそれを考慮した例は少ない.定量性に疑問ある議論を一歩前進させるため,不足データを基に現象を解析する方法を検討したい.本論文では神経回路網を使用した非線形多変量解析のデータ不足状況における挙動を調査し,同回路網解析の適用限界を明らかにすることである.
観測データには測定誤差が含まれる.その統計的性質は確かでない.それを基に神経回路網の適用限界の議論は困難である.この測定誤差問題を回避するため河川モデルを考案し,それと乱数を基に人工河川を定義し,そのデータを基に神経回路網の適用限界を明らかにする方針である.

2 河川水質の離散形式

2. 1 定義

我々は都市河川の中下流域の河川浄化機能に焦点を絞り下記のモデルを考える.
【1】河川の性質は離散的な観測点列の連鎖で表せる.観測点間の水質変化を一指標の変化(ベクトル形式)として表す.
【2】観測点への流入は本流と「支流」の二要素とする.ここで支流とは本流への全流入を加算したものである.
我々は人工河川を取り扱うので,観測点間距離,支流の流入量は区間指定した一様乱数で作成する.支流水に含まれる物質量も同乱数で作成する.水面からの蒸散,地下浸透は単純化のため,ここでは考慮しない.河川に流入する物質量をこのように単純化すると,本流流量と本流水に含まれる物質量が計算できる.
【3】河川は浄化機能を有する.その浄化機能を観測点間距離に依存する.この機能は一次関数で表現できるとする.実河川は上流と下流で河川構造が異なり,水中の物質量も違うと思われる.そこで我々は多摩川の1/4000縮尺の航空写真を羽村堰から河口まで調べた.その結果,日野橋から下流では河川構造の相異は堰の存在以外の顕著な相異は見られなかった.そのような状況ならば河川浄化機能を観測点間距離に比例するとしても良いように思われた.
以上【1】〜【3】の仮定を用いて人工河川水質の変化を観測点で離散的に表現した.それは最も簡単なモデルである.以下,単純離散モデルという.さらに複雑なモデルは第3,4節で研究する.以下にアルゴリズムの詳細を示す.
<1> モデルの各量は無次元である.
<2> 河川はN-1の観測点{2,3,4,...,N}={i}で水質値が定義される.観測点1を本流の始点とする.始点を加え{1,2,3,4,...,N}={Xi}で河川の各点を番号づけする.1側を上流側とする.
<3> 観測点間距離を[0.5, 2]区間の一様乱数(以下,乱数)とする.{di}と表記する.diXi-1とXiの距離である.本区間は多摩川の支流の流量[7]から決定した.
<4> 支流は{2,3,4,...,N}={Bi}, 1<iと番号づけされる.支流は川だけでなく下水道水などを含む本流に流入する水の総称である.BiXiに流入する.
<5> 支流の流量(=本流流入量)を[0, 1]区間の乱数とする.これを{fi}と表現する.1<iである.本流の流量は

である.
<6> 支流Biの水質を[1, 1+ai]区間の乱数とする.{ai}={0, 1/N, 2/N,...,1}である.これを{pi}と表示する.1<iである.このようにすると支流の水質は下流側で悪化することになる.それは多くの河川で見出される状況である.シミュレーションで支流効果を明らかにするため上流側の支流水質の最低値を1とした.
<7> 本流Xi地点の水に含まれる物質量は,

である.1<iである.s1=0とする.ここでkiは観測点Xi-1とXiの間の河川水浄化機能である.ここで仮説【3】「浄化機能は水に含まれる汚染物質量によらない.距離に依存する」が採用されている.kiは普通負数である.またさらに近似してk1=k2=...=K(定数)とすることもできる.観測データ数が少ない(約20以下)では定数が妥当である.

制約条件は0< siである.
<8> 本流の水質基準値{qi}

である.
このアルゴリズムの特徴について議論する.式(3)はベクトル{si}, {fipi}, {di}, 2<=iが与えられればKについて解くことができる.すなわち本流の水質指標値,支流水量,支流水質指標値,観測点間距離の観測があれば,河川浄化機能がまだ働いているか否かを検討できる.最低観測点数は式(3)が与えられていれば2である.しかし式(3)は与えられていないので観測数は多いほど良い.神経回路網は不明の(現象)支配方程式,

をデータから生成し,その偏微分,

から河川浄化係数Kを求めるのが本アルゴリズムの目的である.
「データから生成」できる可能性は以下のように考えられる.神経回路網に入力する一つの説明変数をxiWxij, Wyjkを第1,2層間のニューロン間の結合荷重,第2,3層間の結合荷重,添字i,j,kを第1,2,3層のニューロンを示すとし,第1,2層の最後のニューロンそれぞれ1個をbias neuron,値を1とする.添字ijkは一対一対応しているとする.第2層のニューロンに入力する値をhjとすると,

第2層の出力qjは,第2層のニューロン・エミュレーション関数,シグモイド関数ffとして,

第3層からの出力ykは,

となるから,式(7)で式(3)の項間の+演算子が生成され,ff関数部ではバイアスによりaverage(hj)→0となればffの直線部を使い式(3)の線形結合が再現される可能性がある.式(9)はその直線部の拡大/縮小作用なので再現を助ける.以上一変数で説明したが,多変数でも全く式(7)の作用は同じで+,+,...という結合が生成される.
F(si-1, fipi, Kdi)/∂diが定数でなく関数G(di)={Gi}となる場合はアルゴリズムの一貫性(self consistency)が満足されていないことになる.定数でない度合を距離D

とする.Nはベクトル{Gi}の要素数である.

2. 2 神経回路網のための単純離散モデルデータ

第2.1節のアルゴリズムで作成した神経回路網用単純離散モデルデータの一例(K=-0.2の場合)をTable 1に示す.

Table 1. Neural network learning data for simple discrete river model. The d1,2,3 are descriptors. The d1 is fipi. The term fipi is substance amounts from branches into the main stream. Where the branches are a generic name of branch rivers' water and sewage and rainwater. The d2 is distances between observation points, i.e., di. The d3 is substance amounts at Xi point in main stream, i.e., si-1. The T is teaching data, si.
K=-0.2d1d2d3T
10.73991.48030.00000.4438
20.47951.26530.44380.6703
30.94681.04360.67031.4083
40.54220.62441.40831.8256
50.27200.86821.82561.9240
60.95070.51271.92402.7721
70.05661.61932.77212.5048
81.03080.72412.50483.3908
91.10931.53783.39084.1926
100.14381.14894.19264.1066
111.18381.97094.10664.8962
120.54341.79544.89625.0806
131.08540.98825.08065.9684
141.38160.64955.96847.2201

本流水中のK=0,-0.1,-0.2,-0.3,-0.4の場合の各観測点の物質量(教師データ)をFigure 1に示す.河川浄化係数-K(0-0.4)による変化をFigure 1中に示した.各河川で観測された水質指標(たとえば多摩川のBODデータなど)と類似している.
単位はデータを乱数で作成しているため具体的に[kg]のように書けないので,質量の意味があることを示すために[mass]とした.以下体積を[volume]などとする.


Figure 1. Substance amounts in the water of a virtual river. The vertical axis is the substance amounts that have a dimension of [mass]. The [mass] is a generic-dimension of [kg] etc. The horizontal axis is sampling points, which are numbered. The points are located by the accumulated distances. The left side is the upper stream. Whole data are generated by uniform random numbers. In the figure, digits are the river purification coefficient, K, which are multiplyed by (-).

支流から流入する物質量をFigure 2に示す.


Figure 2. Inflow substance amounts from the branch in the virtual river system. The vertical axis is the amount whose dimension is [mass]. The branch is one at an 'observation' point in the main stream. The branch is a generic name of all inflow. Therefore, we avoid no branch stream. The horizontal axis is the number of the branch, 1<=Bi<=14. We set water quality of the branch is to be clean at upper stream; and we set many kinds of branches flow into the main stream.

我々はFigure 2において,様々な水が流入するように,かつ上流ほど清浄な水が流入するようにした.Table 1のデータからデータを作成したモデル式(2)を神経回路網に与えないで河川浄化係数Kを推測するのが目的である.一種の逆問題である.

2. 3 河川浄化係数Kの計算

第2.1節のアルゴリズムで河川水質のデータを生成し,観測点Xiへの支流からの流入物質量fipiXi-1観測点との距離di,本流Xi-1点の物質量si-1を説明変数(descriptor)とした.本流Xi点の物質量siを教師データとした.神経回路網にとって説明変数は入力データであるのでそのように記す.Table 2に3説明変数,1教師データ,各データ数14の場合の神経回路網学習パラメータを示す.学習方法はback propagation + reconstruction learningを用いた.本計算の目的は,学習が停留状態になった神経回路網から河川浄化係数Kを求める具体的方法を示すことである.

Table 2. Learning parameters for vector data in Table 1. The number of neurons on the second layer is only 3 that are decreased from initial value 8. One of the three neurons is a bias. The bias neuron is on the first layer, too.
# of data(descript., teach.)(14×3, 14)
Neurons on 1,2,3-layers4, 3, 1
Emulation on 2nd-layersigmoid
Emulation on 3rd-layerlinear
# of learning3K
Learning Constant0.2, 0.2
Reconstruction60 times
Erasing factor0.04
Square of BP error0.025


Figure 3. Maximum and average BP error in learning iterations. The vertical axis is the logarithm of the error. The horizontal axis is the learning iterations' number. The learning data are generated by the coefficient, K=-0.2.

Figure 3に学習時の教師データとの差の変化を示す.Figure 3において左端が1回,右端が3K回,150回ごとに点描した.河川浄化係数Kは0〜-0.4まで5段階に変化させたがKの値によらずFigure 3と同様な変化が得られた.


Figure 4. Neural network outputs and the teaching data after learning iterations of 3K. The vertical axis is substance amounts that are calculated by K=-0.2. The dimension is [mass]. The horizontal axis is the observation points, Xi. The left side is the upper stream. The black points are the teaching data, and the bend lines are the outputs. The error line is the difference.

Figures 3, 4から判断すると,観測誤差の無い場合,神経回路網の計算精度はO(-2),2桁である.河川浄化係数Kは0〜-0.4まで5段階に変化させたがKの値によらずFigure 4と同様な変化が得られた.単純離散モデルと神経回路網法は環境解析に十分な精度を備えている.


Figure 5. Partial derivatives of three descriptors and a bias. The vertical axis is index for changes of substance amounts. The changes are calculated by partial derivatives of a neural network on K=-0.2. The dimension is [mass]. The horizontal axis is the observation points, Xi. The left side is the upper stream. The curves d1-3 are partial derivatives of {si-1}, {di}, {fipi}. The B shows a partial derivative of a bias.

河川浄化係数Kは0〜0.4まで5段階に変化させたがKの値によらずFigure 5と同様な変化が得られた.折線d1,d3はX1-14で期待値1である.d2はそのまま河川浄化係数値である(河川浄化項が線形でその偏微分を取ったのであるから,偏微分値=Kである).負の値は観測点間距離が水中の物質量を減じる方向(浄化方向)に作用したことを示す.折線Bは神経回路網の非線形関数機能の変化の一指標である.各説明変数のself consistencyをTable 3に示す.期待値は共に0である.

Table 3. Self consistency of the simple discreate model as for descriptors and a bias. The d1-3 are same as Figure 5.
d1d2 (K)d3B
0.00630.00020.00570.0046

各説明変数の偏微分係数値の変化は少なく神経回路網が正しく式(3)を再現したことを示す.特にO(-4)の河川浄化係数Kの値は信頼できる.Table 4K=0〜-0.4の各djの偏微分係数値を示す.d2が本計算で求める河川浄化係数である.expectationはその期待値である.

Table 4. Partial derivative coefficients for descriptors d1-3 and a bias.
expectation0.00-0.10-0.20-0.30-0.40
d11.071.071.061.061.05
d2 (K)0.02-0.08-0.19-0.29-0.40
d31.021.021.031.031.04
B-1.63-1.28-0.91-0.43-0.13

Table 4は河川浄化係数Kが高精度で計算されたことを示す.神経回路網に学習させたデータに誤差が無く想定モデルに合致したデータで,かつ不足が無いならば±2%の精度で計算可能である.

3 堰効果を導入した離散モデル

3. 1 堰モデル定義

多摩川の航空写真調査から無視できないと考えられる堰の効果を単純離散モデルに導入する.元来,堰は農業用水等の取水目的に作られたが堰は小さなダムであるから,そこに水が滞留し化学反応により水質指標値が変化する.その変化は考慮する必要がある.
堰の性質として1. 水の滞留時間,2. 水温,3. 水深,4. pHなどのデータが公表されていることが望ましいが,そういうデータは入手困難である.その場合,航空写真[8]で堰の規模と水面の状況を調べる.多摩川,上中流域の8堰で調査したところ規模はほぼ同じと判断した.水面の状況から滞留時間,水温についてもおおむね同じであろう.外国の大河川のように著しい相異がある場合は以下の方法は適用が困難であるが「同じ」場合では「各堰の特徴を無視し全堰で同じ」と近似する.そうすると神経回路網のために堰を表現する形式としてscaling[9]により堰の存在非存在のみが有意となる.したがって0/1要素のベクトルとなる.このベクトルの存在は構造活性相関研究分野[10]では,ダミー変数として用いられてきた.構造活性相関では構造を表現する適切な物理化学的測定手段が不明の場合,あるいは従来指標では表せない効果を想定して使用されている.それらは主要な原因を説明する変数として導入されず,また少数の大きな外れ値(outlier)を説明する変数としては用いられていない.
環境問題ではダミー変数は構造活性相関的にも外れ値,突発的事故などによる環境負荷の一時的増大の解析にも使用できる可能性がある.本論文は前者の場合を扱う.後者は今後の課題である.我々の考えている外れ値を付録に示した.
ダミー変数は多用すると危険性である[11]が,注意して使用すれば多摩川,ドナウ川の水質問題については有効であった[12, 13].注意深く使用することが前提であるが,我々はダミー変数は常時有効データ不足である環境科学の研究にもっと採用されて良い方法と信じる.
我々はダミー変数形式の堰データを導入した堰モデルを提案する.アルゴリズムは以下である.
<1>-<6> まで第2.1節と同じである.
<7> 堰を観測点Xi, i=2,6,10,14,...に設定する.
<8> 本流Xi地点の水に含まれる物質は

である.1<iである.s1=0とする.ここでKは観測点Xi-1とXiの間の河川水浄化係数(負が浄化方向,正が汚染方向)である.Lは堰の浄化係数,hiは堰の存在/非存在を示すダミー変数である.制約条件は0<siである.
<9> 本流の水質基準値{qi}は,

である.

3. 2 神経回路網用堰モデルデータ

第3.1節のアルゴリズムで作成した神経回路網用堰モデルデータの一例(K=-0.2の場合)をTable 5に示す.

Table 5. Dam model data for a neural network. These data are generated by uniform random data in algorithm <1-9>. Descriptor d1: inflow substance amounts at observation points, {fipi}. Descriptor d2: distances among observations, {di}. Descriptor d3: substance amounts at upper observations, {si-1}. Descriptor d4: dummy variable for existence of a dam. Teaching data: {si}. K=-0.2, L=-0.15;
d1d2d3d4T
10.73991.48030.000000.4438
20.47951.26530.443810.5203
30.94681.04360.520301.2583
40.54220.62441.258301.6756
50.27200.86821.675601.7740
60.95070.51271.774012.4721
70.05661.61932.472102.2048
81.03080.72412.204803.0908
91.10931.53783.090803.8926
100.14381.14893.892613.6566
111.18381.97093.656604.4462
120.54341.79544.446204.6306
131.08540.98824.630605.5184
141.38160.64955.518416.6201

d1-4の四種類の説明変数(ベクトル)から河川と堰の浄化係数K,Lを神経回路網を使用して数値的に計算することが目的である.これは第1.3節の逆問題よりも説明変数が多く困難な問題である.堰モデルで各Xi点で計算した河川水中の物質量をFigure 6に示す.


Figure 6. Substance amounts calculated by the dam model. The vertical axis is the substance amounts that have a dimension of [mass]. The horizontal axis is sampling points, which are numbered. The points are located by the accumulated distances. The left side is the upper stream. Whole data are generated by uniform random numbers. In the figure, digits are the river purification coefficient, K, which are multiplyed by (-). Ther dam coefficient, L, is fixed as -0.15.

Figure 6において,実際の河川(多摩川)で観測されたデータと類似しているK=-0.1/-0.2がシミュレーション用として適当に思われる.K=-0.4では物質量が2番目の観測点で-0.013となり,わずかに式(5)の制約条件に触れる.神経回路網シミュレーションは可能であるが,実際の河川で採用するのは適当でない.
Table 6K=0〜-0.4の各djの偏微分係数値を示す.d2,4が本計算で求める河川と堰の浄化係数である.expectationはその期待値である.

Table 6. Partial derivative coefficients for descriptors d1-4 and a bias in the dam model. Descriptor d1: inflow substance amounts at observation points, {fipi}. Descriptor d2(K): distances among observations, {di}. Descriptor d3: substance amounts at upper observations, {si-1}. Descriptor d4(L): dummy variable for existence of a dam. "B" is that of a bias neuron.
expectation0/-0.15-0.1/-0.15-0.2/-0.15-0.3/-0.15-0.4/-0.15
d11.001.011.011.011.01
d2 (K)-0.04-0.12-0.22-0.31-0.40
d31.001.001.001.011.01
d4 (L)-0.13-0.13-0.13-0.13-0.14
B-4.25-3.52-2.73-1.96-1.30

Table 6において,d1,3は支流からの流入物質量,本流の上流側から流入する物質量である.河川と堰以外に物質量を増減する要因は無いので,それらの期待値は1である.B欄は神経回路網の非線形関数機能の変化を示す.本表は河川と堰の浄化係数K,Lが高精度で算出されたことを示す.神経回路網に学習させたデータに誤差が無く,想定モデルに合致しデータで,かつ不足が無いならば±4%の精度で河川浄化係数が,±2%の精度で堰の浄化係数が同時に予測可能である.Table 7に各説明変数のself consistencyを示す.d2(K), d4(L)ともO(-4)でモデルの一貫性に問題はない.以上の調査から河川水質データの離散形は環境解析に充分適用できると考えられる.

Table 7. Self consistency of d1-d4 and bias, in case of K=-0.20, L=-0.15.
d1d2 (K)d3d4 (L)B
0.00300.00010.00270.00010.0247

4 地下浸透を導入した離散モデル

前節までのモデルの検証では説明変数に不足は無かった.モデルに多くの要因を説明変数として導入すると現象記述機能は高くなる.一方データを用意することが困難になる.説明変数データが不足した状態のモデル精密化の効果が神経回路網出力にどう影響するかを調査する.

4. 1 地下浸透モデル定義

堰モデルに地下浸透効果を導入する.地下浸透データの実測はほとんど無い.本節では地下浸透データを説明変数に加えずに神経回路網を学習させる.そのようなデータ不足の状態で学習しても神経回路網は自動関数適合機能により学習を完了する.それは「不十分な説明変数で現象を再現する」という状態である.その状態の現象再現性の程度を調べることは,我々が不十分な情報から結論を導く時の「結果の信頼性」を調査することに通じる.環境問題解析はこのような状況にあるのではないかと思われる.地下浸透モデルのアルゴリズムは以下である.
<1>-<7> まで第3.1節と同じである.
<8> 地下浸透量を[0, 0.5]区間の乱数とする.理科年表環境編[14]の日本の河川データを見る限り妥当なオーダである.地下浸透は流出する支流として扱う.{yi}と表記する.
<9> 本流Xi地点の水に含まれる物質は,

である.1<iである.s1=0とする.ここでKは観測点Xi-1とXiの間の河川水浄化係数である.Lは堰の浄化係数,hiは堰の存在/非存在を示すダミー変数である.qiは<10>で定義される本流の水質基準値である(第2.1節,第3.1節の式(4)と同じ).制約条件は0<siである.
<10> 本流の水質基準値{qi}はqi=si/gi, 1<iである.

4. 2 神経回路網用地下浸透モデルデータ

第4.1節のアルゴリズムで作成した神経回路網用地下浸透モデルデータの一例(K=+0.1, L= -0.15の場合)をTable 8に示す.

Table 8. Learning data for the underground-penetration model. These data are generated by uniform random data in algorithm <1-10>. Descriptor d1: inflow substance amounts at observation points, {fipi}. Descriptor d2(K): distances among observations, {di}. Descriptor d3: substance amounts at upper observations, {si-1}.Descriptor d4(L): dummy variable for existence of a dam. Teaching data is T, {si}. The descriptor for {qi-1yi} is not given.
d1d2 (K)d3d4 (L)T
10.73991.4803000.8879
20.79491.14920.887911.4106
30.09190.68481.410601.5015
40.25820.86821.501501.8440
50.48661.72441.844002.4889
60.20621.48512.488912.5818
71.03361.53782.581803.6001
80.02020.71373.600103.3069
91.04371.18853.306904.1506
101.03530.98824.150615.0930
110.44661.96175.093005.3559
120.15310.62785.355905.5334
130.57651.85775.533406.1387
140.07241.47726.138716.0167

4. 3 説明変数不足の場合の地下浸透モデルの解析

Table 8のデータは地下浸透に関する説明変数(式(13)のqi-1yi項)が不足している.それとTable 9の学習パラメータを用いて河川の神経回路網解析を行いFigure 7の結果を得た.

Table 9. Learning parameters for the underground-penetration model. The number of neurons on the second layer is 3 that are decreased from initial value 8. One of the three neurons is a bias. The bias neuron is on the first layer, too.
# of data(descript., teach.)(14×4, 14)
Neurons on 1,2,3-layers5, 3, 1
Emulation on 2nd-layersigmoid
Emulation on 3rd-layerlinear
# of learning3K
Learning Const.0.2, 0.2
Reconstruction60 times
Erasing factor0.04
Square of BP error0.034*
* in case of K=+0.1, L=-0.15.


Figure 7. Neural network outputs and the teaching data after learning iterations of 3K. The vertical axis is substance amounts that are calculated by K=+0.1, L=-0.15. The dimension is [mass]. The horizontal axis is the observation points, Xi. The left side is the upper stream. The black points are the teaching data, and the bend lines are the outputs. The error line is the difference.

河川浄化係数Kは-0.1〜-0.3まで5段階に変化させたがKの値によらずFigure 7と同様な変化が得られた.説明変数不足にもかかわらず神経回路網の出力値は期待値を精度良く再現する.神経回路網に自動関数適合機能が存在するためである.これは環境問題解析に重要な問題で,注意して解析しないと誤った結論を導き易い.Figure 8に同神経回路網から得られた各Xi点の偏微分係数を示す.


Figure 8. Changes of partial derivatives of four descriptors. The vertical axis is an index for changes of substance amounts. The changes are calculated by partial derivatives of a neural network on K=+0.1, L=-0.15. The dimension is [mass]. The horizontal axis is the observation points, Xi. The left side is the upper stream. The curves d1-4 are partial derivatives of {si-1}, {di}, {fipi}, and dummy variable for existence of a dam.

河川浄化係数Kは0.1〜0.3まで5段階に変化させたがKの値によらずFigure 8と同様な変化が得られた.説明変数d1,d3はX1-14で期待値1である.Table 10に各K, L値の期待値と神経回路網の計算した値を示す.

Table 10. River and dam purification coefficients calculated by the underground-penetration model.
expectation0.1/-0.150.0/-0.15-0.1/-0.15-0.2/-0.15-0.3/-0.15
d10.990.990.990.991.00
d2 (K)0.01-0.08-0.16-0.25-0.33
d2(K+ qi-1yi)0.07-0.03-0.13-0.23-0.33
d30.960.960.960.960.97
d4 (L)-0.13-0.14-0.15-0.15-0.16
B -4.08-3.31-2.62-1.90-0.97

Table 10において,河川浄化係数値(d2, K)がK>-0.30で正しくない.しかし堰の係数(d4, L)は正しい.河川浄化係数値は地下浸透項の平均値(-0.0287)を加算した値(d2, K+qi-1yi)にほぼ等しい.シミュレーションではなく実測定値で解析する場合,データ不足の場合は正しくない結果となる恐れがある.それをself consistency値で判定できるか否かを調べる.同値をTable 11に示す.

Table 11. Self consistency of the underground-penetration model.
KLd1d2(K)d3d4(L)B
+0.1-0.150.00290.00000.00260.00010.0482
0.0-0.150.00240.00000.00240.00010.0268
-0.1-0.150.00220.00010.00220.00010.0159
-0.2-0.150.00250.00010.00260.00010.0102

Self consistencyの値ではK=+0.1の場合のバイアス・ニューロンに関する偏微分係数値が異常をしめしている.これは神経回路網で不足する説明変数を補う何らかの関数機能変化が起こっていることを示す.しかし,他の偏微分係数値はTable 7と比較して特に異常ではない.K=+0.1の場合は不足が著しいため検出できたと思われる.ゆえにself consistency値だけではデータ不足状態を検出するには不十分である.この結果は神経回路網解析の限界と思う.バイアス・ニューロンに関する偏微分係数をFigure 9に示す.


Figure 9. Changes of partial derivatives for a bias neuron. The vertical axis is amplitude of partial derivatives calculated by K=+0.1, L=-0.15. The horizontal axis is the observation points, Xi. The left side is the upper stream.

地下浸透効果は各Xi点ごとに違い,河川浄化係数にパターン的に類似している.そのためd2が本来神経回路網に入力されるべき地下浸透量を包含するように神経回路網の中で作用したと考える.その証拠がTable 10のd2(K+ qi-1yi)欄とFigure 9K=0.1の曲線である.確かにK=0.1のときd2の値は最大誤差を示している.一方,堰の効果はTable 12に示したように相関係数が小さく影響が少ないと判断する.

Table 12. Correlations between river/dam purifications and underground-penetration terms.
K+0.100.00-0.10-0.20-0.30
d2 vs qi-1yi-0.16-0.17-0.18-0.20-0.23
d4 vs qi-1yi0.030.03-0.00-0.04-0.13

4. 4 観測点で異なる河川浄化係数kiの計算

説明変数不足の場合の神経回路網の出力の特徴をさらに調査するため,単純離散モデルを再検討する.第2.1節では式(3)により河川水質のデータを生成し,観測点Xiへの支流からの流入物質量fipiXi-1観測点との距離di,本流Xi点の物質量si-1を説明変数とし,本流Xi点の物質量siを教師データとし,それらから河川浄化係数Kを求めた.未定係数K以外は全てベクトルの要素で与えられている.従って,それらの各要素を満足する平均的な値としてKが求まる.
本来の式(2)では観測点の河川浄化係数kikidiという両項未定形である.故に3説明変数{si-1, fipi, kidi}の学習データでは原則としてkiを求めることはできない.それでもその3変数ベクトルにより神経回路網を学習し,その偏微分係数からkiを求めようとすると,

となる.神経回路網は非線形ではあるが複数のデータから尤もな未定定数を求める方法の一つであるので,平均化されたkiが求まる.しかし式(7)は個別のkiを求められる可能性を示唆する.ただし可能性であるので数値計算で検証する.神経回路網のki算出機能について追跡容易にするため,{ki}を特徴ある数列,等差数列とする.河川浄化係数Kを使用して,

とする.このとき式(7)のaverage(ki)項は項数がN+1なので,

ここで改めて

と書く.
式(16)より神経回路網の学習をデータ数N+1個ではなく要素mを除去したN個で行うと,average(ki)mN+1の場合より除去したkm分だけ変化する.これをdmとすると,

dmが高精度に計算できるのなら,

である.負号は基準の取り方による.
式(20)は小さい値dmN倍するので,実測の誤差を含むデータでは実行困難である.無誤差のシミュレーションならば可能性がある.このような観測集合の部分集合でモデルの性質を調べる方法は,モデルの正当性をその出力から判断する一方式としてleave-one-out法[15]と呼ばれて来た.式(20)の方法は神経回路網の出力値ではなく偏微分係数から算出される項の性質に関する点が本来のleave-one-out法とは異なる.
式(20)は小値dmN倍しているので除去要素mの単独計算だけからは精度良くkmが求まらない.ここで{km}の等差数列の性質を利用し,除去要素mの集合{m; m =0,1,...,N}を考え各々のkmを求め,平均変化(勾配という)を求めた.式(20)のaverage項を使い,

となる.ゆえに{km}の勾配2K/Nが求められる.3変数の尤度方程式が複数あることを利用して,4説明変数の情報を得ることになるので,一種の平均量が求められることを限度と考える.
ここでdmをどの説明変数について調べれば良いかを議論する.常識的には式(2)のkidi項であるが,神経回路網に入力するのは値だけであり,その値がkidiから計算されたのものであることを提示していない.従って,神経回路網は他の説明変数と区別することが出来ず,どの説明変数の偏微分係数に表れるかは不明である.ここで「神経回路網は複数データの最尤定数求解法」であることを利用する.従って偏微分係数は本来ならば定数になるので,そうなっていない係数がその候補である.
環境問題解析では有効観測数は多くないので観測データの部分集合化は計算精度低下をもたらす.leave-one-out法の拡張leave-n-out法は不可能な場合が多いが,それが可能ならばTable 13を利用すれば勾配に関する情報は多くなる.この単純離散形式の拡張型を双未定係数モデルと呼ぶ.
神経回路網は隠れ層のニューロン数を多くすれば学習点に関し教師データと出力値の差は0に漸近する.そのため神経回路網を使用したデータ解析において外れ値(outlier)を見出すのは難しくその議論は少ない.しかも本論文の神経回路網は隠れ層を第2層のみとし,そのニューロン・エミュレーション関数をsigmoid関数とし,かつ第3層の同関数を線形関数とした(Tables 2, 9).かつreconstruction learningにより第2層のニューロン数を2まで減少させ,神経回路網の非線形fitting能力を限定した(Figure 3参照).さらにleave-one-out法を使用している.従って,外れ値の存在が疑われる「観測データ」の場合,それが与える影響を除去しなければならない場合がある.ただし外れ値の議論は本論文の趣旨とは少し離れるので付録Aに記した.

4. 5 河川浄化係数kiの計算例

第4.4節の方法を検証するため,第2.1節のアルゴリズムで河川水質のデータを生成し,観測点Xiへの支流からの流入物質量fipiXi-1観測点との距離di,本流Xi-1点の物質量si-1を説明変数とした.

Table 13. Leave-n-out method expression for the slope of arithmetic progression.
leave-n-outexpression
12K/N*
26K/(N-1)
312K/(N-2)
420K/(N-3)
530K/(N-4)
642K/(N-5)
756K/(N-6)
* N is the number of data.

河川浄化係数kiについては第4.4節の式(8)を用いた.このときK=-0.1, N=14であった.本流Xi点の物質量siを教師データとした.Table 14に3説明変数(kidiは1変数である),1教師データ(各データ数14)の場合の神経回路網学習パラメータを示す.Figure 10K=定数の場合との相異を図示した.定数との相異は図示すると僅かである.神経回路網の学習方法はback propagation+reconstruction learningを用いた.本計算の目的は,学習が停留状態になった神経回路網から第4.4節の方法で河川浄化係数{ki}の勾配の平均値average(ki)= -0.1*2/(14-1)=-0.015を求めることである.

Table 14. Learning data to calculate variable river-purify coefficients {ki}. The d1,2,3 are descriptors. The d1 is fipi. The term fipi is substance amounts from branches into the main stream. Where the branches are a generic name of branch rivers' water and sewage and rainwater. The d2 is the multiplyed values, kidi, which are calculated by distances and river-purify coefficients. They are unknown variables; however, they are processed as a descriptor here. The d3 is substance amounts at Xi point in main stream, i.e., si-1. The T is teaching data, si.
d1d2d3T
10.73991.48030.00000.7399
20.47951.26530.73991.1999
30.94681.04361.19992.1146
40.54220.62442.11462.6279
50.27200.86822.62792.8465
60.95070.51272.84653.7578
70.05661.61933.75783.6649
81.03080.72413.66494.6177
91.10931.53784.61775.5377
100.14381.14895.53775.5224
111.18381.97095.52246.4030
120.54341.79546.40306.6426
131.08540.98826.64267.5456
141.38160.64957.54568.7973


Figure 10. The difference between variables and a constant for the river purification. The vertical axis is the substance amount, whose dimension is [mass/volume]. The horizontal axis is ordered observation points, where the left side is upper stream. Total river purification power is same for the both conditions; i.e., K=-0.1=Ski/N.

学習が停留状態になった神経回路網の出力値とその期待値,その差をFigure 11に示す.図は不完全な学習データであっても神経回路網が高精度で現象を再現することを表している.それは神経回路網の非線形fitting機能であるが,これで現象が説明できたと考えると危険である.真に説明できたか否かは偏微分係数を調べる必要がある.Figure 12に各説明変数とバイアスの偏微分係数を示す.


Figure 11. The difference between neural network's outputs and the expectations. The diamond-shaped black points are expectations, and the NN curve is outputs, which show the fitting-ability of the neural network. The figure shows fitting-ability even if the learning data are incomplete.


Figure 12. Partial derivatives for 3-descriptors and bias. The d1-3s are in Table 11. The B is partial derivative for a bias on the first layer. The vertical axis is amplitude of partial derivatives, whose unit is [mass/volume]; where the denominator is uniform scaled derivatives in [0.05, 0.95]. The bias is fixed values; however, it is same characters to descriptors for a neural network. Therefore, we calculate the bias distribution to estimate neural network's function.

Figure 12において,各説明変数の値は「神経回路網は複数データの最尤定数求解法」であるためには一定である必要がある.Self consistencyをTable 15に示す.

Table 15. Self consistency of extended discrete river model.
d1d2d3B
0.00450.00000.00410.0141

Self consistencyが不十分な係数はバイアス・ニューロンの係数である.この変化は神経回路網の関数機能がデータにより変化している度合を示す.そこでは不足する説明変数を神経回路網の内部で補っている可能性を示している.その内部補正量の大きさは,神経回路網の出力値をyi(添字iはデータの順序番号である)とするとき,

である条件から求められる.定数値は不明であるので,

と仮定し,

とする.ベクトル{Pi}の要素値の傾きを最小二乗法で直線にfittingして調べる.それをFigure 13に示す.計算された値は-0.0014で期待値-0.015に符号は一致する.式(23)の仮定は傾きには関係しないので,この結果は受け入れることができる.ただし最小二乗法の直線fit時のR2=0.79は理論値1に近くない.それはこの方法の限界を示すように思われる.
以上leave-one-out法とバイアスの偏微分係数から説明変数不足の場合でも「データのみから不足する説明変数の性質を逆算できる」可能性を示した.それが可能であるのは「モデルが推測でき,かつそれが正しい場合」である.


Figure 13. Change of {Pi} that is derived from eq. (24). The change relates the slope of {ki}; this is a hypothesis. Calculated slope is -0.0014, the expectation is -0.015.

説明変数d2(kidi項に相当する)のleave-one-out法の偏微分係数の変化と式(20)からはkmが直接求められる可能性がある.説明変数d2の偏微分係数は各観測点でほとんど変わっていない(Table 15のd2)ため,その可能性がある.計算結果をFigure 14に示す.河川浄化係数が各観測点ごとに違い,かつその変化がおおむね直線的であることは正しく算出されている.しかし傾きが0.0044で期待値-0.015とは違う.この結果は第4.4節の式(20)方法の限界を示している.


Figure 14. Coefficient {ki} estimated by descriptor-d2 from leave-one-out method and eq. (20). The vertical axis is the river purify coefficient, [mass/distance]. The horizontal axis is the observation points. The lefit side is upper stream.

5 結論

現代の都市部を流域とする河川の水質について離散形の四モデルについて考察した.その離散モデルで記述された人工河川の水質を神経回路網により解析可能であるか否かを研究した.神経回路網は非線形現象の解析に有用であるが,その機能が学習により獲得できるためには説明変数と観測値に不足があってはならない.不足とは観測データが本来あるべきことが自明であるのに何らかの理由で不明な場合と,現象を説明する観測自体がされていないか,そういう観測自体の存在不明の三種類がある.前者を欠測という.欠測を含んだデータの情報処理は近年かなり研究されているが後二者はほとんど検討されていない.我々は環境問題の真の理解のために後者を回避できない問題と考える.同時にそれは神経回路網出力の検定に関わる問題である.本論文の目的はこの観測データ不足問題を考察することであった.
本論文では神経回路網を使用した非線形多変量解析のデータ不足状況における挙動を調査し同回路網解析の適用限界を明らかにするため,観測データに含まれる測定誤差を除外することにした.測定誤差の統計的性質は確かでない.それを基に神経回路網の適用限界の議論は困難である.この測定誤差問題を回避するため河川モデルを考案し,乱数を基に人工河川を定義し,そのデータを基に神経回路網の適用限界を明らかにする方針を採用した.
我々の提示した河川モデルは単純離散モデル,堰モデル,地下浸透モデル,双未定係数モデルである.単純離散モデル,堰モデルについてはデータが完全な場合を計算し河川と堰の水質浄化力を誤差±2〜4%の精度でシミュレーションできた.このレベルが神経回路網の本来の能力と我々は考える.環境問題を議論するのに必要十分である.
地下浸透モデル,双未定係数モデルでは,あえて必要な説明変数データを省略してデータ不足な場合を作った.このデータ不足状態であっても神経回路網は精度良く河川水中の物質量をシミュレーションする.しかし偏微分値を計算すると,地下浸透モデルでは河川の浄化機能が正しく計算できていないことが明らかになった.その理由は説明変数にない地下浸透効果を河川浄化機能が補完したためである.両者の物理化学的性質はモデルでは類似している.双未定係数モデルでは河川浄化機能の平均的変化を神経回路網のバイアス部が補完した.その補完機能を利用し偏微分係数から河川浄化機能の平均的変化を逆算できることを示した.ただし理論上の可能性だけである.数値計算では符号のような大まかな特徴だけが逆算できた.以上の結果は与えられていない説明変数の効果がデータのみからある程度は推定できることを示す.それが可能であるのは「モデルが推測でき,かつそれが正しい場合」であるが,このような神経回路網の未知の使用法が明らかになったことは環境問題の議論を定量化する観点からは有用である.

参考文献

[ 1] 大垣眞一郎監修, 河川環境管理財団編, 「河川と栄養塩類 管理に向けての提言」foundation of river and watershed environment management, 技報堂出版 (2005), ISBN4-7655-3403-0.
[ 2] 日本水環境学会編, 日本の水環境3 関東・甲信越編, 技報堂出版 (2000), ISBN4-7655-3167-8.
[ 3] 末次忠司,河原能久,賈仰文,倪广恒, 都市河川流域における水・熱循環の統合解析モデルの開発, 土木研究所資料, 第3713号 (2000), 河川部都市河川研究室 ISSN 0386-5878.
[ 4] 賈仰文,倪广恒,河原能久,末次忠司, 都市河川流域の水循環解析と雨水浸透施設の効果の評価, 水工学論文集, 44, 151-156 (2000).
[ 5] Goloka B. Sahoo, Chittaranjan Ray, Jack Z. Wang, Stephen A. Hubbs, Rengao Song, Jay Jasperse, Donald Seymour, Use of articial neural networks to evaluate the effectiveness of riverbank filtration, Water Research, 39, 2505-2516 (2005).
www.elsevier.com/locate/watres
[ 6] 渡辺美智子,山口和範, EMアルゴリズムと不完全データの諸問題, 多賀出版, ISBN4-8115-5701-8.
東洋大,渡辺教授のWeb資料:
http://stat.eco.toyo.ac.jp/~michiko/em/emohp(0)/emohp.ppt
[ 7] 公共用水域水質測定結果データ集:
http://www2.kankyo.metro.tokyo.jp/kansi/mizu/sokutei/sokuteikekka/suishitu.htm. (CSV形式)
東京都下水道局HP: http://www2.kankyo.metro.tokyo.jp/kansi/mizu/sokutei/sokuteikekka/suishitu.htm.
東京都下水道局事業概要 平成16年版第4章「流域下水道主要事業の展開」http://www.gesui.metro.tokyo.jp/gijyutou/jg16/jigyougaiyou16/no4.pdf. (PDF形式)
[ 8] 国土地理院(Geographical Survey Institute)
http://mapbrowse.gsi.go.jp/airphoto/indexmap200k/5339/533934.html.
1/4000〜1/20000大縮尺航空写真(無料)の一例
http://www.ikutoko.com/
[ 9] 青山智夫,神部順子,長嶋雲兵, 欠測データ集合を扱う神経回路網法CQSAR:Compensation Quantitative tructure-Activity Relationshipsの開発, J. Comput. Chem., Jpn., 投稿中(2005.10.27).
[10] 構造活性相関懇話会, 薬物の構造活性相関 ドラッグデザインと作用機作研究への指針 化学の領域増刊122号, 南江堂, 東京 (1979).
[11] 青山智夫,神部順子,長嶋雲兵, 階層型神経回路網出力の検定法, J. Comput. Chem., Jpn., 投稿中(2005.10.13).
[12] J. Kambe, Y. Yuan, T. Aoyama, U. Nagashima, Neural network analysis of water pollution for a main river, Tamagawa, in Tokyo metropolis, Proceedings of International Conference on Control Automation and Systems 2004, TP04-4 (2004).
[13] J. Kambe, T. Aoyama, U. Nagashima, Comparison with Water Quality of main Rivers in the world, based on OECD reports, Proceedings of International Conference on Control Automation and Systems 2005, FE10-4 (2005).
[14] 大島康行,浅島誠,高橋正征,原沢英夫,松本忠夫編, 理科年表環境編, 丸善, 東京 (2003), ISBN4-621-07335-4.
国立天文台編, 理科年表環境編 第2版, 丸善, 東京 (2006.1.182003), ISBN4-621-07641-8.
[15] T. Sekiya, in T. Fujita, ed., Structure-Activity Relationship and Drug Design, Kagaudojin, Kyoto (1986).
[16] 市川新, 都市河川の環境科学, 培風館, 東京 (1980), 3043-4140-6955.
[17] Environmental Performance and Information Division, OECD Environment Directorate, "OECD Environmental data, compendium 2002".

A 付録 外れ値

A. 1 教師データの外れ値の検出

神経回路網計算において,外れ値検出は関数適合機能と関係する.第2層sigmoid関数,第3層線形関数の場合,神経回路網の機能Fは,

なる有限sigmoid関数展開である.ffは式(8)に同じであるが,パラメータを明示すると,

である.パラメータaiはBP学習により第1,2層間のニューロン間結合から自動的に算出される.qiは第1層のバイアス・ニューロンと第1,2層間のニューロン間結合から同じく自動的に算出される.式(1)のB項は第2層のバイアス・ニューロンと第2,3層間のニューロン間結合から自動的に計算される.式(22)のFを使い任意の2ベクトルX={x1, x2,..., xN}, Y={y1, y2,..., yN}を対応づけると,

である.e={e1, e2,.... eN}において各要素は平均値average(e)=mを中心とした正規分布,

をなすと仮定する.標本分散sを,

とする(回帰分析なのでN-2).eiの検定統計量Tiを,

とする.

なる関数は,中心mからn離れたeiが正規分布全要素に対し比率1-2P(n)=P'(n)であることを示す.ゆえにP'(Ti)の小さいeiが外れ値の候補である.普通はP'(Ti) <0.01(or 0.05)なるei要素を外れ値とする.この基準で(観測値の)外れ値を除外して神経回路網を使用すれば本論文の議論は外れ値が存在する場合も適用可能と考える.ただし本論文の乱数から作成したデータには外れ値は無いので本節の処理を施す必要はない.M個の多変量解析の場合式(A3)は,

となる.
式(A4)〜(A7)の処理は同じである.正規分布でなくとも対称分布ならば上記の考え方はそのまま踏襲できる.

A. 2 Leave-one-out法における神経回路網の関数機能の非同一性の検出

Leave-one-out法においては学習データ集合から1データを除外する.するとそのデータ集合に外れ値が無くとも,また同一の初期値をニューロン間の結合に採用し,同一の回数学習を行っても,神経回路網の関数機能が全て同じとはならない.それをどのように調べるのか検討する.
学習データ集合のm(<=N)要素を除外したとする.神経回路網に入力する説明変数(複数)を{x}とする.出力を{y}とする.この学習データから生成された関数をFm(x)とする.

こうしてN-1要素のベクトル{ym, j}を得る.1データを除去しない場合の出力ベクトル(N要素)を{Y0j}と書く.神経回路網の関数機能に関係する標本分散に相当する次の量を得る.

smもまたN要素のベクトルである.従って,

となるから,正規分布と仮定した場合のその値の正規分布と仮定した場合のTm値で示された関数機能の「位置」P'(Tm)を求めうる.
m番目のデータを除外した場合,著しく関数機能が他と違った場合は,その値が「外れ値」である可能性が高い.

A. 3 外れ値を除外しない場合

A.1-2の方法で検出した外れ値は,通常現象とは無関係な異常値として解析から除外する.
河川水質の予測方法はいくつか知られている.その中に多変量解析を使い測定データ(training data)を学習して,別の測定データ(checking data)の再現性を評価し,水質を決定している要因を取捨選択するGMDH (Group Modeling of Data Handling)法[16]がある.この方法は本論文のように少ない観測値から予測する目的で開発された方法である.このGMDH法でも外れ値の所は予測できないとしている.
我々は一般的な場合,外れ値除外は妥当と考えている.しかし,重金属イオン濃度の水質解析では除外できない場合がある.(欧州第二の国際河川,ドナウ川の支流の1994-5年頃のCd2+, Cu2+, Pb2+イオン濃度値などがその場合である可能性がある[17].)そのとき外れ値をダミー変数により支配方程式の中に取り入れる方法がある.式(11)に取入れた場合,

である.ここでKは観測点Xi-1とXiの間の河川水浄化係数(負が浄化方向,正が汚染方向)である.Lは堰の浄化係数,hiは堰の存在/非存在を示すダミー変数である.Mは未知の新たにダミー変数として取入れた仮定の「現象」miの係数である.{mi}は外れ値の所だけが1,他が0のベクトルである.制約条件は0< siである.
式(A14)を適用して観測値と計算値の差の標準偏差値が急激に0に接近した場合は,ダミー変数に具体的意味があると考え,変数に該当する事件が起こっていないか調査する必要がある.事件が無ければ,外れ値と考え除外するのが妥当であるが,もし「ある」のなら,それは外れ値を利用した観測データ不足問題の一解決法となる.


Return