A deconvolution method for scintillator gamma-ray spectrum analysis based on convex optimization

Journal Pre-proof A deconvolution method for scintillator gamma-ray spectrum analysis based on convex optimization Bin Liu, Hongrun Yang, Huanwen Lv, Futing Jing, Xilong Gao, Mingfei Yan

PII: DOI: Reference:

S0168-9002(20)30008-5 https://doi.org/10.1016/j.nima.2020.163399 NIMA 163399

To appear in:

Nuclear Inst. and Methods in Physics Research, A

Received date : 18 November 2019 Revised date : 4 January 2020 Accepted date : 4 January 2020 Please cite this article as: B. Liu, H. Yang, H. Lv et al., A deconvolution method for scintillator gamma-ray spectrum analysis based on convex optimization, Nuclear Inst. and Methods in Physics Research, A (2020), doi: https://doi.org/10.1016/j.nima.2020.163399. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier B.V.

*Manuscript Click here to view linked References

Journal Pre-proof

1

A deconvolution method for scintillator gamma-ray spectrum

2

analysis based on convex optimization

3 4 5 6

Bin Liu a, Hongrun Yang a, Huanwen Lv a, Futing Jing a, Xilong Gao a, Mingfei Yan b

7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42

Abstract: A deconvolution method based on convex optimization is developed and applied to

a

Science and Technology on Reactor System Design Technology Laboratory, Nuclear Power Institute of China, Chengdu, 610213, China RIKEN, Wako, Saitama 351-0198, Japan

of

b

gamma-ray spectrum deonvolution of simulated complicated spectra of NaI(Tl) measurements. The gamma-ray deconvolution problem is transformed into a minimization problem with inequality

pro

constraints. Then the Lagrange dual method is applied to transform the primal Lagrange problem into the dual problem. The interior-point method is applied and the logarithm barrier function is applied to construct the Newton system for iteration search. PCG (Preconditioned Conjugate Gradient) is applied to do the inexact search of the Newton system. The proposed method is applied to deconvolution of 4 typical gamma-ray spectra: a) characteristic gamma-ray spectrum of 137Cs; b) characteristic gamma-ray spectrum of

60

Co; c) characteristic gamma-ray spectrum of

152

Eu; d) mixed fission products in the

re-

primary loop of PWR (Pressurized Water Reactor), including the isotopes of Kr, Xe, I and Cs. Compared with the ML-EM method, improved accuracies and high resolution can be obtained with the proposed method. Moreover, the proposed method costs about 1/3 computational time of the ML-EM method.

Keywords: NaI(Tl) dectector; Gamma-ray spectrum deconvolution; Lagrange dual method; Truncated

urn al P

Newton method

1. Introduction

The sodium iodide NaI(Tl) is widely used in the radiation detection field for its high efficiency and low cost

[1-2]

. However, use of NaI(Tl) is sometimes limited for its well-known low energy resolution,

especially in gamma-ray spectrum measurement. The low energy resolution comes from the following three main reasons: statistical fluctuation in light generation in the crystal for gamma-ray absorption, statistical fluctuation in the number of collected charges in the anode of the photomultiplier tube, and the electronic noise from both the electronics and the detector leakage current [3-4]. As a result of the low energy resolution, there might be overlap between close photo-peaks, which makes it hard to identify radio-nuclides.

The deconvolution techniques with a predetermined detector response are effective to resolve overlapping indications for improving resolution. Many numerical algorithms have been developed for

Jo

gamma-ray spectrum deconvoltuion, such as least square method Maximum Likelihood Expectation Maximization (ML-EM) Gold algorithm

[9]

, Richardson-Lucy method

[10]

[5]

, regularization method

[7]

, Maximum Entropy Method (MEM)

, and metaheuristic algorithm

methods, the ML-EM algorithm is considered to be the most prominent method

[11]

[7, 11]

[6]

,

[8]

,

. Among these

. The mentioned

methods can effectively be used for gamma-ray spectrum deconvolution. However, resolution of the deconvoluted spectrum is still needed to be enhanced, especially for the complicated gamma-ray spectrum. Therefore, a deconvolution method based on convex optimization is proposed in this paper. The proposed method firstly transforms the gamma-ray spectrum deconvolution problem into a convex 1

Journal Pre-proof

53

transform the primal Lagrange problem into the dual problem. The interior-point method is applied and the logarithm barrier function is applied to construct the Newton system for iteration search. PCG (Preconditioned Conjugate Gradient) is applied to do the inexact search of the Newton system. The proposed method is applied to gamma-ray spectrum deconvolution of simulated complicated spectra of NaI(Tl) measurements. Performance of the proposed method is compared with the ML-EM method. 2. Theory 2.1. Gamma-ray spectrum deconvolution A measured gamma ray spectrum

of

52

optimization problem with inequality constraints. And then the Lagrange dual method is applied to

O  E  can be represented by the first Fredholm integral

equation [7],

pro

43 44 45 46 47 48 49 50 51



O  E    R  E , E  s  E   dE 

54

0

(1)

R  E, E is the detector response function, E is the amount of energy deposited in the

where

56

detector, E’ is the incident photon energy, and

57

discretized as,

s  E is the incident γ-ray spectrum. Eq. (1) can be

re-

55

n

Oi   i   Rij s j , i  1,..., m

58

(2)

urn al P

j 1

59

where M is the number of channels in the multi-channel analyzer (MCA), N is the number of

60

mono-energetic γ-lines in the

61 62

detecting an incident γ-ray with an energy of channel i in channel j. Eq. (2) can be expressed in matrix form as,

 O1   1   R1,1 O     R  2    2    2,1           Om   m   Rm,1

63

66 67 68 69 70

R1,2 R2,1 R2,1

R2,1

R1,n   s1  R2,1   s2        R2,1   sn 

(3)

For the convenience of discussion, Eq. (3) is written as,

O  Rs

where

(4)

O  O   .

Jo

64 65

s  E spectrum, ε is the noise term and Rij is the probability of

2.2. Maximum Likelihood Estimation using Expectation Maximization (ML-EM ) The ML-EM algorithm minimizes the Kullback-Leibler information divergence based on Bayes formula [12]. The Kullback-Leibler information divergence is,

 KL  O, Rs    Oi  log Oi  log  Rs i  m

71

i 1

2

(5)

Journal Pre-proof 72

Minimization

 KL

is equivalent to maximizing the likelihood function

J s ,

m

J  s    Oi log  Rs i

73

(6)

i 1

Then the joint probability mass function and the marginal probability mass function can be constructed as,

p X ,Y   j , i; s   Rij s j

76

where X and Y are the random variables constructed.

And according to the Bayes theory, the Q function of the k-th iteration needed for the E–step (Expectation) of ML-EM algorithm is, m

n

Q  s O; sk    rOi log Rij  log s kj  pˆ ijk

80

i 1 j 1

where

Ni  rOi , and r is the positive integer introduced and, pˆ  p X Y  j i ; s k ij

82



k

Rij s kj

R s

k il l

urn al P

The M-step (Maximize likelihood) is the maximization of Q, which yields the final iterative formulation of ML-EM, shown in Eq. (10),

m

s kj 1  s kj 

85

i 1

OiRij

s

k j

Rij

2.3. The deconvolution method based on convex optimization The method is initially developed by S. J. Kim et al

[13]

, which is used to solve the large scale basis

pursuit de-noising (BPDN) problem. Method proposed in this paper is a degradation version of the algorithm developed by S. J. Kim et al. The derivation is illustrated as follows. Initially, problem stated by Eq. (4) can be transformed into an optimization problem with nonnegative constraints of s,

min Rs  O

where

94 95

Jo

92

93

(10)

n

j 1

86 87 88 89 90 91

(9)

n

l 1

83 84

(8)

re-

81

(7)

of

77 78 79

pY  i; s    Rs i

pro

74 75



2 2

denotes the

L2

2 2

(11)

subject to s  0

norm1.

The method developed by S. J. Kim et al is an effective method to solve large scale linear problem Ax=y. The objective function is the combination of L1 norm and L2 norm, which can be written by Eq.

1

As to a vector

x   x1 , x2 ,

, xn  , the L p

norm of x is defined as

equals 1 and 2, the Lp norm becomes L1 norm and L2 norm.

3

x

p

  x12  x22 

xn2 

1 p

. While p

Journal Pre-proof 96

(12),

97

Ax  y 2   x 1

98 99 100

The L1 norm term is used to ensure the sparsity of the reconstructed results. However, for the

102 103

gamma-ray spectra deconvolution problem, there is no need to ensure sparsity of the reconstructed results. So the objective function of gamma-ray spectra deconvolution problem is just the L2 norm

Rs  O 2 . 2

Problem stated by Eq. (11) is a typical convex optimization problem with inequality constraints. Introducing a new variable

z  Rs  O , problem stated by Eq. (11) can be equivalent to, min f  s   zT z

104

subject to z  Rs  O

Eq. (13) can be considered as the primer problem of the gamma-ray spectrum deconvolution. Its Lagrangian is,

L  s, z, v   zT z  vT  Rs  O  z 

107

v

108

where

109

Eq. (15), which can be obtained at the point where

reL

equals zero.

g    inf L  s, z,    1 4 T  T O

urn al P

g    p * for any  

m1

113

Eq. (11). So the minimization problem of Eq. (12) can be equivalent to maximizing

118 119 120 121 122 123 124 125

g   .

The dual feasible point can be constructed as,

  2  Rs  O

(16)

where θ is a constant. The constant θ can be used as the step length of the iteration and it is determined by the backtracking line search method, which will be illustrated in the following part. The difference between the primal objective value of s and the associated lower bound

g   is

called the duality gap, which is described as Eq. (17),

Jo

116 117

L  s, z , v  .

, where p * is the optimal solution of problem stated by

Therefore,

115

(15)

According to the Lagrange dual method [14], the dual function gives the lower bound of

112

114

(14)

is the scalar consist of the associating Lagrange multiplier. The dual function is written as

110 111

(13)

pro

105 106

(12)

of

101

2

  Rs  O 2  g   2

(17)

The duality gap is always nonnegative and becomes zeros at the optimal point. So the duality gap is adopted as the termination condition. The target duality gap

 target used in this paper is 1×10-10.

And then the interior-point method is applied to solve the dual problem. The logarithmic barrier function is defined as,

 s   4

1 n  ln  si  t i 1

(18)

Journal Pre-proof Define a convex function shown in Eq. (19),

  s, t   t  z T z    s  

128

where t varies from 0 to ∞. And for each

129

Then a set called the central path

130

corresponding

131

ˆ

exists and be unique.

can be obtained. For the point on the central path, the

t   , sˆ is the optimal solution of Eq. (11).

is dual feasible. While



in this paper. The Newton system is written in Eq. (20),

H s   g

132 133 134

t  0 , assume sˆ  arg min   s, t 

sˆ t  0

Newton’s method is used to minimize

(19)

of

127

pro

126

(20)

where H is the Hessian and g is the gradient of the current iteration. The Hessian can be written as,

H   2s  s, t 

 t 2s Rs  O 2   2s ln  s  2

re-

135

(21)

 2tRT R  D 0  t  0 D  where,

2 2 D  diag  2 ,..., 2  sn   s1

urn al P

136 137 138

(22)

The gradient is,

g   s  s, t   2 s1   2tR z   ...   2 sn 

139

147 148

Solving the Newton system exactly is computationally costly. So the search direction is always computed as an approximate solution. The inexact search of the Newton system is also called the Truncated Newton method

[15]

. In this paper, the PCG (Preconditioned Conjugate Gradient) method is

applied to solve the Newton system approximately. The PCG method is the CG method with a preconditioned matrix to improve the ill condition of the system

Jo

140 141 142 143 144 145 146

(23)

T

[16].

The preconditioner is a symmetric

and positive definite matrix used to handle and illness of the problem. The preconditioner used in this paper is,

 2tdiag  RT R   D 0  P  0 D  

The update rule of t is,

5

(24)

Journal Pre-proof  max  min 2n  , t ,    min t    min  t ,

149

151

The step length θ is determined by the basic backtracking algorithm backtracking algorithm, the step size is taken as

According to the

while Eq. (26) is satisfied.

  s  s     s   s

152

Finally, the flow chart of the proposed method is shown in Fig. 1.

(26)

re-

pro

153

   , 0    1

[16].

of

150

(25)

154

164 165 166 167 168 169 170

urn al P

Fig. 1 Flow chart of the deconvolution method based on convex optimization 3. Simulation results

3.1 The response function matrix

The response function matrix of NaI(Tl) is obtained by Monte Carlo simulation. Model of the NaI(Tl) detector of the work Davood Alizadeh is adopted in this paper

[11]

. The NaI scintillator is a cylinder

with 7.62 cm diameter and 7.62 cm height, which is coated with 0.185 cm MgO and 0.05 cm Al, shown in Fig. 2. A point gamma-ray source is positioned at the distance of 10 cm from the detector front surface to produce energy distribution of pulses in the NaI(Tl) scintillator.

Jo

155 156 157 158 159 160 161 162 163

Fig. 2. MCNP modeling of the NaI scintillator

Response functions are determined by pulse-height tally of MCNP code. The energy of the gamma-rays is varied from 115 keV to 3080.5 keV in 3 keV steps, making the detector responses binned into 1024 channels. Gaussian broadening of the response function is determined by FT8GEB card of MCNP code, in which the FWHM is defined as:

FWHM  a  b E  cE 2 6

(27)

Journal Pre-proof where a, b and c are coefficients. Coefficients in Davood Alizadeh’s paper are also adopted, which are determined by fitting the experimental data

Four gamma-ray spectra are used to test the deconvolution method: a) characteristic gamma-ray spectrum of 137Cs; b) characteristic gamma-ray spectrum of 60Co; c) characteristic gamma-ray spectrum of

152

Eu; d) a typical gamma-ray spectrum of mixed fission products in the primary loop of PWR

(Pressurized Water Reactor), including the isotopes of Kr, Xe, I and Cs. The characteristic gamma-ray energies and the corresponding probabilities are from the evaluated data base ENSDF.

of

The gamma-ray spectra and their simulated responses are shown in Fig. 3. The Simulated responses are normalized for better visualization. As is shown in Fig. 3, the simulated response of the gamma spectra are contaminated by Compton continuum, escaping peaks and Gaussian broadening of the

pro

photo-peaks.

(a) 137Cs

195 196 197 198

(c) 152Eu

(d) Mixed fission products

Fig. 3. The test gamma-ray spectra and their simulated responses

3.2 The deconvolution results

3.3.1 deconvoluted with ML-EM algorithm

The ML-EM algorithm is applied for deconvolution of the four spectra. The deconvoluted results are shown in Fig. 4. Accuracies of the deconvoluted results are evaluated by RMSE (Relative Mean Square

Jo

194

(b) 60Co

urn al P

184

186 187 188 189 190 191 192 193

. The coefficients a, b and c used in this paper are

-0.0072, 0.0679 and 0.1609 respectively.

183

185

[11]

re-

171 172 173 174 175 176 177 178 179 180 181 182

Error), shown in Eq. (28),

RMSE 

s  sexact sexact

2

(28)

2

where s is the deconvoluted spectrum and sexact is the test spectrum. As is shown in Fig. 4, the deconvoluted spectra of

137

Cs and

60

Co agree well with the exact solution. While the test spectrum

becomes complicated, the deconvoluted spectra deviate from the exact solutions, such as

152

Eu and the

mixed fission products. The RMSE of the deconvoluted results of the test spectra a), b), c) and d) are 7

Journal Pre-proof 199

4.60×10-3, 1.70×10-2, 5.73×10-2 and 2.81×10-1 respectively.

of

200 (a) 137Cs

(b) 60Co

pro

201

210

(d) Mixed fission products

Fig. 4. Deconvoluted spectra with ML-EM method 3.3.2 Deconvoluted with the method based on convex optimization

The deconvoluted spectra with the method based on convex optimization are shown in Fig. 5. Improved accuracies can be obtained with the proposed method. The RMSE of the deconvoluted spectra with the proposed method are 5.49×10-16, 1.64×10-15, 9.15×10-13 and 2.90×10-3 respectively.

(a) 137Cs

(b) 60Co

(c) 152Eu

(d) Mixed fission products

Jo

211

(c) 152Eu

urn al P

203 204 205 206 207 208 209

re-

202

212 213 214

Fig. 5. Deconvoluted with the method based on convex optimization 8

Journal Pre-proof During the deconvolution, ML-EM method performs 1×106 iterations while the proposed method performed 1×104 iterations. The computations were carried on the computer with the Intel i5-4590 CPU with a maximum frequency of 3.3 GHz. The computational costs of the test spectra a), b), c) and d) with ML-EM method are 1.54×103 seconds, 1.55×103 seconds, 1.53×103 seconds and 1.55×103 seconds. While the computational time with the method based on convex optimization are 5.14×102 seconds, 5.17×102 seconds, 5.15×102 seconds and 5.16×102 seconds. The computational time of the proposed method is about 1/3 of the ML-EM method.

of

3.3.4 Resolution of the deconvoluted spectra with the proposed method

Resolution of the deconvoluted spectrum is highly depends on deconvolution accuracy. To evaluate the resolution the proposed method can achieve, the test spectrum with two gamma-ray energy points

pro

of 3.0085 keV and 3.0115 keV is deconvoluted. The energy interval is 3 keV. The deconvoluted spectrum is shown in Fig. 6.

re-

215 216 217 218 219 220 221 222 223 224 225 226 227

228

249 250 251

(a) Whole spectrum

(b) Local enlargement

urn al P

Fig. 6. Resolution of the proposed method

The RMSE of the deconvolution is 3.7811×10-15. As is shown in Fig. 6, due to the high accuracy of the proposed method, the best resolution the MCA (Multiple Channel Analyzer) allows can be achieved.

4. Conclusion and discussion

A deconvolution method based on convex optimization is proposed and applied to gamma-ray spectrum deonvolution of simulated spectra of NaI(Tl) measurements. Gamma-ray deconvolution problem is initially transformed into a minimization problem with inequality constraints. Then the Lagrange dual method is applied to transform the constrained problem into a non-constrained one. And interior-point method is applied and the logarithm barrier function is applied to construct the Newton system for iteration search. PCG algorithm is applied to do the inexact search of the Newton system. Moreover, the deconvolution method is applied to deconvolution of 4 typical gamma-ray spectra: a) characteristic gamma-ray spectra of 137Cs; b) characteristic gamma-ray spectra of 60Co; c) characteristic

Jo

229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248

gamma-ray spectra of 152Eu; d) mixed fission products in the primary loop of PWR (Pressurized Water Reactor), including the isotopes of Kr, Xe, I and Cs. Compared with the ML-EM method, improved accuracies can be obtained with the proposed method. Moreover, computational time of the proposed method is about 1/3 of the ML-EM method.

Reference [1] Experimental methods of nuclear physics. Atomic energy press, 2011 (in Chinese). [2] B. D. Milbrath, B. J. Choate, J. E. Fast, et al. Comparison of LaBr3: Ce and NaI(Tl) scintillators for radio-isotope

9

Journal Pre-proof identification devices. Nucl. Instrum. Methods Phys. Res. A, 572(2007)774-784. [3] T. Burr, M. Hamada. Radio-isotope identification algorithms for NaI γ spectra. Algorithms, 2 (2009) 339-360. [4] A. Likar, T. Vidmar. A peak-search method based on spectrum convolution. J. Physics. D: Appl. Phys. 36 (15) (2003) 1903-1915. [5] L. Bouchet. A comparative study of deconvolution methods for gamma-ray spectra. Astron. Astronphys. Suppl. 113(1995)167-183. [6] M. Morhac, V. Matousec. High-resolution boosted deconvolution of spectroscopic data. J. Comput. Appl. Math, 235(2011)1629-1640.

of

[7] L. J. Meng, D. Ramsden. An inter-comparison of three spectral-deconvolution algorithms for gamma-ray spectroscopy. IEEE Trans. Nucl. Sci., 47(2000) 1329-1336.

[8] J. L. Arcos. Gamma-ray spectra deconvolution using maximum-entropy methods. Nucl. Instr. Methods Phys. Res. A, 369(2)

pro

(1996) 634-636.

[9] P. Bandzuch, M. Morhac, J. Kristiak. Study of the van cittert and gold iterative methods of deconvolution and their application in the deconvolution of experimental spectra of positron annihilation. Nucl. Instrum. Methods Phys. Res. A, 384(2-3) (1997) 506-515.

[10] Zeyu Wang, Pin Gong, Xiaobing Tang, et al. Spectrometry analysis algorithm based on R-L deconvolution and fuzzy interference. Applied Radiation and Isotopes, 153(2019) 1088171-1088176.

Inst. And Methods in Phys. Res. A, 915(2019)1-9.

re-

[11] Davood Alizadeh, Saleh Ashrafi. New hybrid metaheuristic algorithm for scintillator gamma ray spectrum analysis. Nuclear.

[12] Curtis R. Vogel. Computational Methods for Inverse Problems. Society for Industrial and Applied Mathematics, 2002. [13] S. J. Kim, S. Boyd, D. Gorinevsky. An Interior-Point Methods for Large-Scale L1-Regulairzed Least Squares. IEEE Journal of selected topics in signal processing, 1(4), (2007) 606-617.

urn al P

[14] Stephen Boyd, Lieven Vandenberghe. Convex Optimization. Cambridge, U. K.: Cambridge University Press, 2004. [15] Stephen G. Nash. A survey of truncated-Newton methods. Journal of Computational and Applied Mathematics, 124(2000) 45-59.

[16] Jorge Nocedal, Stephen J. Wright. Numerical Optimization. Springer Series in Operations Research, 2006, Springer.

Jo

252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277

10

Journal Pre-proof *Declaration of Interest Statement

Declaration of interests ☐ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Jo

urn al P

re-

pro

of

☒The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

Journal Pre-proof

Jo

urn al P

re-

pro

of

Bin Liu: Methodology, Programing. Hongrun Yang: Supervision. Huanwen Lv: Validation. Futing Jing: Writing and Reviewing. Xilong Gao: Writing draft preparation. Mingfei Yan: Monte Carlo simulation.