LB2045

ベルトールドのLB2045はNaIを使ったガンマ線スペクトルメータである。

本体だけで計測することも,PCを接続して専用ソフトJ-Gammaで計測することもできる。理由は不明だが,両者でcalibration factorが異なる。本体用の値として,ファームウェアバージョン2.20未満では35程度,2.20以上では83程度の値が校正シートに書かれているが,いずれも 137Cs について,前者はBq/cps,後者はBq/L/cps単位のようである(マリネリ容器の容量0.42Lだけ異なる)。J-Gamma用のcalibration factorは固定である。詳細は不明だが,効率を660keVでは0.025,800keVでは0.019,600keVでは0.027としているようである。660keVで効率0.025ということは,137Cs のcalibration factorが 1/0.85/0.025 = 47 に相当する。不可解な仕様だが,ベルトールドジャパンではこれを説明できないらしい。

校正用線源の測定

まず十分ウォーミングアップしてから付属の9kBqの 137Cs 線源でエネルギー校正する。これはすぐ完了し,本体での表示は662keVにぴったり合うように見える。しかし,この状態でRS232C経由でPCにデータを送って,PC側のソフトJ-Gammaで表示すると,かなりずれている。この不具合はベルトールドジャパンでも認識していて,修正される予定とのことである(といいながらちっとも修正されない)。

J-Gammaを使わなくても,PCに取り込まれた結果のテキストファイルをRでフィットすればよい。付属の9kBqの 137Cs 線源を測定した結果が 120615155505.TXT である。このデータ部分を取り出してプロットすると,次のオレンジ色のグラフが得られる(黒部分は後述):

校正用線源のスペクトル

662keVのピークがきれいに見える。この半分のエネルギーを対称の中心として二つ見える崖はCompton edge(コンプトン端(たん))と後方散乱ピークである。正確なピーク位置はおもな放射性核種参照。次の計算で「511」と書いた値は電子質量(keV単位)である(Particle Data Groupの最新の値は 0.510998928 ± 0.000000011 MeV)。

1/(2/511+1/661.657)         # 184.3233 後方散乱
661.657-1/(2/511+1/661.657) # 477.3337 Compton

上のグラフを描くには次のようにする:

sp = scan("http://oku.edu.mie-u.ac.jp/~okumura/stat/data/120615155505.TXT",
          skip=39, nlines=1024, comment.char=";") # スペクトルを読む
sp = floor(sp * 601 + 0.5) # cpsからカウント(整数)に変換。601は秒数
ch = 0:1023  # チャンネル番号
par(mgp=c(2,0.8,0))  # Rのグラフィックオプションを好みでいじる
plot(ch, sp, type="l", xlab="Channel", ylab="Counts", xlim=c(0,400), col="#f39800")

横軸はチャンネル番号(0〜1023)のままにした。400以上は特にフィーチャーがないので省略した。

まずはコベル法でピークカウントを積分してみる。出力はコメント(# で始まる行)で示す。

sum(sp[301:379]) - 79 * mean(sp[c(280:300, 380:400)])
# [1] 244608
sum(sp[311:369]) - 59 * mean(sp[c(280:310, 370:400)])
# [1] 244520.4

1次式+正規分布でフィットしてみる。

# nlmeパッケージをインストールしていない場合:
# options(repos="http://cran.ism.ac.jp") # 例えば
# update.packages()
# install.packages("nlme")
library(nlme)
data = data.frame(ch=ch, sp=sp)
r = 300:380
fit = gnls(sp ~ c + d*(ch-340) + e*dnorm((ch-m)/s),
           data=data, subset=r,
           start=list(c=100,d=-3,e=30000,m=341,s=10),
           weights=varPower(fixed=0.5),
           control=list(nlsTol=1e-4))
points(ch[r], fitted(fit), type="l")
a = coef(fit)
a
#            c            d            e            m            s 
#   129.183659    -3.278402 28038.778015   341.573936     8.731706 
sum(a['e']*dnorm((ch-a['m'])/a['s']))
# [1] 244826.4

どの方法で計算しても,ピーク面積は244500〜244800カウント程度である。最後のピークフィットの値を使うと,ピークの中心チャンネル a['m'] は 341.57 で,cps(counts per second)値は 244826.4 / 601 = 407.365 である。

ここで上の 120615155505.TXT の中を読めばわかるように,本体ではcps値を 504.8 と認識している。ROI(region of interest)すなわち積分範囲が広すぎるのかもしれない。

このcps値に calibration factor を掛けたものが 137Cs のBq値である。ここで,LB2045のファームウェアのバージョンが2.20より古いものは calibration factor として35程度の値を持っているが,この単位は Bq/cps である。より新しいものでは calibration factor が83程度であるが,この単位は Bq/L/cps で,LB2045のマリネリ容器の容量 0.42L で割った値になっている。

本機の場合,calibration factor が 83.0 であり,重量として(わざと変な値)4kg を設定したので,本体の 504.8cps を使って計算すると,

> 504.8 * 83 * 0.42 / 4
[1] 4399.332

となるべきであるが,本体表示値は 4397Bq/kg である。ほぼ一致していると見てよかろう(微妙な違いが出る理由は不明)。一方,上のフィット値を使えば

> 407.365 * 83 * 0.42 / 4
[1] 3550.186

となり,やや小さい。

この時点で,LB2045の 137Cs の662keVでの効率は0.025程度であるという裏情報を得た。137Cs の1Bqは1秒間に0.851個の662keVガンマ線を放出するので,これは 1 / 0.851 / 0.025 = 47 (Bq/cps) に相当する。この値を使えば,

> 407.365 * 47 / 4
[1] 4786.539

となり,今度は少し大きめに出る。

以上をまとめると,Cs-137 の calibration factor は 34Bq/cps ないし 47Bq/cps 程度である。

飯舘村の砂(失敗)

成功した解析は下にあるが,失敗例も載せておく。

データは 120615162546.TXT である。

sp = scan("http://oku.edu.mie-u.ac.jp/~okumura/stat/data/120615162546.TXT",
          skip=39, nlines=1024, comment.char=";")
sp = floor(sp * 600 + 0.5)
data = data.frame(ch=ch, sp=sp)

単純に2次式+三山でフィット。

fit = gnls(sp ~ c + d*(ch-350) + b*(ch-350)^2
           + e1*dnorm((ch-m1)/s1)
           + e2*dnorm((ch-m2)/s2)
           + e3*dnorm((ch-m3)/s3),
           data=data, subset=250:450,
           start=list(c=200,d=-5,b=0.01,e1=10000,e2=10000,e3=10000,
                      m1=300,m2=340,m3=400,s1=10,s2=10,s3=10),
           weights=varPower(fixed=0.5),
           control=list(nlsTol=1e-4)) # 収束しないときは1e-3に直す
a = coef(fit)
a
#             c             d             b            e1            e2 
#  2.112697e+02 -4.976858e+00  3.468806e-02  9.032938e+03  1.008448e+04 
#            e3            m1            m2            m3            s1 
#  6.078507e+03  3.124179e+02  3.430158e+02  4.070545e+02  1.340072e+01 
#            s2            s3 
#  7.821025e+00  1.015221e+01
plot(ch, sp, type="l", xlab="Channel", ylab="Counts", xlim=c(0,500), col="#f39800")
points(data$ch[250:450], fitted(fit), type="l")
points(data$ch[250:450], type="l",
       a['c']+a['d']*(data$ch[250:450]-350)+a['b']*(data$ch[250:450]-350)^2)
sum(a['e1']*dnorm((data$ch-a['m1'])/a['s1']))
# [1] 121047.8
sum(a['e2']*dnorm((data$ch-a['m2'])/a['s2']))
# [1] 78870.94
sum(a['e3']*dnorm((data$ch-a['m3'])/a['s3']))
# [1] 61710.29
飯舘村の砂のスペクトル1

あまり合わない。

飯舘村の砂(成功)

データは 120615162546.TXT である。

sp = scan("http://oku.edu.mie-u.ac.jp/~okumura/stat/data/120615162546.TXT",
          skip=39, nlines=1024, comment.char=";")
sp = floor(sp * 600 + 0.5)
ch = 0:1023
par(mgp=c(2,0.8,0))
plot(ch, sp, type="l", xlab="Channel", ylab="Counts", xlim=c(0,500), col="#f39800")

マイナーなピークも入れる(より正確な値はおもな放射性核種参照)。

Cs-134Cs-137
563keV569keV605keV796keV802keV662keV
8.4%15.4%97.6%85.5%8.7%85.1%

ただし,パラメータが多くなりすぎないように,適宜まとめる。ピーク幅はエネルギーの平方根に比例すると仮定,効率曲線はベキ乗則に従うとした。また,和田先生のご指摘により,ピークの係数を s で割った形にした。

r = 260:450
data = data.frame(ch=ch, sp=sp)
fit = gnls(sp ~ (c + d*(ch-350) + b*(ch-350)^2
                   + (e1*(0.08338/0.9226)*dnorm((ch-0.9314*m1)/(0.9226*s)) +
                      e1*(0.15373/0.9276)*dnorm((ch-0.9415*m1)/(0.9276*s)) +
                      e1*(0.9762/0.956)*dnorm((ch-m1)/(0.956*s)) +
                      e2*0.851*dnorm((ch-m2)/s) +
                      e1*(0.8546/1.0967)*dnorm((ch-m3)/(1.0967*s)) +
                      e1*(0.08688/1.101)*dnorm((ch-1.00765*m3)/(1.101*s))) / s)
                * (ch/350)^u,
           data=data, subset=r,
           start=list(c=200,d=-3,b=0.02, # 初期値
                      e1=10000,e2=10000, # 収束しないときは適宜直す
                      m1=313,m2=342,m3=407,s=9,u=-1.2),
           weights=varPower(fixed=0.5),
           control=list(nlsTol=1e-4)) # 収束しないときは1e-3に直す
a = coef(fit)
plot(ch, sp, type="l", xlab="Channel", ylab="Counts", xlim=c(0,500), col="#f39800")
points(ch[r], fitted(fit), type="l")
points(ch[r], type="l",
       (a['c']+a['d']*(ch[r]-350)+a['b']*(ch[r]-350)^2)
       * (ch[r]/350)^a['u'])
a
#             c             d             b            e1            e2 
#  1.726192e+02 -3.278716e+00  2.336625e-02  7.757909e+04  1.103418e+05 
#            m1            m2            m3             s             u 
#  3.134416e+02  3.415751e+02  4.069815e+02  9.027592e+00 -1.188280e+00 
飯舘村の砂のスペクトル2

図はよく合う。

効率曲線はエネルギーの -1.2 乗くらいということになる。137Cs の量は,カウントを600秒で割ってcpsに直して calibration factor を掛けて

sum(a['e2']*0.851*dnorm((ch[r]-a['m2'])/a['s'])/a['s']*(ch[r]/350)^a['u']) / 600 * 83 * 0.42
# [1] 5621.019
sum(a['e2']*0.851*dnorm((ch[r]-a['m2'])/a['s'])/a['s']*(ch[r]/350)^a['u']) / 600 * (1/0.851/0.025)
# [1] 7579.108

つまり 5621Bq ないし 7579Bq である。134Cs と 137Cs の比は

a['e1']/a['e2']

で 0.7030798 となり,Cs-137校正のベクレルメータの補正でも計算できるように,2012-06-15(大量放出から1.25年)の値として妥当である。

フィッティングの詳細は

summary(fit)

で調べる。関係のある出力を書き出しておく:

Coefficients:
       Value Std.Error  t-value p-value
c     172.62    5.4547   31.646       0
d      -3.28    0.0669  -48.996       0
b       0.02    0.0009   25.938       0
e1  77579.09  359.9130  215.550       0
e2 110341.80  664.3885  166.080       0
m1    313.44    0.0799 3923.023       0
m2    341.58    0.0615 5557.321       0
m3    406.98    0.0672 6056.817       0
s       9.03    0.0361  250.060       0
u      -1.19    0.0336  -35.325       0

 Correlation: 
   c      d      b      e1     e2     m1     m2     m3     s     
d  -0.072                                                        
b  -0.684 -0.512                                                 
e1 -0.484  0.030  0.263                                          
e2 -0.392 -0.089  0.321  0.124                                   
m1  0.083  0.012 -0.019  0.049 -0.323                            
m2 -0.184  0.105  0.086  0.205 -0.156  0.399                     
m3  0.234  0.042 -0.270 -0.083 -0.110  0.031 -0.020              
s  -0.540 -0.085  0.346  0.306  0.321 -0.216 -0.024 -0.099       
u   0.042 -0.748  0.402  0.028  0.206 -0.240 -0.252 -0.088  0.121

各パラメータの相対誤差は t-value の逆数(Std.Error / Value)である。u を固定値と考えれば,134Cs と 137Cs の誤差はほぼ e1e2 の誤差と考えてよい。両者の和の誤差は,

v = vcov(fit)
sqrt(v['e1','e1'] + 2*v['e1','e2'] + v['e2','e2'])

で 793.7626 となる(e1e2 の誤差の和 1024.302 より小さいのは両者に正の相関があるからである)。

別の土

データは 120704151211.TXT である。これは上の例に比べて汚染も少ないので効率曲線はエネルギーの -1.2 乗に固定して計算した。

sp = scan("http://oku.edu.mie-u.ac.jp/~okumura/stat/data/120704151211.TXT",
          skip=38, nlines=1024, comment.char=";")
sp = floor(sp * 601 + 0.5)
data = data.frame(ch=ch, sp=sp)
r = 260:450
fit = gnls(sp ~ (c + d*(ch-350) + b*(ch-350)^2
           + (e1*(0.08338/0.9226)*dnorm((ch-0.9314*m1)/(0.9226*s)) +
              e1*(0.15373/0.9276)*dnorm((ch-0.9415*m1)/(0.9276*s)) +
              e1*(0.9762/0.956)*dnorm((ch-m1)/(0.956*s)) +
              e2*0.851*dnorm((ch-m2)/s) +
              e1*(0.8546/1.0967)*dnorm((ch-m3)/(1.0967*s)) +
              e1*(0.08688/1.101)*dnorm((ch-1.00765*m3)/(1.101*s))) / s)
           * (ch/350)^(-1.2),
           data=data, subset=r,
           start=list(c=200,d=-3,b=0.02,e1=10000,e2=10000,
                      m1=313,m2=342,m3=407,s=9),
           weights=varPower(fixed=0.5),
           control=list(nlsTol=1e-4)) # 収束しないときは1e-3に直す
a = coef(fit)
plot(ch, sp, type="l", xlab="Channel", ylab="Counts", xlim=c(0,500), col="#f39800")
points(ch[r], fitted(fit), type="l")
points(ch[r], type="l",
       (a['c']+a['d']*(ch[r]-350)+a['b']*(ch[r]-350)^2)
       * (ch[r]/350)^(-1.2))
summary(fit)
# Coefficients:
#        Value Std.Error   t-value p-value
# c     8.4041   0.73743   11.3965  0.0000
# d    -0.0626   0.00514  -12.1765  0.0000
# b     0.0004   0.00012    3.1630  0.0018
# e1 1185.0545  35.78182   33.1189  0.0000
# e2 1819.2397  64.36106   28.2662  0.0000
# m1  316.7075   0.46358  683.1753  0.0000
# m2  343.3805   0.33962 1011.0819  0.0000
# m3  409.2899   0.42369  966.0059  0.0000
# s     8.7932   0.22063   39.8554  0.0000
sum(a['e2']*0.851*dnorm((ch[r]-a['m2'])/a['s'])/a['s']*(ch[r]/350)^(-1.2)) / 600 * 83 * 0.42
# [1] 92.11347
sum(a['e2']*0.851*dnorm((ch[r]-a['m2'])/a['s'])/a['s']*(ch[r]/350)^(-1.2)) / 600 * (1/0.851/0.025)
# [1] 124.2013
a['e1']/a['e2']
# 0.651401 
v = vcov(fit)
sqrt(v['e1','e1'] + 2*v['e1','e2'] + v['e2','e2'])
# [1] 80.55856
別の土のスペクトル

上の計算のように 137Cs は 92Bq ないし 124Bq である(サンプルの重量 0.432kg で割れば Bq/kg 単位になる)。相対誤差(1σ)は e2t-value の逆数で3.5%ほど。

みかげさんちのサンプル

みかげさんのSPViewerのLB2045サンプル(白米)を使った例。43201秒(半日)計測。これはよく見ないと山が見えないほどのもの。しかも511にも山がある。

白米の例

上と同じようにできるが,511の山を回避するため,下限を少し上にした。

sp = scan("http://www.mikage.to/radiation/spviewer/lb2045_120122220029.txt",
          skip=35, nlines=1024, comment.char=";")
sp = floor(43201 * sp + 0.5)
data = data.frame(ch=ch, sp=sp)
r = 270:450
fit = gnls(sp ~ (c + d*(ch-350) + b*(ch-350)^2
           + (e1*(0.08338/0.9226)*dnorm((ch-0.9314*m1)/(0.9226*s)) +
              e1*(0.15373/0.9276)*dnorm((ch-0.9415*m1)/(0.9276*s)) +
              e1*(0.9762/0.956)*dnorm((ch-m1)/(0.956*s)) +
              e2*0.851*dnorm((ch-m2)/s) +
              e1*(0.8546/1.0967)*dnorm((ch-m3)/(1.0967*s)) +
              e1*(0.08688/1.101)*dnorm((ch-1.00765*m3)/(1.101*s))) / s)
           * (ch/350)^(-1.2),
           data=data, subset=r,
           start=list(c=200,d=-3,b=0.02,e1=10000,e2=10000,
                      m1=313,m2=342,m3=407,s=9),
           weights=varPower(fixed=0.5),
           control=list(nlsTol=1e-4)) # 収束しないときは1e-3に直す
a = coef(fit)
plot(ch, sp, type="l", xlab="Channel", ylab="Counts", xlim=c(0,500), col="#f39800")
points(ch[r], fitted(fit), type="l")
points(ch[r], type="l",
       (a['c']+a['d']*(ch[r]-350)+a['b']*(ch[r]-350)^2)
       * (ch[r]/350)^(-1.2))
summary(fit)
# Coefficients:
#        Value Std.Error  t-value p-value
# c   167.4171   3.66841  45.6374  0.0000
# d    -0.1892   0.02699  -7.0107  0.0000
# b     0.0015   0.00060   2.4351  0.0159
# e1 1421.5306 137.08970  10.3693  0.0000
# e2 1589.7410 192.99491   8.2372  0.0000
# m1  304.6055   0.97338 312.9351  0.0000
# m2  333.3311   0.97054 343.4482  0.0000
# m3  395.4143   1.08879 363.1687  0.0000
# s     8.9602   0.63998  14.0008  0.0000
sum(a['e2']*0.851*dnorm((ch[r]-a['m2'])/a['s'])/a['s']*(ch[r]/350)^(-1.2)) / 43201 * 83 * 0.42
# [1] 1.158605
sum(a['e2']*0.851*dnorm((ch[r]-a['m2'])/a['s'])/a['s']*(ch[r]/350)^(-1.2)) / 43201 * (1/0.851/0.025)
# [1] 1.562206
a['e1']/a['e2']
# 0.89419
v = vcov(fit)
sqrt(v['e1','e1'] + 2*v['e1','e2'] + v['e2','e2'])
# [1] 303.4506

つまり 137Cs は 1.15Bq ないし 1.56Bq,重量 0.3123kg で割れば 3.7Bq/kg ないし 5.0Bq/kg,その1σ相対誤差は e2t-value の逆数で約12%,つまり 0.45Bq/kg ないし 0.61Bq/kg である。134Cs / 137Cs 比は 0.89419,したがってCs合計は 137Cs の 1.89419 倍で,その1σ相対誤差は 303.4506 / (1421.5306 + 1589.7410) で約10%,つまり 0.71Bq/kg ないし 0.95Bq/kg である。

結論として,半日測れば誤差を 1Bq/kg 以下にできる。誤差の3倍を検出限界とすれば,3Bq/kg くらいは測れる。


奥村晴彦

Last modified: 2012-08-15 21:30:30