ノート
線形変換と交差検証法による空間線量率予測
2017 年 16 巻 3 号 p. 70-76
詳細
2017 年 16 巻 3 号 p. 70-76
環境中のγ線線量値を長期にわたり予測するため，放射性核種の物性と担体の拡散を考慮した空間線量値の時間変化を直線化する方法を提示した．現象の変化が直線化できると予測は線形補外となる．交差検証法で過去のデータの一部を用いて予測を行い実測との対応から予測精度を見積る．将来の予測計算も，過去の処理と同一である．方法自体の精度のオーダは同じである．本処理の可能性を示す指標も示した．
福島県伊達市下小国集会所のモニタリング・ポストのデータでテストすると，線量予測式と実測値の決定係数は0.9以上，1年先の予測値では0.7であった．観測点の気象的擾乱，周辺環境の変化は決定係数の低下として現れる．
我々は，福島原子力発電所事故により他の市町村に避難し，帰宅時期を考慮している人を想定した．
福島第一原発の事故により放射性物質が東日本各地に飛散し，降雨などにより沈着した [1]．原子力規制委員会(Nuclear Regulation Authority: NRA)は2012年4月頃より各地に順次γ線量用のモニタリング・ポストを設置し，測定線量値をデータベースとして公開している [2]．測定精度は1 [nSv/h]，10分間隔の測定である．一般の者が入手できる最も精度の高い環境放射線データである．蓄積データ量は各地点で5年分弱，約24万データである．このデータから将来の線量値を予測できれば，地域の将来の計画等に有用と考える．
空間線量率はモニタリング・ポスト(Monitoring Post: MP)の検出器を通過する300∼800 [keV]領域のガンマ線束(flux)の仕事率[(J/kg)/s] = [Gy/s]である．この値を人体のモデルを決めて，検出器のファームウェアでGy/Sv変換を行い生物学的影響の大きさである線量当量[Sv/h]にする．その値をNRAが公表している．検出器の地表面からの距離は現在，大部分のMPで1 [m]である．違う場合はその旨表示されている．MPの設置位置の自然放射線の寄与は除かれている．
ガンマ線束は検出器を中心とした球殻内の放射性核から発する．ガンマ線の半減到達距離は核種毎に違う．IAEAは134Cs, 137Csについて71.4, 69.7 [m]としている [3]．ゆえに検出器近傍の放射性核の分布が空間線量値に寄与する割合が大きい．ガンマ線が検出器に仕事をするのと同様に周辺の物質にも仕事をする．ガンマ線の透過率という．IAEAは半値層(厚)として材質ごとに値を公表している [3]．環境中で注目すべき物質は水である．水は軽い物質であるがガンマ線を減衰させる．ゆえにMPの近くの地面が積雪状態であると空間線量率は一時的に低くなり融雪後元に戻る．放射性物質は移動していない．このような測定データの取り扱いは注意が必要で，考慮しないと予測に誤差をもたらす．
検出器周辺の放射性核の分布からガンマ線束を算出するのはガンマ線のモンテカルロ・トレーサ(シュミレータ)を使う [4]，そのとき計算域内の遮蔽物の物性値と位置情報が必要で，空間内の諸物のモデリングが必要である．そのような多種類の空間データを必要としない，空間線量率の経時変化を追跡できる実験式への要求は存在する．放射性物質の質量濃度変化M (t)に関する実験式は幾つか公表されている．Garland, Pomeroy [5]とHatano [6]は，
M(t) = M0exp(-λt), (1)
M(t) = M0t-4/3exp(-λt), (2) を提案し，Kinase [7]は，
M(t) = M0{fexp(α1t)+(1-f)exp(α2t)} ×{x exp(λ134t)+exp(λ137t)}/(x+1), (3) を提案している．ここでM0は初期の沈着濃度[Bq/kg]，λは対象となる1個の放射性核の半減期, tは時間である．fは早い環境中の減衰線分比(0<f<1)，α1はその環境半減期である．α2は遅い環境中の成分の半減期である．λ134, λ137は134Cs,137Csの半減期である [8]．xは134Cs,137Csの存在比である．福島原発事故後の調査では放射性物質濃度の減衰には2つの要素が認められる [9]ので，現時点では式(3)が精度が良い．しかしf, α1, α2は不確実性のあるパラメータである．これを克服するため区間データとしモンテカルロ・シミュレーションで空間線量率の取りうる範囲を計算する．ゆえに，測定結果は決定係数(R2)ではなく，「実測値のほんどがパラメータの不確かさを考慮した推定値の変化の範囲内にある」という結果である [7]．我々は利便性を考慮して，区間データではなく単一値とその信頼性値を求める方向を目的とする．
以上は過去の放射線量に対するfittingに関する検討である．Fittingの式を予測に利用することは可能である．その予測時のR2は統計学の標準的処理 [10]から決定したいと我々は考える．すなわち，現象の線形化処理と交差検証法(leave-one out法)を採用する．線形化を採用する理由は回帰線とR2の定義である [11]．線形化可能な範囲は現象の非線形性に依存する．第一段階としてそれを決定する．
核種が環境中に放出された日を起点0 [d]とする．本論文では2011年3月15日0時(JST)とする．起点からの経過日をt [d]とする．両核種の合計の沈着量をM0とする．両核種の存在量M (t)は，
M(t) = M0{xexp(λ134t)+exp(λ137t)}/(1+x), (4) である．単位は[Bq/volume]である．VolumeはSI単位系では[m3]である．放射性Csの特徴として土壌中の粘土層と強く結合するので，今回の事故のように降雨等で沈着した時は地表面に数[cm]の層状に分布する．そのためIAEAは土壌汚染ではvolumeの単位を[m2]としている [3]．我々もこの慣行に従う．ゆえに単位は[Bq/m2]とする．134Csと137Csの無限平面均一分布から放出されるγ線エネルギーの強さは5.4: 2.1×10−6 [mSv/h]/[kBq/m2]である [3]．ゆえに線束量E (t)は，
E(t) = KM0{5.4x exp(λ134t)+2.1exp(λ137t)}/(1+x), (5) である．単位は[Sv/h]である．Kは目的地点のγ線遮蔽環境を表すスカラー定数である(物理的背景は付録A)．
放射性核種は土壌中に存在し，移流拡散方程式に従い移動する．その速度を決定する移流拡散係数は定数ではなく時間，土壌の物性，気象，地下水の影響下にある．それらを全て考慮するのは現実的ではない．そこで放射性核種を動かす媒体の運動が小さいことを利用して移流場の小さい拡散方程式の一般解を考察する．拡散方程式は変数分離形で解け,
M(t,x,y,z) = ΣLM0(L)GL(x,y,z; kx,ky,kz)exp(-λt), (6) である．添字Lは初期状態の空間的境界条件を1つのパラメータで表現したものである．{kx,ky,kz,λ}定数は拡散係数に関係する量である．拡散係数が土壌中の位置(L)で異なるのならkxはkx (L)である．式(6)で注目すべきは時間方向にはexp (-λt)となることである．ただし，
ΣLM0(L) GL(x,y,z; kx,ky,kz)が付随しているから単純な指数関数ではないが，近似として式(1)が成立する場合もあり得る．これが初期の研究の正当性の根拠である．
式(5)を放射性核の含有する媒体の拡散を考慮して拡張する．Σ演算子と複数の指数関数の和の形式にすると，パラメータ増によるfitting能力を向上させる方向である．予測を目指すのならOccamの剃刀の原理 [12]に反する．Σ演算子とパラメータを明示的に使用しないのなら，
E(t) = KM0exp{β(t)t}{5.4xexp(λ134t)+2.1exp(λ137t)}/(1+x), (7) である．ここでβ(t)は放射性核種の物理的崩壊式を環境中の線束量に合致させる補正関数として式(7)に導入されている．β(t)は空間線量の実測から決定する．ゆえに，ベクトルであり，決定法は以下である．
式(7)を線量値obs (t)に適用し差分D (t)を得る．
D(t) = obs(t)-E(t), (8)
2つの時刻t = {t0,t1}においてD (t0),D (t1)を得る．D (t0) = 0となるようにKを設定してD (t1)の計算にそれを使用する．式(7)から，
D(t1) = exp{β(t1)t1}, β(t1) = Ln{D(t1)}/t1. (9)
式(9)はβ(t)関数がベクトルとして得られることを示している．ただし，分母に時間の項があるため短期間の測定では誤差が無視できない．影響を少なくするには数多くの測定と最小二乗法などの適用が必要である．
ここでβ(t)関数の一例を示す．現在，空間線量率の将来予測を必要とする人物は避難先から帰郷するか否かの境界の線量地域に居住されていた方と思う．ゆえに環境基準 [13]の0.23 [μSv/h]前後の地域の例を計算するのが適切と思われる．我々は福島県伊達市下小国集会所{140.579188E, 37.768476N}を選んだ．そこは2012年4月1日には0.63 [μSv/h]であった空間線量値が2017年12月10日には0.20 [μSv/h]になっている．放射性核の崩壊現象だけの計算では0.28 [μSv/h]になる．その差がβ(t)関数の効果である．
式(9)による例をFigure 1に示す．
Vector representation of β(t) function.
Vertical axis is the value of β(t) function in unit [1/d].
Horizontal axis is the elapsed day from March 15, 2011. Dots correspond to β(t) function per 15 [d], and straight line is the regression one.
式(9)から予測される通りβ(t)関数値は分散する．直線回帰式は
β(t) = − 1.439×10−8t − 9.905×10−5, R2 = 0.0458, (10) である．R2値が小さいため，同関数を「直線」としてでは無くt = 「定義域の最大t」の時の値(定数)と見なすべきであろう．ゆえに，β(2190) = − 1.306×10−4．半減期で表すと14.5[y]である．
ここで，有効数字を4∼5桁で計算しているが，元データの有効数字が2桁であるから，計算結果はそれ以上にはならない．しかし多段の計算のため丸め誤差を累積させないように2桁以上余分に計算する．以下，この方針に従う．
Figure 1の点の分布には段差構造は無いので，同地では放射性物質が環境的擾乱により大きく移動拡散した事実は無い．念のために段差や大きな変動がβ(t)関数に現れる地点を探した所，原発の近くでその例を発見したので付録Bに示す．その場所では環境的擾乱が空間線量値を大きく変動させるので予測は難しい．さらに長期間の測定が必要である．
下小国地区では定数としたが，定数となるのは拡散方程式の初期の質量分布が均一か一点である時の一般解の時間項の部分である．質量分布の均一性を保証する半径は放射性核のガンマ線の半減到達距離の2∼4倍である．ゆえに200 m程度である．この範囲の均一性が要求される．
結局，式(7)は，
E(t) = KM0exp(βt){5.4xexp(λ134t)+2.1exp(λ137t)}/(1+x), (11) に簡約化される．福島原発事故の場合x = 1(付録A)なので，1/2，M0をKに繰り込んで，
E (t) = Kexp (βt){5.4exp (λ134t)+2.1exp (λ137t)}としても良い．
将来の空間線量率を予測する場合，関数E (t)が直線に近い方が高精度に予測できる．一般的には自然対数loge(t): = Ln (t)を取り，Ln{E (t)} =
Ln(K)+Ln[{5.4exp(λ134t)+2.1exp(λ137t)}exp(βt)], (12) とする．記号“: = ”の意味は「と定義する」である．この式はtに関して直線的に変化しない．
Ln{E(t)}F(t) = At+B, (13) となるような関数F (t)の性質を調べる．まず定義域を設定する．“t = 0”においてE (0) = 1ゆえにB = 0，Aをt = 4000 = tLでE (tL) = A (tL)となるように決定する．以下tの単位を[d]とする．tLはF (t)が単射の範囲で定義すること．
F (t)の定義域は0 < t ≤ tLである．ゆえにLn{E (t)} < 0, A < 0である．
F(t) = At/[Ln[(1/7.5){5.4exp(λ134t)+2.1exp(λ137t)}exp(βt)]], t≠0. (14)
以下，βを半減期10 [y]に対応させる(他の値は付録Cにある)．βが決まると，式(14)からF (t)の関数形が定まる．式(14)を使いやすいように多項式とすると，7.5はt = 0における規格化のためである．
F(t) = −1.70110×10−12t3 +1.73945×10−8t2 +4.90167×10−5t+6.34876×10−1, R2 = 0.999996, (15) である．代数的逆関数は，
F(t)−1 = 1/F(t) = 3.43738×10−12t3 − 1.52480×10−8t2 −1.38128×10−4t +1.57844, R2 = 0.999990, (16) である．時間変数の多項式であるがF (t)は比であり，単位は無い．このtは事故後の経過日である．
時刻t1(> 0)の観測値をdose (t1)とする．普通は“ t1 ”を観測開始日とする．F (t)の定義域，t1日までの減衰量と測定値から，
Cr = Ln[dose(t1)/Ln[(1/7.5){5.4exp(λ134t1) +2.1exp(λ137t1)}exp(βt1)]], (17) を導入する．観測値dose (t)は，
Linear{dose(t)} = Ln{dose(t)}F(t)+Cr, (18) の処理で直線化できる．定数Crの単位はdose (t)と同じである．式(18)ではバイアス項として働く．
この変換は二段階変換である．観測値を直接，直線化する一段の変換も可能であるが，逆変換が二値をとる関数になる．使用法が制限される(プログラミングではIF文の多階層化構造となる)ので，このような二段階変換を採用した．
このLinear ()関数の世界では測定値の時間変化は式(11)で表される物理過程は直線化されている．線形最小二乗法で傾き(A)とバイアス部(B)を求め，Linear{dose (t)} = At+Bとし，tを将来に延長すれば線量予測が可能である．現象が放射性の核崩壊と拡散だけならばB = 0であるが，観測値に適用するとB≠0となる．式(11)とは異なる物理過程が存在する．
予測値は式(11)のLinear ()関数での予測値である．それから元の観測のdose (t)を求めるには，
dose(t) = [Linear{dose(t)}-Cr]F(t)−1, (19) とする．
原子力委員会の各地のγ線量値データベースは10分毎の測定である．これをdose10(td, tm)と書く．添字dは[d]，mは10分毎のカウンタ値(0∼143)を示す．欠測が存在するのでm添字の数は24×6 = 144と限らないのでM (d)とする．欠測でもカウンタは離散的な時刻に対応するので値が増大する．
dose(t) = Σm dose10(td,tm)/M(d), (20)
~tm = Σtm/M(d). (21)
カウンタ平均値~tmが72 ± 5から外れる場合(偏在欠測データ)は精度的に良くない測定なので欠測とした．式(21)によって平均化した線量値を予測計算に使用する．
ここで欠測について説明する．欠測を，観測機器の不調によるデータが無い場合と，降雪時の二種類とした．降雪時を除く理由は，地上の積雪が土中の放射性Csのγ線の一部を遮蔽し，線量値を低下させるためである．
ここで式(11)以外の物理過程を考察する．土中の水層(又は水分)は積雪と同じくγ線を遮蔽する．ゆえに土中の水分が季節変化すると，その変化は線量値に反映される [14]．反映される度合いは観測地点により異なる．その度合いが大きい程，式(11)に基づく予測精度を低下させる．本論文で取り上げた観測点は2014年以降にその度合いが大きくなった．
以上のように我々は考えていたが，最近の研究によると土中の水分ではなく植物活動に原因があるとされる．近く論文が公表される．植物の代謝は気温，地温に関係する．
測定の期間を単位[d]で[t1,tL]とする．“t1 ”は0.5年以下の半減期の核種が少なくなる365 [d]以上とする．環境中の放射性Csのγ線強度を調べるため，周辺環境の年周期を無視できない．また予測精度を調べられる期間が3回以上は必要として(tL-t1) > 365*n, n = 4とする．予測精度は部分集合から残りの集合の値を検定する交差検証法であるので線形化が成されていればn = 1にこだわる必要はない．
福島県伊達市下小国集会所の近くの森林で2016.3.30∼4.2大規模な山火事があり，福島大学の現地調査ではリター層が灰になり土壌流出が眼視的に考えられる状況であった [15]．福島原発事故由来の放射性Csがリター層に多く存在することは1990年代の報告 [16]と異なるが，最近の報告 [17]にある．山火事が放射性物質を再拡散させることがチェルノブイリ事故で報告されている [18]．この山火事の近隣で0.2 [μSv/h]を超える線量値を示す地点は下小国集会所である(それ以下の低い空間線量率では測定誤差に影響が大きい)．
予測は測定開始日t1: = 2012.3.31から，2012.12.10までのデータから1年先の365 [d]を予測し，測定値と比較した(欠測日と降雪日を除く)．下小国集会所の降雪日は671∼693, 712∼716, 768∼770, 1011∼1020, 1062∼1085, 1417∼1424, 1770∼1787 [d]である．4節の結果から式(15)の係数は付録Cの15 [y]のものに変更した．
この処理を，データを1年分ずつ増やして2016.12.9まで4年分計算した結果をTable 1に示し，1年先の予測精度を決定係数で示した．また現在入手できる2016.12.9までのデータから同様の予測法で今後の3年間の線量値を予測した(Figure 2)．
| Observation Period | Effective**/total days | Fitting R2 | Prediction R2 | |
| 1 | B*~2012.12.10 | 255/255 | 0.872 | 0.759 |
| 2 | B~2013.12.10 | 583/620 | 0.956 | 0.784 |
| 3 | B~2014.12.10 | 875/985 | 0.982 | 0.413 |
| 4 | B~2015.12.10 | 1217/1350 | 0.986 | 0.636 |
| 5 | B~2016.12.9 | 1559/1715 | 0.909 | Figure 1 |
*) B is the start day that is 2012.3.31.
**) Effective days is that no-measurement and snow days are subtracted from total days.
R2 is calculated for the effective days in Linear{dose (t)} set.
Prediction R2 is calculated among predicted dose rates and observations. The observation period is from 2012.3.30 to 2016.12.9.
Prediction and observation curves from 2012 to 2019.
The vertical axis is the ambient dose rate [μSv/h] and the horizontal is elapsed days from March 15, 2011. The blue zigzag are observation dose rates, and the red is fitting (until 2016) and prediction (after 2017).
Table 1より，式(11)は決定係数0.9以上で線量値を追跡できるが，1年間の予測に適用すると同係数は0.7前後である．年度により気象環境が異なるので，その影響のためであろう．複雑な気象環境の影響を受ける環境中の放射線量を0.7前後で予測できることは一つの成果と考える．
Figure 2では2016年末(現在)から3年分の予測曲線を示した．気象環境の予測は考慮されていない．予想曲線は実測データの存在する部分は追跡できているし，スムーズに延長されている．
Table 1の項番4と5を比較すると，2016年3月末の山火事の影響が決定係数の低下に表れている．放射性Csは森林のリター層に存在している．そこが灰化したため式(11)で表現できない土壌中の水分布の構造が変化し決定係数を低くしたと考える．ただし2016年末現在，予測線からの分散は大きくはなったが，全体としては離脱していないので放射性Csを含む土壌流出は空間線量計では検出できないレベルである．
ここで，2014年12月11日∼2015年12月10日の予測R2が0.413と低下している理由を考察する．
Figure 3はFigure 2の2014年8月7日∼2015年7月16日分を拡大表示したものに相当する．その図から考察すると観測値の分散が大きいことと長期間の欠測(2014年2月4日∼4月10日)が特徴である．欠測の前後で観測値が計算値の前後に振れている点は積雪と土壌中の水の影響を想像させる．
Fitting and observations of low R2 period in Table 1.
Blue dots are observations. Red curve is calculated by Eq. (11).
MPの近くの積雪データが公表されている地点は福島気象台である．そこでは2014年2月5日∼3月1日まで積雪があり最大深37 cmである，3月6∼8日2∼3 cm，3月21∼22日2 cmである．積雪が消えてもしばらく線量値は元のレベルに戻らないので，この欠測はやむを得ない．ゆえに有効なデータ数が減少し，土中の水分の影響による分散が大きくなったことが原因である．この気象的外乱は予測不能である．
長期間にわたる環境中のγ線線量値を解析するため，放射性核種の物性と担体の拡散を補正し，線量値の時間変化を直線化する式を提示した．式の単位の対応についても説明した．提示式の処理は線量値を自然対数変換し，それに時間補正乗数を作用させる二段階変換である．補正乗数は11年分の表現形を示した．
現象の時間変化値が直線化できると予測は線形補外となる．過去のデータの一部からの予測は実測と比較対応が可能であり，予測精度も計算できる．将来の予測では実測は無いが，予測計算が過去の処理と同一であるから，方法自体の精度も同じである．予測が外れるのは降雪などの気象要因による．
本論文では福島県伊達市の下小国集会所の線量予測に上記の方法を適用し，実測を決定係数0.9以上で表す式から1年先の予測を行い，決定係数0.7を得た．気象環境の影響を受ける環境中の1年先の放射線量を決定係数0.7前後で予測できたことは一つの成果である．
同観測点の近くで2016年3月末に山火事があった．その影響が決定係数の低下に表れた．ただし2016年末現在，予測線からの分散は大きくなったが，全体としては予測曲線からは離脱していない．ゆえに放射性Csを含む土壌の流出は線量計では検出できない程度と考える．
本研究は平成28年度放射性物質環境動態・環境および生物への影響に関する学際共同研究費(筑波大学)により成されました．アイソトープ環境動態研究センター・環境動態予測部門長の浅沼順教授の支援に感謝いたします．
