欠測データを含むデータの解析が可能な多階層型ニューラルネットワークシミュレーション(CQSAR)法を用いた河川の上流および中・下流を示す水質パラメータの抽出
− 東京多摩川の水質データ(1994〜2002)を用いて −

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


Return

1 はじめに

近年,様々な環境保全問題が注目されている.特に河川の水質は,人の生活に密着する重要な問題である.先に我々は第1報 [1]において,1997年から1999年の3年間の17地点, 12種類の水質データを,主成分分析とニューラルネットワークシュミレータNeco[2 - 7]に実装された3層パーセプトロン型ニューラルネットワークの入力パラメータの偏微分係数解析[8, 9]を用いて解析した.17地点の場所とサンプルの分析機関をTable 1に示した.また,17のサンプル採取地点の海からの距離を大まかに捉えるために,測定地点の概略をFigure 1に示した.12種類の水質データの内容はTable 2に示した.


Figure 1. Schematic Map of Sampling Points along the Tamagawa River, Tokyo, Japan. See Table 1 for details of sampling points.

Table 1. Sampling points and agency of analysis
Sampling PointsAgency
1和田橋東京都
2調布橋建設省
3羽村堰東京都
4永田橋建設省
5拝島原水補給点東京都
6拝島橋建設省
7日野橋建設省
8関戸橋建設省
9是政橋建設省

Sampling PointsAgency
10多摩川原橋建設省
11多摩水道橋建設省
12砧下取水点東京都
13第三京浜多摩川橋建設省
14調布取水点東京都
15田園調布堰上建設省
16六郷橋建設省
17大師橋建設省

Table 2. Abbreviations and explanation of chemical index of water analysis
測定項目略称備考
上流・中流・下流の指標pos.上流:1,中流:2,下流:3.
水素イオン濃度pH海水や石灰岩の影響を受けてpHが高くなる場合がある.
溶存酸素量DO量が少ない方が河川は汚染されている.
生物化学的酸素要求量BOD好気性微生物の有機養分分解の際に消費する酸素消費量.
化学的酸素要求量COD酸素を消費する物質の含有量を示す指標.
全窒素T-N生物の生育にとって欠くことのできない栄養塩分.
全リンT-P生物の生育にとって欠くことのできない栄養塩分.
塩化物イオンCl-海水の影響を受けない順流域の人為的汚染状況の指標.
アンモニア性窒素NH4-Nし尿による汚染にはアンモニウムイオンが検出される.
亜硝酸性窒素NO2-N硝酸イオンの還元あるいはアンモニアの酸化によって生ずる.
硝酸性窒素NO3-N一般の淡水には0.1〜1mg/l以上含まれる.
リン酸性リンPO4-P生活排水の影響を受けた河川のリン濃度は大きい.
電気伝導率COND清澄な河川では低く,汚濁の進んだ場合逆に高い値を示す.

主成分分析により得られた第一主成分と第二主成分を用いてサンプルの採取地点をKahunen-Loeveプロット(K-L Plot)すると,その平面上で採取地点を上流,中流および下流域に分類することができた.第一主成分では上流と中流・下流域の2つのグループに分けられる.第二主成分では上流・中流と下流域の2つのグループに分けられる.第一主成分の分離が第二主成分による分離に比べ,より鮮明である.このことより多摩川では上流から中流域の水質の変化に比べて中流から下流域の水質の変化の方が緩やかであることが言えた.そのため多摩川の水質保全には,中流域での水質汚染原因を取り除くことが有効であることが判った.
このように水質に関する化学的パラメータデータを用いた多変量統計解析的手法を用いることで,河川の環境保全の指針が得られることが示唆された.しかしながら,第一主成分と第二主成分は,12種類のパラメータの複雑な線形結合となっており,主成分分析では,12種類のパラメータのうち上流・中流・下流の分類に大きな寄与をするパラメータの抽出はできなかった.
偏微分係数解析からはCl-, COND, NH4-Nなどの微係数が大きいことが判った.これらの変数は主に上流・中流と下流を分類するものである.また,これらに続く絶対値の大きさを持つT-P, NO3-N, PO4-Pは,中流と下流および上流を分ける変数となっており,上流を分類するパラメータは,pHとDOであることが判った.
本研究では,東京都環境局の「公共用水域および地下水の水質測定結果」[10]に掲載されている1994年から2002年までの多摩川の水質データをすべて用いた解析をおこなった.残念なことに2000年から2002年までのデータには欠測があり,従来の方法では十分な解析が困難であった.
そこで青山,長嶋によって新たに開発された,多階層型ニューラルネットワーク(複数)で欠測データを含む学習を行い,その因果関係を解くアルゴリズム(以下,CQSAR)[11, 12]を適用し, 12種類の水質パラメータから,上流および中・下流の分類に重要な役割を果たすパラメータの抽出を行い,そのパラメータの解析を試みた.具体的には12種類の水質パラメータを説明変数とし,上流・中流・下流(pos.)を目的変数として, 3層の階層型ニューラルネットワークのネットワークの重みを解析することで,上流および中・下流の分類に重要なパラメータの抽出し,そのパラメータの経年変化や位置的変化を解析する.

2 データ

1994年から1999年までの多摩川のデータは17地点 (Figure 1),12種類のパラメータが完全に揃っているが, 2000年から2002年までのデータには欠測がある.
Table 3 に1994年から1996年の3年分のデータを示した. Table 4に1997年から1999年の3年分のデータを, Table 5に2000年から2002年のデータをそれぞれ示した. Table 5の陰の部分は, 欠測であることを示している.

Table 3. Values of chemical index of water analysis from 1994 to 1996

Table 4. Values of chemical index of water analysis from1997 to 1999

Table 5. Values of chemical index of water analysis from 2000 to 2002

Shaded squares indicate defective data.

3 欠測データを含むデータの解析法(CQSAR)

階層型ニューラルネットワークを用いた多変量解析は,因果関係が不明で説明変数と観測値が非線形関係にあるような現象の分析に用いられる.化学研究における主な用途は薬理活性解析,環境問題解析などであるが,それらには,しばしば欠測問題が存在する.多変量解析では欠測は単に説明変数の減少に止まらず,複数の説明変数の中の一部分が欠測であっても,それが他の説明変数の該当部分を使用できなくする悪影響を派生する.そのために多変量解析が実施不可能になる場合がある.
多変量解析において欠測を取り扱う方法はいくつか知られている.
方法1:欠測のあるデータを除く.
これは,一番簡単で広く使われている方法であるが,欠測以外の情報を捨ててしまうことになる.情報自体が少ない場合,情報の損失の影響は大きい.
方法2:欠測データに他のデータの平均値を当てはめる.
方法1を用いた場合の最大の欠点である「情報の損失」は回避されるが,欠測データとして本来あるべきデータと当てはめられる平均値の差が大きい場合,解析結果に悪影響を与える場合が多い.
方法3:欠測データを何らかの方法で予測し,それを欠測データの代わりに用いて解析に用いる.
この方法の解析結果は欠測データの予測精度に依存する.通常は方法2よりよい結果を与える.この方法は欠測データ予測精度が解析結果に大きく反映する.
方法4:欠測を持つデータに対する統計モデルを仮定して,欠測データ周辺の類似性から欠測データを予測する.
この方法としては,EM(Expectation Maximization)法[13, 14]が良く知られている.EM法は不完全データの問題を完全データの枠組みで逐次的に解く方法であり,その操作は以下の通りである.
(1) 手元にある観測データYobsだけでは解を導くことが困難な問題に対して,解決が容易なレベルまでデータを完全化し,(欠損データ部分Ymisが存在するとして)問題を定式化する.(2) パラメータの推定値Iを導くことが目的であればYmisに暫定的な値を埋め込みI0を求める.(3) I0を用いてYmis を改良する.(4) Inの値が収束するまで前の2つのステップを繰り返す.
ここで我々の注目する化学,薬理学,環境科学データについてYmisの存在は確かである.そのため,操作(1),(2)は可能である.しかし(3),(4)は問題を多く含む.たとえば,そのような結果を繰り込んで観測値を補正することは物理,化学の分野では禁止されている.このような禁制は薬理活性問題の化合物の置換基物性値と生理学的活性の測定の関連を考えれば当然である.逆に統計学では「確率論的モデルを仮定し,その現れた現象を観測データ」と考えている.そのような場合には,(3),(4)は適用可能である.言い換えると化学,薬理学,環境科学データについては,ナイーブなEM法の適用はもともと禁じられている.
つまり薬理活性問題や環境問題では,現象に関係がある複数種類の測定値を用い(仮定),それらを説明変数として現象の再現を試みる.この際,関係の無い説明変数の適用,過剰のパラメータの採用,不適切な関数などの存在など手段が適切でないならば測定値は高精度で再現されても,現象について得られる情報は少ないか全く無いということは広く知られている.そのため,このようなシミュレーションの有意性判定は現象の総合的な知識が必要である.そのため,EM法をそのまま薬理活性問題,環境問題に適用することは難しい.
CQSAR[11, 12]は,方法3つまり欠測データを何らかの方法で予測し,それを欠測データの代わりに用いて解析に用いる方法に基づくアルゴリズムであり,実測の50%以上にも及ぶ欠測の悪影響を取り除き,多変量解析を可能にする階層型ニューラルネットワークを用いた新しい手法である.以下では方法1とその精度を比較する.

3. 1 CQSARのアルゴリズム

CQSARのアルゴリズムは,大まかに2つのステップから構成される.CQSARの第1ステップは,欠測部を各説明変数および観測値を独立に補完することである.この時,補完法の選択は主に観測数に依存し,観測値の分布と連続的変化の有無に依存する.我々の経験から,観測数が5個未満ならば線形最小2乗法を用いる.観測数が5個以上7個未満で現象の非線形性が明らかな場合,例えば現象が単調で次第に飽和するような場合ならば,単一sigmoid関数,また,極値の存在が明らかで,その部分の確定が重要な薬理学分野の問題の場合は,可変絶対値関数fittingで補完する.ただし環境科学問題のように過去の測定データが利用できる場合は,その測定データを関数の離散表現として新たにパラメータを付加した関数fittingで補完する.観測値の数が7以上で説明変数の変化が連続的である場合,ニューラルネットワークと等差数列を使う補完法が良い.この補完法は現象の支配方程式を要求しないので容易に使用でき,さらに精度も優れている.
第2ステップは第1ステップで完全化したデータのニューラルネットワークを用いた定量的構造活性相関QSARである.この際,補完データには誤差が相当量含まれている可能性が高いので,その悪影響を除くため,通常の学習に用いられるバックプロパゲーション学習法と再構築学習法[15]などのsimulated annealing法の併用が望ましい.
CQSARの適用制限条件は「観測数がある程度(経験的には7以上)無いと精度が保証できない」ことである.またその偏微分係数の精度を要求するのならば,経験的に11以上の観測数が必要である.したがって広範囲の問題に適用可能である.

3. 2 CQSARの精度

簡単な例を用いて,CQSARの精度を示す.CQSARのみならず計算法の適用限界を数値的に調べる場合,誤差のない理想的なデータが必要である.環境および化学分野で一般的な測定に類似の初等関数の等間隔サンプリングにより,その種のベクトルを生成した.説明変数には,x1/2, x, x2, (x - a)2, (x - a)(x - b)(x - c)の五種類の関数値を与え,観測値にb - (x - a)2 の関数を選んだ. ここで0 < x, a, b, c < 1である.ニューラルネットワークでは増加と減少は学習上同等であるので,増加のみを選んだ.サンプリング数は環境および化学分野で頻繁に見出される測定値数{18}とした.
欠測のパターンをFigure 2に示す.パターンAは,前半の6データに欠測が存在する.パターンBは中央部の6データ,パターンCは後半の6データに欠測が存在する.パターンDはすべてのデータに欠測が存在する.パターンDの欠測データの割合は17%である.パターンEは,すべてのデータに欠測が存在し,50%のデータが失われている.パターンDおよびEでは,欠測データがある場合に広く用いられている方法1:欠測のあるデータを削除する方法は適用できない.


Figure 2. Patterns of defects.

用いた階層型ニューラルネットワークの構造は3層型であり,ニューロン数はバイアスニューロンを含み,それぞれ6, 4, 1である. 中間層のニューロン数4は,欠損のないデータを再構築学習法を用いて学習させ,10から減少させることによって得た.中間層のユニットを3, 4, 5としたニューラルネットワークの情報量基準AIC[16]はそれぞれ, 95.2, 73.3, 92.1 であり.この場合4が適当であることが判る[17].
本例でのCQSARによる欠測データの補完は,データ数が7以上であるので,ニューラルネットワークと等差数列を使う補完法[11]を用いておこなった.
leave-one-out 法およびCQSARで予測された値を学習したニューラルネットワークのAICを計算し学習の検定を行った.leave-one-out 法での誤差の平均標準偏差とAICをTable 6に方法1の結果とともに示した.誤差の平均標準偏差をみると,パターンAでは,CQSARは方法1に比べ10倍程度性能が良いことが判る.パターンB, Cでもわずかであるが良い.CQSARも方法1も,パターンBに比べA,Cが悪いのは,それらが外挿に対応しているからである.パターンD, Eに関して,CQSARがパターンAとほぼ同等の精度を持っていることは驚きである.言うまでもなくパターンD, Eに対し,方法1は適用できない.

Table 6. Averaged standard deviation of error(ASDE) and AIC for each defects pattern
Pattern of defectsCQSARApproach1:remove defects
ASDE (ratio*)AIC(ratio*)ASDEAIC
A: defects in the first block0.0151(8.6)74.2(1.18)0.130187.6
B: defects in the center0.0046(1.4)67.3(1.06)0.006371.6
C: defects in the terminal0.0343(3.2)75.5(1.11)0.110784.2
D: dispersed defects of 17%0.0132(¥)71.1(¥)- impossible -
E: dispersed defects of 50%0.0195(¥)74.4(¥)- impossible -
* ratio=Approach1/CQSAR

CQSARで補完を実施したニューラルネットのAICは,方法1よりも小さな値を示しており,より妥当な学習をしていることがうかがえる.また,AICの値が小さいほど誤差の平均標準偏差が小さく,両者には何らかの相関があるように見える.

4 CQSARの多摩川の水質データへの適用

欠測データを含む多摩川の水質データについて,欠測を含むデータを除く方法1とCQSARの適用可能性を比較した.Table 7に各年度のデータの欠測の割合と適応可能性を示した.Table 7のindexesとobservation pointsの意味を説明するためにデータの構造の概要をFigure 3に示した.

Table 7. Ratios of defect parts and scopes of two methods
yearsdefectsindexes/observation pointsapproach-1CQSAR
19940.00000.0000/0.0000practicablepracticable
19950.00000.0000/0.0000
19960.00000.0000/0.0000
19970.00000.0000/0.0000
19980.00000.0000/0.0000
19990.00000.0000/0.0000
20000.02450.2941/0.0833
20010.14710.7059/1.0000impossible
20020.20100.7059/1.0000


Figure 3. Structure of data.

1994年から2000年までのデータではすべてのデータが存在するため方法1でも学習可能である.ところが2001年と2002年では垂直方向にすべてのデータを除かなければならないので方法1の適用は不可能であった.一方,CQSARでは,2000年から2002年においても欠測データを補完するため,学習可能であった.

5 多摩川のデータでの上・中・下流の3つの分類の検定

Tables 3, - 5に示した測定データを見れば下流の16-17地点のCl-とCONDが極端に大きいことから,中流と下流の分類を与えるパラメータはCl-とCONDであることが直ちに判る.ところが第1-6地点を上流,7-15を中流とする分類を与えるパラメータは明白ではない.そのため,CQSARを用いて,1)第1-7地点を上流,8-15を中流,16-17を下流とした場合,2)第1-6地点を上流,7-15を中流,16-17を下流とした場合,3) 第1-5地点を上流,6-15を中流,16-17を下流とした場合の3つの場合について学習の誤差およびAICを指標にその分類の妥当性を評価した.
それぞれのネットワーク構造は再構築学習法により最適化し,AICを用いて確認した.この時,再構築学習は10万回行い,50回ごとに0.05以下の結合を消去した.消去はそれぞれの場合,合計2000回行われネットワーク構造の最適化が行われた.また最適化された中間層ニューロン数に1を足した(または引いた)ネットワークで学習を行い,AICを計算して中間層ニューロン数の妥当性を検証した.例えば2)の場合,再構築学習法では中間層ニューロン数が2となったので,中間層のニューロン数を,1, 2, 3として学習を行いそれぞれのAICを計算すると, 66.2, 28.2, 32.6となり,中間層ニューロン数2が妥当であることが判った.
最終的なネットワーク構造は,1)の場合が(12, 3, 1), 2), 3)の場合が(12, 2, 1)であった.
それぞれの場合の学習誤差の最大と平均およびAICをTable 8に示す.1-6を上流とする場合が,一番エラーの少ない学習をしており.AICも一番小さい.この結果よりTables 3, - 5にあるように,第1-6地点を上流,7-15を中流,16-17を下流とすることは妥当であることが判る.

Table 8. Conformation of assignment to upper region
Assignment of upper regionMax. erroraverage errorAIC
1)      1-70.07700.0070434.2
2)      1-60.00260.0009828.2
3)      1-50.00440.0011128.5

6 CQSARによる学習と再構築学習法を用いた結果

先の偏微分係数解析を用いた解析[1]は,上流・中流・下流を分ける明確なパラメータ抽出はできなかった. CQSARと再構築学習法を用いて学習を行うと,ネットワーク構造が最小化され中間層ニューロン数が2で十分であることが判った.通常のバックプロパゲーション学習法では5ニューロン以下では学習が上手くいかず,6ニューロン必要であった.これが前回の論文[1]で偏微分係数で特定のパラメータが抽出できなかった理由である.言い換えれば,本例の場合は,適切な学習には再構築学習法のようなシミュレーテイッドアニーリング法とバックプロパゲーション法の組み合わせが必要であるということである.
最適化されたニューラルネットワークの結合行列をTable 9に示した.再構築学習法による結合の淘汰により,上流と中・下流を分ける主なパラメータとして,DO(溶存酸素量)とT-P(全リン)およびPO4-P(リン酸性リン)が抽出された.ただしPO4-Pの影響は小さい.以下DOとT-Pについて考察する.

Table 9. Weight matrix of optimized neural network by reconstruction learning
The first layer - the second layer
pHDOBODCODT-NT-PCl-NH4-NNO2-NNO3-NPO4-PCOND
10-1800035000020
2022000-380000-50

The second layer - the third layer
12
14350

Figure 4に1997年から2001年までの各測定地点でのDOの変化を示す.ここでは,データ中の最大値(第5地点,1997年)を1.0,最小値(第17地点,1997年)を0.0としてスケールしてある.DOは値が大きいほど清浄であることを示す.これを見ると確かに第6地点と第7地点の間にターニングポイントがあることがわかる.第5地点を見ると1997年で非常に清浄であった水質が2001年にかけて徐々に汚染が進んでいることが判る.反対に第10地点は清浄さが回復している.第3地点は1997年から清浄化が進んだが,2001年にはまた1997年レベルまで汚染が元に戻ってしまった.逆に第8地点では,1997年から汚染が進み,2001年ではまた1997年レベルまで清浄化が回復した.


Figure 4. The changes of DO-index at 17 observation points between 1997 and 2001.

この第3地点と第8地点,第10地点を除くと1997年から2001年にかけてほぼすべての流域において水質汚染が進んでいることが判った.
また各観測地点のT-P変化をFigure 5に示す.T-Pは値が大きいほど汚染が進んでいることを示している.ここでも最大値と最小値を1.0, 0.0としてスケールしてある.これを見るとDOと同様第6地点と第7地点の間にターニングポイントがあることが判る.


Figure 5. The changes of T-P-index at 17 observation points between 1997 and 2001.

ここで判ることはでT-Pは上流すなわち第1地点から第6地点に関しては全く情報を与えないが,中・下流域で値が大きくなっている.中流域および下流域のT-Pの変化は,1997年から1999年にかけて汚染が進み,その後汚染は減少し,2001年では1997年位までは回復していることを示している.
新聞の記事 (2004年3月21日朝日新聞)によると「羽村堰の流水効果で多摩川の水質は年々良くなっている.」というコメントは,T-Pを見る限り正しいようであるが,DOではそのような傾向は見られない.今後堰の清浄化の効果を見積もるための方法の開発が必要である.
水質汚染という観点からみると,多摩川では上流から中流域の水質の変化に比べて中流から下流域の水質変化の方が緩やかであることが判った.このことにより多摩川の水質保全には,中流域での汚染原因を取り除くことが重要であることが示唆された.

7 まとめ

欠測データがあり,十分な解析が困難であった多摩川の水質データに,欠測データを推測して補完する方法であるCQSARを適用して解析を行った.ニューラルネットワークによる多摩川にデータでの上流および中・下流の分類の検定で,第6地点を上流と中流の境界とすべきであることが判った.
CQSARでは,欠測データを補完し,学習するため,欠測を持つ2000年から2002年のデータ解析も可能 となった.
本研究により水質検査項目の12種類のパラメータのうち,主なパラメータとしてDOとT-Pが抽出された.
DOはすべての流域において水質汚染が進んでいることを示した.T-Pは上流に関して全く情報を与えないが,中・下流で浄化傾向にあることを示した.T-Pの結果は朝日新聞の記事の記述にある「羽村堰の流水効果で多摩川の水質は年々良くなっている.」を一部支持した.
今後堰の効果を見積もるための河川モデルの構築が必要である.
水質汚染という観点からみると,多摩川では上流から中流域の水質の変化に比べて中流から下流域の水質変化の方が緩やかであることが判った.このことにより多摩川の水質保全には,中流域での汚染原因を取り除くことが重要であることが示唆された.

参考文献

[ 1] 神部順子, 福田朋子, 長嶋雲兵, 青山智夫, J. Chem. Software, 8, 27 (2002).
[ 2] 井須芳美, 長嶋雲兵, 細矢治夫, 青山智夫, J. Chem. Software, 2, 76 (1994).
[ 3] 井須芳美, 長嶋雲兵, 細矢治夫, 大島茂, 坂本曜子, 青山智夫, J. Chem. Software, 3, 1 (1996).
[ 4] Isu, Y., Nagashima, U., Aoyama,T., Hosoya, H., J. Chem. Info. Comp. Sci., 36, 286 (1996).
[ 5] 藤谷康子, 小野寺光永, 井須芳美, 長嶋雲兵, 細矢治夫, 青山智夫, J. Chem. Software, 4, 19 (1998).
[ 6] 田島澄恵, 松本高利, 田辺和俊, 長嶋雲兵, 細矢治夫, 青山智夫, J. Chem. Software, 6, 115 (2000).
[ 7] 福田朋子, 田島澄恵, 斎藤久登, 長嶋雲兵, 細矢治夫, 青山智夫, J. Chem. Software, 7, 115 (2001).
[ 8] Aoyama, T., Ichikawa, H., Chem. Pharm. Bull., 39, 372 (1991).
[ 9] Aoyama, T., Wang, Q., Zhu, H., Nagashima U., IPSJ-SIG notes, 2000-HPC-84, 7 (2000).
[10] 公共用水域および地下水の水質測定結果,東京都環境局(1994-2002年度).
[11] Bitou, K., Yuan, Y., Aoyama, T., Nagashima, U., Proceedings of International Conference on Control Automation and Systems 2003, CD-ROM (TA06-04),10.22-25, 2003.
[12] 青山智夫, 神部順子, 長嶋雲兵, J. Comp. Chem. Jpn (2006), 投稿中.
[13] 渡辺美智子,山口和範, 「EMアルゴリズムと不完全データの諸問題」, 多賀出版, ISBN4-8115-5701-8.
[14] 渡辺美智子,Web資料:
http://stat.eco.toyo.ac.jp/~michiko/em/emohp(0)/emohp.ppt
[15] Aoyama, T., Ichikawa, H., Chem. Pharm. Bull., 39, 1222 (1991).
[16] Akaike, H., IEEE Trans. on Automatic Control, AC-19, 716 (1974).
[17] 栗田多喜夫, 電子情報通信学会論文誌, J73-D-II, 1872 (1990).


Return