CN107290769A - 核电厂点源体源组合的复合辐射源强逆推方法及系统 - Google Patents

Abstract

本发明公开了一种核电厂点源体源组合的复合辐射源强逆推方法及点源体源组合的复合辐射源系统，该方法中，通过探测器获得多个位置的剂量率，将归一化辐射源强在空间上进行离散，利用射线跟踪方法计算出光学距离，结合材料、积累因子等信息开展方程组系数的计算，逆推出源强；然后对探测器位置进行剂量率计算，通过把测量值和计算值进行线性回归分析，计算出标准偏差、斜率、截距等关键参数，进而计算出品质因数来衡量每次计算结果的可接受程度；同时提出加权迭代方法，降低不确定度较大的探测器引入的误差，利用迭代的方式多次重复上述步骤直到品质因数达到预设定值，得到期望的辐射源强信息。

Description

核电厂点源体源组合的复合辐射源强逆推方法及系统

技术领域

本发明涉及核电厂中辐射源强度的计算方法和系统，具体涉及一种核电厂点源体源组合的复合辐射源强逆推方法及点源体源组合的复合辐射源强逆推系统。

背景技术

核电厂的放射性来自压力容器中燃料组件的活性区域，辐射源主要由裂变产物、锕系元素和腐蚀活化产物组成。系统运行中，辐射源随冷却剂流经一回路主系统(包括压力容器、主泵、稳压器、主管道等)、化学容器控制系统等，辐射源分布在冷却剂及相关设备表面。辐射源本身放射性强，工作人员在核电厂正常运行时的日常活动所受的剂量占年总剂量的20％左右，而在核电站大修期间，工作人员所受到的剂量要占到年总剂量的80％，在核电厂大修期间主要是通过缩短工作人员在辐射区的停留时间来减少受照剂量。

核电厂中辐射源分布比较广泛，尤其是经过长时间的使用和大修后，根据工程经验越来越难以推断出各个位置辐射源的辐射强度，所以在很多数据的计算时，尤其是依据辐射源强度的计算中，因为难以获得准确的基础信息而极大地影响了准确性和实用性，同时，在国内目前的防护设施和手段不是很完备的情况下，也使得工作人员被照射的风险大大增加。

现有技术中，需要推算辐射源强度时，一般采用源项分析法，首先，根据放射性物质的产生和消失途径确定其产生项(如流入项，衰变产生项等)和消失项(如过滤项，泄漏项等)，并明确各项的物理模型，然后根据上述各项对放射性物质建立核子浓度平衡方程(组)，最后联立方程(组)求解，然而在这些计算过程中存在大量的简化和近似计算，所以其结果往往与真实数值差距较大，在实际应用时存在很多的障碍，另外，在考虑到辐射源本身对人体造成的危害、核电厂内部复杂的几何结构、放射性核素准确信息难以获取、核电厂探测器测量值的不确定度等因素时，上述方法在安全性、准确性等方面都存在问题，亟待改进或提出新的辐射源强度获取途径。

由于上述原因，本发明人对现有的计算源强信息的方法做了深入研究，根据经验，通常核电厂中部分放射性组件可以简化为点源，还可以简化为简单的体源，如球体源、圆柱体源、长方体源等，并且对该体源进行离散化处理，从而根据离散化处理后的信息进行源强逆推得到辐射源强度信息，从而设计出一种能够解决上述问题的核电厂点源体源组合的复合辐射源强逆推方法及点源体源组合的复合辐射源强逆推系统。

发明内容

为了克服上述问题，本发明人进行了锐意研究，设计出一种核电厂点源体源组合的复合辐射源强逆推方法及点源体源组合的复合辐射源强逆推系统，该方法及系统可以在充分保障人体辐射安全的情况下，得到核电厂内部复杂几何空间结构下的体源源强数据；该方法中，在核电厂中的预定位置放置探测器，并且还在该位置放置带有屏蔽的探测器，进而获得辐射源放出伽马射线的平均能量；另外，在核电厂中还设置有多个监测核电厂辐射值的探测器，以获得部分采点的剂量率，利用点核积分和加权最小二乘法结合的方式，同时将归一化辐射源强在空间上进行离散，利用射线跟踪方法判断每个体源放出的γ射线在空间的穿行路程并计算出光学距离，结合材料、积累因子等信息开展方程组系数的计算，进而逆推出源强；然后得到对探测器位置处剂量率的计算值，把测量值和计算值进行线性回归分析处理，获得标准偏差、斜率、截距等关键参数，进而得到能表示物理含义的品质因数，该品质因数能够衡量每次计算结果的可接受度；同时提出一种加权迭代方法，降低不确定度较大的探测器引入的误差，利用迭代的方式多次重复上述步骤直到品质因数满足预设的条件，进而得到期望的辐射源强和辐射场结果的不确定度，从而完成本发明。

具体来说，本发明的目的在于提供以下方面：

(1)一种核电厂点源体源组合的复合辐射源强逆推方法，其特征在于，该方法包括如下步骤：

步骤一，用探测器探测核电厂内的剂量率D₁,D₂,D₃…D_i，

步骤二，根据探测到的剂量率信息，建立如下式(一)所示的含有辐射源强度的超定方程组，

其中，所述超定方程组的系数矩阵a_i,j通过下式(二)和(三)得到，

在式(三)中，当辐射源为点源时，p＝0，L＝1，M＝1，N＝1；当辐射源为体源时，p＝1；

步骤三，通过最小二乘法计算步骤二中的超定方程组得到辐射源强度信息，所述辐射源强度为下式(四)

S_j,0＝(a_j,i·a_i,j)^-1·a_j,i·D_i (四)

其中，D_i表示第i个探测器探测得到的剂量率；j表示辐射源的个数；m表示辐射源个数能达到的最大值；S_j表示第j个辐射源的强度；S_j,0表示初始计算未进行迭代的第j个辐射源的强度；a_i,j表示系数矩阵，是第j个辐射源对第i个探测器的剂量响应系数；BD(E,L(μ(E),r₀→r_p)表示积累因子，是E和L(μ(E),r₀→r_p)的函数；L(μ(E),r₀→r_p)表示光学距离，是μ(E)和r₀→r_p的函数；μ(E)表示截面/线性衰减系数；r₀→r_p表示辐射源到探测点的距离；C(E)表示通量-剂量转换因子，是E的函数；E表示能量，是核电厂中辐射源发出的伽玛射线的平均能量；表示离散源强，其中的L、M和N分别表示体源在三维坐标上离散后三个坐标轴上的离散标号；

优选地，在步骤三之后，所述方法还包括如下步骤，

步骤四，根据步骤三中得到的辐射源强度信息计算探测器位置处的剂量率，D′₁,D′₂,D′₃…D′₁；

步骤五，对探测器探测到的剂量率信息和计算得到的探测器位置处的剂量率信息进行线性拟合，得到拟合后的两者关系的线性方程，进而得到拟合参数，所述拟合参数包括：平均不确定度、拟合优度和对应的权重矩阵；

步骤六，将步骤五中得到的新的权重矩阵迭代至步骤二中的超定方程组，得到加权的超定方程，进而重复步骤二、步骤三和步骤四，直至获得期望的辐射源强度信息；

其中，D′_i表示计算出的第i个探测器位置处的剂量率。

(2)根据上述(1)所述的核电厂点源体源组合的复合辐射源强逆推方法，其特征在于，所述离散源强通过下式(五)获得：

其中，S_U(L)、S_V(M)和S_W(N)分别表示体源在三维坐标上离散后U坐标轴上的源强权重因子、V坐标轴上的源强权重因子和W坐标轴上的源强权重因子；

优选地，当所述体源为圆柱体源时，S_U(L)、S_V(M)和S_W(N)分别通过下式(六)、(七)和(八)获得：

其中，η_1,1、η_1,2、η_2,1、η_2,2、η_3,1和η_3,2都表示余弦分布常数，R表示圆柱体源的半径，Z表示圆柱体源的高度，表示圆柱体源的角度；

优选地，当所述体源为球体源时，S_U(L)、S_V(M)和S_W(N)分别通过下式(九)、(十)和(十一)获得：

其中，η_1,1、η_1,2、η_2,1、η_2,2、η_3,1和η_3,2都表示余弦分布常数，R表示球体源的半径，θ表示球体源的水平角度，表示球体源的垂直角度；

优选地，当所述体源为长方体源时，S_U(L)、S_V(M)和S_W(N)分别通过下式(十二)、(十三)和(十四)获得：

其中，η_1,1、η_1,2、η_2,1、η_2,2、η_3,1和η_3,2都表示余弦分布常数，x表示长方体源的长度，z表示长方体源的高度，y表示长方体源的宽度。

(3)根据上述(1)所述的核电厂点源体源组合的复合辐射源强逆推方法，其特征在于，对于所述核电厂中辐射源发出的伽玛射线的平均能量E，其测算方法包括如下子步骤：

子步骤1，在核电厂内部选取预定位置，该预定位置距离辐射源的距离为t，在该预定位置放置探测器，收集所述探测器探测到的剂量率I₀，

子步骤2，取回所述探测器，在其外部包覆屏蔽层后放置在所述预定位置，收集所述探测器探测到的剂量率I；

或者，取回所述探测器，在预定位置放置屏蔽体，再将所述探测器放置在屏蔽体内，收集所述探测器探测到的剂量率I；

子步骤3，根据子步骤1和步骤2得到的I和I₀，通过下式(十五)计算包覆层或屏蔽体的质量衰减系数μ，

I/I₀＝BDe^-μt (十五)

子步骤4，根据子步骤3的计算结果，得到辐射源发出的伽玛射线的平均能量E。

(4)根据上述(1)所述的核电厂点源体源组合的复合辐射源强逆推方法，其特征在于，计算所述光学距离L的方法包括如下子步骤，

子步骤a，跟踪伽马射线从辐射源到探测点的穿行过程，记录伽马射线穿过辐射区域的顺序，

子步骤b，分别计算每个辐射区域的距离，结合每个辐射区域材质的线性减弱系数，最后求出总的光学距离L。

(5)根据上述(1)所述的核电厂点源体源组合的复合辐射源强逆推方法，其特征在于，使用最小二乘法处理步骤二中的超定方程组，并获得辐射源强度信息的过程包括如下子步骤：

子步骤3-1，将超定方程组用矩阵的形式表示为AX＝b；

子步骤3-2，求该矩阵的法方程A^TAX＝A^Tb，即X＝(A^TA)^-1A^Tb；

子步骤3-3，用对称矩阵的三角分解法解法方程，记G＝A^TA，其中，G为对称矩阵；

子步骤3-4，利用三角分解法解出G＝LDL^T，其中L是小三角矩阵，D为对角矩阵；

子步骤3-5，解下三角矩阵方程组:LY₁＝A^Tb；

子步骤3-6，解对角矩阵方程组：DY₂＝Y₁；

子步骤3-7，解上三角矩阵方程组：L^TX＝Y₂。

(6)根据上述(2)所述的核电厂点源体源组合的复合辐射源强逆推方法，其特征在于，在步骤五中，通过下式(十六)进行线性拟合，

其中，表示估计的剂量率；表示估计的斜率， 表示估计的截距，

n表示探测器个数i能达到的最大值，表示计算出的探测器位置处剂量率的平均值，表示探测器探测到的剂量率的平均值。

(7)根据上述(6)所述的核电厂点源体源组合的复合辐射源强逆推方法，其特征在于，在步骤五中，根据不确定度得到权重函数，再通过权重函数获得权重矩阵W，所述权重矩阵W通过下式(十七)得到，

其中，f表示拟合不确定度， 表示平均拟合不确定度，f_i表示第i个探测器位置的拟合不确定度；；表示权重函数。

(8)根据上述(6)所述的核电厂点源体源组合的复合辐射源强逆推方法，其特征在于，

在步骤六中，当S_i＞0，且品质因数M达到最大值时停止加权迭代，并输出辐射源强度信息，此时输出的辐射源强度信息即为所述期望的辐射源强度信息；

其中，每次执行步骤六时都相应地得到一个品质因数M，所述品质因数M通过下式(十八)得到，

其中，R²表示拟合优度，

(9)一种核电厂点源体源组合的复合辐射源强逆推系统，其特征在于，该系统用于执行上述(1)-(8)所述的核电厂点源体源组合的复合辐射源强逆推方法。

(10)根据上述(9)所述的核电厂点源体源组合的复合辐射源强逆推系统，其特征在于，该系统包括探测器、伽玛射线平均能量计算模块和辐射源强度计算模块；

所述探测器有多个，包括预定位置探测器和核电厂辐射值监测探测器，

所述预定位置探测器设置在核电厂辐射区域内与辐射源之间距离确定的预定位置，且在所述预定位置探测器外部任选地包覆有可拆卸的屏蔽层；

所述预定位置探测器用于将探测到的辐射剂量率信息传递至伽玛射线平均能量计算模块，

所述核电厂辐射值监测探测器分布在核电厂的辐射区域中，用于将分别探测到的核电厂中剂量率信息传递至辐射源强度计算模块，

所述伽玛射线平均能量计算模块用于计算伽玛射线的平均能量E，

所述辐射源强度计算模块用于计算核电厂中辐射源强度。

本发明所具有的有益效果包括：

(1)根据本发明提供的核电厂点源体源组合的复合辐射源强逆推方法能够在充分保障人体辐射安全的情况下，得到核电厂内部复杂几何空间结构下的体源源强数据；

(2)根据本发明提供的核电厂点源体源组合的复合辐射源强逆推方法通过多次迭代计算，确保最终得到的体源源强信息更为贴紧真实值，具有很高的工程应用价值。

附图说明

图1示出根据本发明一种优选实施方式的整体工作流程图。

具体实施方式

下面通过附图和实施例对本发明进一步详细说明。通过这些说明，本发明的特点和优点将变得更为清楚明确。

在这里专用的词“示例性”意为“用作例子、实施例或说明性”。这里作为“示例性”所说明的任何实施例不必解释为优于或好于其它实施例。

根据本发明提供的核电厂点源体源组合的复合辐射源强逆推方法，该方法包括如下步骤：

步骤一，接收电厂内的探测器探测到的剂量率信息D₁,D₂,D₃…D_i，为了探测上述多个剂量率，需要用到多个探测器，本发明中，可以向电厂中放置多个探测器，也可以直接利用核电厂中已经存在的探测器，核电厂中已经存在的探测器为核电厂辐射值监测探测器，还可以结合使用上述两种方式，对于上述探测器所在的位置要求是：辐射源与所述位置之间没有屏蔽体，本发明中所述剂量率为辐照剂量率。本发明中，所述探测器的数量大于核电厂中辐射源的数量。

步骤二，根据探测到的剂量率信息，建立含有辐射源强度的超定方程组，所述超定方程组为下式(一)，

步骤三，通过最小二乘法计算步骤二中的超定方程组得到辐射源强度信息，所述辐射源强度为下式(四)

S_j,0＝(a_j,i·a_i,j)^-1·a_j,i·D_i (四)

所述超定方程组的系数矩阵a_i,j是在将体源在辐射空间坐标上离散后通过下式(二)和(三)得到的，

在式(三)中，当辐射源为点源时，p＝0，L＝1，M＝1，N＝1；当辐射源为体源时，p＝1；即所述辐射源可以是体源也可以是点源。

本发明中，由于一个核电厂中需要测算辐射源强度的区域很多，需要测算的辐射源也有很多，在不同的区域或者针对不同的辐射源，可以使得选择将其拟化成点源或者体源，并且都可以通过上述式(三)予以测算，当辐射源为点源时，p＝0，L＝1，M＝1，N＝1；当该辐射源是体源时，p＝1。当辐射源是点源时，可以不对辐射源做离散处理，即

经过步骤三之后，已经能够得到辐射源的强度信息，但是该强度信息可能并不够准确，所以通过下述步骤继续计算，以便获得更为贴近真实值的辐射源强度信息；

步骤四，根据步骤三中得到的辐射源强度信息计算探测器位置处的剂量率，D′₁,D′₂,D′₃…D′_i；

步骤五，对探测器探测到的剂量率信息和计算得到的探测器位置处的剂量率信息进行线性拟合，得到拟合后的两者关系的线性方程，进而得到拟合参数，所述拟合参数包括：平均不确定度、拟合优度和对应的权重矩阵；本发明中所述的权重矩阵可以是内权重矩阵或者外权重矩阵，其获得方法是一致的，区别在于外权重矩阵的话不确定度不是由系统计算得到的，而是由操作者手段输入的探测器误差范围。

步骤六，将步骤五中得到的权重矩阵迭代至步骤二中的超定方程组，得到加权的超定方程，进而重复步骤二、步骤三和步骤四，直至获得期望的辐射源强度信息；

本发明中，D表示探测器探测到的剂量率；D_i表示第i个探测器探测得到的剂量率；i表示探测器的个数；j表示辐射源的个数，m表示辐射源个数能达到的最大值；S表示辐射源的强度；S_j表示第j个辐射源的强度；S_j,0表示初始计算未进行迭代的第j个辐射源的强度；a_i,j表示系数矩阵，是第j个辐射源对第i个探测器的剂量响应系数，本发明中即表示点源对探测器的响应系数也表示体源离散后对探测器的响应系数；BD(E,L(μ(E),r₀→r_p)表示积累因子，是E和L(μ(E),r₀→r_p)的函数；L(μ(E),r₀→r_p)表示光学距离，是μ(E)和r₀→r_p的函数，即，光学距离是能量和实际距离的函数；μ(E)表示线性衰减系数；r₀→r_p表示辐射源到探测点的距离；C(E)表示通量-剂量转换因子，是E的函数；E表示能量，是核电厂中辐射源发出的伽玛射线的平均能量；D′_i表示计算出的第i个探测器位置处的剂量率；表示离散源强，其中的L、M和N分别表示体源在三维坐标上离散后三个坐标轴上的离散标号。其中，所述探测点表示探测器的位置，更准确的说是探测器上接收到辐射信息的位置。

本发明中所述的离散源强通过下式(五)获得：

其中，S_U(L)、S_V(M)和S_W(N)分别表示体源在三维坐标上离散后U坐标轴上的源强权重因子、V坐标轴上的源强权重因子和W坐标轴上的源强权重因子；

优选地，当所述体源为圆柱体源时，S_U(L)、S_V(M)和S_W(N)分别通过下式(六)、(七)和(八)获得：

其中，η_1,1、η_1,2、η_2,1、η_2,2、η_3,1和η_3,2都表示余弦分布常数，所述余弦分布常数默认值为零，可以根据实际情况进行设置，Z表示圆柱体源的高度，表示圆柱体源的角度，所述圆柱体源的角度就是圆柱的旋转角度；即，

优选地，当所述体源为球体源时，S_U(L)、S_V(M)和S_W(N)分别通过下式(九)、(十)和(十一)获得：

其中，η_1,1、η_1,2、η_2,1、η_2,2、η_3,1和η_3,2都表示余弦分布常数，所述余弦分布常数默认值为零，可以根据实际情况进行设置，R表示球体源的半径，θ表示球体源的水平角度，表示球体源的垂直角度；即，

优选地，当所述体源为长方体源时，S_U(L)、S_V(M)和S_W(N)分别通过下式(十二)、(十三)和(十四)获得：

其中，η_1,1、η_1,2、η_2,1、η_2,2、η_3,1和η_3,2都表示余弦分布常数，所述余弦分布常数默认值为零，可以根据实际情况进行设置，x表示长方体源的长度，z表示长方体源的高度，y表示长方体源的宽度；即

本发明中所述的积累因子是本领域中常用的专业名词，可参照本领域中的通常含义进行解释和计算，本发明中给出其一般情况下的计算公式如下：

其中K^x的拟合公式如下：

K(E,x)＝cx^a+d[tanh(x/X_k-2)-tanh(-2)]/[1-tanh(-2)]；

其中E为光子能量，MeV；_x为源点到计算点的距离；b为一个平均自由程处的积累因子；a,c,d,X_k为经验参数，累因子系数选择时，可以选择采用对数差值方式，即：

a(E_a)＝{a(E₁)·[log(E₂)-log(E_a)]+a(E₂)·[log(E_a)-log(E₁)]}/[log(E₂)-log(E₁)]

在一个优选的实施方式中，所述核电厂中辐射源发出的伽玛射线的平均能量E的测算方法包括如下子步骤：

子步骤1，在核电厂内部选取预定位置，该预定位置距离辐射源的距离为t，在该预定位置放置探测器，收集所述探测器探测到的剂量率I₀，

子步骤2，取回所述探测器，在其外部包覆屏蔽层后放置在所述预定位置，收集所述探测器探测到的剂量率I；

或者，取回所述探测器，在预定位置放置屏蔽体，再将所述探测器放置在屏蔽体内，收集所述探测器探测到的剂量率I；

子步骤3，根据子步骤1和步骤2得到的I和I₀，通过下式(十五)计算包覆层或屏蔽体的质量衰减系数μ，

I/I₀＝BDe^-μt (十五)

子步骤4，根据子步骤3的计算结果，查表得到辐射源发出的伽玛射线的平均能量E。所述查表的表可以是材料截面表，该表记载在ANSI/ANS 6.4.3，“Gamma-rayAttenuation Coefficients and Buildup Factor for Engineering Materials”,American Nuclear Society，1991.的16-67页。在本发明中，所有的用辐射能量，都用上述平均能量计算，如果核电厂中不同区域其能量差异较大，可以考虑对该区域单独测算，即单独测算平均能量，单独测算辐射源强度。

在一个优选的实施方式中，计算所述光学距离L的方法包括如下子步骤，子步骤a，跟踪伽马射线从辐射源到探测点的穿行过程，记录伽马射线穿过辐射区域的顺序，即通过射线跟踪法算算出辐射源到探测点的距离r₀→r_p，其中r₀表示辐射源的位置，表示探测点的位置r_p。子步骤b，分别计算每个辐射区域的距离，结合每个辐射区域材质的线性减弱系数，最后求出总的光学距离L。

具体来说，计算γ射线穿行路程时，用组合几何方法描述空间，并且将不同介质的空间划分为不同的区域。分别求出伽马射线与每一个基本体的交点与入口之间的距离Di和与出口之间的距离Do。求出每个区域中带“+”和“-”的所有基本体编号plus和minus，该过程可包括如下的六个步骤，

(1)每条线的起点r₀所在区域号Ipstart的确定：

倘若某个区域中没有“-”基本体，那么该区域中，所有的"+"基本体都必须满足起点r₀位于所有的“+”基本体中，那么可以认为射线的起始区域为该区域；倘若该区域中有“-”基本体，那么该区域中所有的"+"基本体都必须满足起点r₀位于所有的“+”基本体中，并且所有的"-"基本体都必须满足不包含此射线的起点r₀，那么认为射线的起始区域为该区域。

(2)每条线的终点r_p所在区域号Ipend的确定：

同样的，倘若某个区域中没有“-”基本体，那么该区域中，所有的"+"基本体都必须满足终点r_p位于所有的“+”基本体中，那么可以认为射线的终止区域为该区域；倘若该区域中有“-”基本体，那么该区域中所有的"+"基本体都必须满足终点r_p位于所有的“+”基本体中，并且所有的"-"基本体都必须满足不包含此射线的终点r_p，那么认为射线的终止区域为该区域。

(3)每条线的起点r₀所在区域号对应的区域出口距离Zo的确定：

若γ射线起始区域号中没有“-”基本体，起始区域中所有的"+"的基本体取出口距离Do中最小者为该γ射线起始区域的出口距离。若γ射线起始区域中有“-”基本体，首先起始区域中所有"+"的基本体取出口距离Do中最小的，然后所有"-"基本体取进口距离Di中最小的，取两者的最大值为该γ射线起始区域的出口距离。

(4)射线经过的每个区域的编号IP的确定：

终点不在最外层的区域情况下，若区域号中没有“-”基本体，首先进行相邻子区域的判断，对于所有的“+”基本体，基本体进口距离小于等于区域的进口距离并且小于基本体的出口距离时(Di<＝Zin<Do)，该区域为上一个区域的相邻区域，求出相应的区域编号IP；若γ射线区域号中有“-”基本体，对于所有的“+”基本体，基本体进口距离小于等于区域的进口距离并且小于基本体的出口距离(Di<＝Zin<Do)，并且对所有“-”基本体，基本体进口距离大于区域进口距离或者基本体出口距离小于等于区域进口距离(Di>Zin或Do<＝Zin)时，该区域为上一个区域的相邻区域，求出相应的区域编号IP。

(5)射线经过的每个区域的进口距离Zi和出口距离Zo的确定：

若上述求出的相邻区域中没有“-”基本体，该区域出口距离即是所有"+"的基本体出口距离Do中最小的，该区域的进口距离即为上一区域的出口距离；若上述求出的相邻区域中有“-”基本体，先求出所有"+"的基本体取出口距离Do中最小者，在求出所有"-"的基本体进口距离Di中最小者，然后取两者的最大值为该区域的出口距离，该区域的进口距离即为上一区域的出口距离。

(6)终点在最外层的区域情况下，首先找到最外层所有的基本体编号aa，对于区域中“+”基本体包含基本体aa,“-”基本体不包含基本体aa的区域时，寻找是否存在射线经过区域中”-”基本体的进口距离大于区域的进口距离(Di(k,minus(i,m))>Zi(k,n))的“-”基本体，如果存在，区域出口距离取所有"-"基本体进口距离Di中最小者，如果不存在，区域出口距离为射线长度。

跟踪到r_p点所在区域IPend并且射线出口距离等于射线长度时终止。从而得到γ射线穿行路程。

再进行射线穿行区域的次数与每次穿行距离的计算：

若射线的穿行区域编号不为0，那么该区域γ射线的穿行距离等于区域进口距离减掉区域出口距离，γ射线穿行次数加1；若γ射线的穿行区域编号为0，停止跟踪。

利用γ质量衰减系数和区域介质的材料获得γ截面μ_n；

通过上述记录γ射线经过区域时的穿行过程，分别求出每一区域的光学距离然后再求和，即：其中，N表示辐射区域的数量，该数量主要是由厂房内部环境决定的。

在一个优选的实施方式中，使用最小二乘法处理步骤二中的超定方程组，并获得辐射源强度信息的过程包括如下子步骤：

子步骤3-1，将超定方程组用矩阵的形式表示为AX＝b；

子步骤3-2，求该矩阵的法方程A^TAX＝A^Tb，即X＝(A^TA)^-1A^Tb；

子步骤3-3，用对称矩阵的三角分解法解法方程，记G＝A^TA，其中，G为对称矩阵；

子步骤3-4，利用三角分解法解出G＝LDL^T，其中L是小三角矩阵，D为对角矩阵；

子步骤3-5，解下三角矩阵方程组:LY₁＝A^Tb；

子步骤3-6，解对角矩阵方程组：DY₂＝Y₁；

子步骤3-7，解上三角矩阵方程组：L^TX＝Y₂。

其中，X＝(A^TA)^-1A^Tb与S_j,0＝(a_i,j ^T·a_i,j)^-1·a_i,j ^T·D_i相对应，通过所述子步骤3-1至子步骤3-7得到辐射源强度的计算数值，本发明中所述的最小二乘法为本领域中通用的超定方程解算方法。

在一个优选的实施方式中，在步骤五中，通过下式(十六)进行线性拟合，

其中，表示估计的剂量率；表示估计的斜率， 表示估计的截距，

n表示探测器个数i能达到的最大值，表示计算出的探测器位置处剂量率的平均值，表示探测器探测到的剂量率的平均值。

在一个优选的实施方式中，在线性拟合之后，还要分别得到线性拟合的平均不确定度、拟合优度、品质因数、加权函数和对应的权重矩阵，所述品质因数代表本次迭代计算的可信度。在步骤五中，根据不确定度得到权重函数，再通过权重函数获得权重矩阵W，所述权重矩阵W通过下式(十七)得到，

其中，f表示拟合不确定度，f_i表示第i个探测器位置的拟合不确定度；表示平均拟合不确定度， 表示权重函数。

在一个优选的实施方式中，在步骤六中，获得所述期望的辐射源强度信息的判断条件是当S_i＞0时，且品质因数M达到最大值，即当S_i＞0，且品质因数M达到最大值时停止加权迭代，并输出辐射源强度信息，该辐射源强度信息就是最终得到的期望的辐射源强度，也是最接近真实值的辐射源强度。

本发明的目的在于获得最接近真实值的辐射源强度，而步骤三得到的辐射源强度的可信度比较低，其与真实值之间的误差会比较大，所以为了提高该数值的准确性，即获得最接近真实值的辐射源强度，本发明中给出了步骤四至步骤六的加权迭代过程，并最终设定了迭代终止的条件，以在保证结果准确的情况下尽量减少工作量，缩短作业时间，提高数据获取的效率。根据上述加权迭代方法和迭代终止的判断标准。另外，本发明中得到的辐射源强度比源项分析的方法获得的辐射源强度更准确，更为贴近真实值，能够保证获得值与真实值在一个数量级之内。在一个优选的实施方式中，每次执行步骤六时都相应地得到一个品质因数M，所述品质因数M通过下式(十八)得到，

其中，R²表示拟合优度，

在一个优选的实施方式中，超定方程组的矩阵形式见下式(十九)

其中，ε表示每个探测器引入的误差；考虑物理涵义，实际上每个探测点处引起的误差可以认为是辐射源引起的，那么上述方程简化为下式(二十)，

进而能够发现，系数矩阵a_i,j等价于第j个辐射源对第i个探测器的剂量响应系数，其中，探测器的剂量响应系数是利用点核积分技术计算的，所述点核积分技术计是本领域中常规的计算方法。

根据本发明提供的一种核电厂点源体源组合的复合辐射源强逆推系统，该系统用于执行本发明上文中所述的核电厂点源体源组合的复合辐射源强逆推方法。

优选地，该系统包括探测器、伽玛射线平均能量计算模块和辐射源强度计算模块；

所述探测器有多个，包括预定位置探测器和核电厂辐射值监测探测器，

所述预定位置探测器设置在核电厂辐射区域内与辐射源之间距离确定的预定位置，且在所述预定位置探测器外部任选地包覆有可拆卸的屏蔽层；所述预定位置距离辐射源的距离可以在后续的计算中作为已知量出现；

所述预定位置探测器用于将探测到的辐射剂量率信息传递至伽玛射线平均能量计算模块，用以计算伽玛射线平均能量；

所述核电厂辐射值监测探测器分布在核电厂的辐射区域中，分别位于本发明中所述的关键位置，用于将分别探测到的核电厂中剂量率信息传递至辐射源强度计算模块，

所述伽玛射线平均能量计算模块用于计算伽玛射线的平均能量E，

所述辐射源强度计算模块用于计算核电厂中辐射源强度。

实验例：

以大亚湾核电站1号机组核岛内NB281房间作为实验对象，该房间为核岛控制区内用来放置装有放射性废水收集桶的地方，废水收集桶为一圆柱形容器，内部放射性液体源强为0.7586E+10MeV/cm³.s(或者4.2898E+14/s)。将该圆柱形容器的上半部分简化为2个点源，下半部分简化为1个圆柱体源，在废水收集桶中间部分每隔50cm设置一个探测器，共五个探测器，每个探测器获得的探测值分别为2.032mSv/hr、0.685mSv/hr、0.255mSv/hr、0.1446mSv/hr、0.0929mSv/hr，即为本发明中的D₁,D₂,D₃,D₄,D₅，根据本发明提供的平均能量获取方法及系统得到平均能量为1.3MeV，采用本发明提供的源强逆推方法及系统，得到点源源强为1.0758E+14/s和1.0757E+14/s，体源源强为2.1363E+14/s。

从最终的结果可知，得到的体源源强和点源源强之和与该辐射源的辐射强度真实值基本一致，所以可以说明本发明提供的方法及系统能够获得贴近真实值的辐射源强度信息。

以上结合了优选的实施方式对本发明进行了说明，不过这些实施方式仅是范例性的，仅起到说明性的作用。在此基础上，可以对本发明进行多种替换和改进，这些均落入本发明的保护范围内。

Claims (10)

1.一种核电厂点源体源组合的复合辐射源强逆推方法，其特征在于，该方法包括如下步骤：

步骤一，用探测器探测核电厂内的剂量率D₁,D₂,D₃…D_i；

步骤二，根据探测到的剂量率信息，建立如下式(一)所示的含有辐射源强度的超定方程组，

其中，所述超定方程组的系数矩阵a_i,j通过下式(二)和(三)得到，

在式(三)中，当辐射源为点源时，p＝0，L＝1，M＝1，N＝1；当辐射源为体源时，p＝1；

步骤三，通过最小二乘法处理步骤二中的超定方程组，得到如下式(四)所示辐射源强度信息，

S_j,0＝(a_j,i·a_i,j)^-1·a_j,i·D_i (四)；

其中，D_i表示第i个探测器探测得到的剂量率；j表示辐射源的个数；m表示辐射源个数能达到的最大值；S_j表示第j个辐射源的强度；S_j,0表示初始计算未进行迭代的第j个辐射源的强度；a_i,j表示系数矩阵，是第j个辐射源对第i个探测器的剂量响应系数；BD(E,L(μ(E),r₀→r_p)表示积累因子，是E和L(μ(E),r₀→r_p)的函数；L(μ(E),r₀→r_p)表示光学距离，是μ(E)和r₀→r_p的函数；μ(E)表示截面/线性衰减系数；r₀→r_p表示辐射源到探测点的距离；C(E)表示通量-剂量转换因子，是E的函数；E表示能量，是核电厂中辐射源发出的伽玛射线的平均能量；表示离散源强，其中的L、M和N分别表示体源在三维坐标上离散后三个坐标轴上的离散标号；

优选地，在步骤三之后，所述方法还包括如下步骤，

步骤四，根据步骤三中得到的辐射源强度信息计算探测器位置处的剂量率，D′₁,D′₂,D′₃…D′_i；

步骤五，对探测器探测到的剂量率信息和计算得到的探测器位置处的剂量率信息进行线性拟合，得到拟合后的两者关系的线性方程，进而得到拟合参数，所述拟合参数包括：平均不确定度、拟合优度和对应的权重矩阵；

步骤六，将步骤五中得到的权重矩阵迭代至步骤二中的超定方程组，得到加权的超定方程，进而重复步骤二、步骤三和步骤四，直至获得期望的辐射源强度信息；

其中，D′_i表示计算出的第i个探测器位置处的剂量率。

2.根据权利要求1所述的核电厂点源体源组合的复合辐射源强逆推方法，其特征在于，所述离散源强通过下式(五)获得：

其中，S_U(L)、S_V(M)和S_W(N)分别表示体源在三维坐标上离散后U坐标轴上的源强权重因子、V坐标轴上的源强权重因子和W坐标轴上的源强权重因子；

优选地，当所述体源为圆柱体源时，S_U(L)、S_V(M)和S_W(N)分别通过下式(六)、(七)和(八)获得：

其中，η_1,1、η_1,2、η_2,1、η_2,2、η_3,1和η_3,2都表示余弦分布常数，R表示圆柱体源的半径，Z表示圆柱体源的高度，表示圆柱体源的角度；

优选地，当所述体源为球体源时，S_U(L)、S_V(M)和S_W(N)分别通过下式(九)、(十)和(十一)获得：

其中，η_1,1、η_1,2、η_2,1、η_2,2、η_3,1和η_3,2都表示余弦分布常数，R表示球体源的半径，θ表示球体源的水平角度，表示球体源的垂直角度；

优选地，当所述体源为长方体源时，S_U(L)、S_V(M)和S_W(N)分别通过下式(十二)、(十三)和(十四)获得：

其中，η_1,1、η_1,2、η_2,1、η_2,2、η_3,1和η_3,2都表示余弦分布常数，x表示长方体源的长度，z表示长方体源的高度，y表示长方体源的宽度。

3.根据权利要求1所述的核电厂点源体源组合的复合辐射源强逆推方法，其特征在于，对于所述核电厂中辐射源发出的伽玛射线的平均能量E，其测算方法包括如下子步骤：

子步骤1，在核电厂内部选取预定位置，该预定位置距离辐射源的距离为t，在该预定位置放置探测器，收集所述探测器探测到的剂量率I₀，

子步骤2，取回所述探测器，在其外部包覆屏蔽层后放置在所述预定位置，收集所述探测器探测到的剂量率I；

或者，取回所述探测器，在预定位置放置屏蔽体，再将所述探测器放置在屏蔽体内，收集所述探测器探测到的剂量率I；

子步骤3，根据子步骤1和步骤2得到的I和I₀，通过下式(十五)计算包覆层或屏蔽体的质量衰减系数μ，

I/I₀＝BDe^-μt (十五)

子步骤4，根据子步骤3的计算结果，得到辐射源发出的伽玛射线的平均能量E。

4.根据权利要求1所述的核电厂点源体源组合的复合辐射源强逆推方法，其特征在于，计算所述光学距离L的方法包括如下子步骤，

子步骤a，跟踪伽马射线从辐射源到探测点的穿行过程，记录伽马射线穿过辐射区域的顺序，

子步骤b，分别计算每个辐射区域的距离，结合每个辐射区域材质的线性减弱系数，最后求出总的光学距离L。

5.根据权利要求1所述的核电厂点源体源组合的复合辐射源强逆推方法，其特征在于，使用最小二乘法处理步骤二中的超定方程组，并获得辐射源强度信息的过程包括如下子步骤：

子步骤3-1，将超定方程组用矩阵的形式表示为AX＝b；

子步骤3-2，求该矩阵的法方程A^TAX＝A^Tb，即X＝(A^TA)^-1A^Tb；

子步骤3-3，用对称矩阵的三角分解法解法方程，记G＝A^TA，其中，G为对称矩阵；

子步骤3-4，利用三角分解法解出G＝LDL^T，其中L是小三角矩阵，D为对角矩阵；

子步骤3-5，解下三角矩阵方程组:LY₁＝A^Tb；

子步骤3-6，解对角矩阵方程组：DY₂＝Y₁；

子步骤3-7，解上三角矩阵方程组：L^TX＝Y₂。

6.根据权利要求1所述的核电厂点源体源组合的复合辐射源强逆推方法，其特征在于，在步骤五中，通过下式(十六)进行线性拟合，

其中，表示估计的剂量率；表示估计的斜率， 表示估计的截距，

n表示探测器个数i能达到的最大值，表示计算出的探测器位置处剂量率的平均值，表示探测器探测到的剂量率的平均值。

7.根据权利要求6所述的核电厂点源体源组合的复合辐射源强逆推方法，其特征在于，在步骤五中，根据不确定度得到权重函数，再通过权重函数获得权重矩阵W，所述权重矩阵W通过下式(十七)得到，

其中，f表示拟合不确定度， 表示平均拟合不确定度，f_i表示第i个探测器位置的拟合不确定度；表示权重函数。

8.根据权利要求6所述的核电厂点源体源组合的复合辐射源强逆推方法，其特征在于，

在步骤六中，当S_i＞0，且品质因数M达到最大值时停止加权迭代，并输出辐射源强度信息，此时输出的辐射源强度信息即为所述期望的辐射源强度信息；

其中，每次执行步骤六时都相应地得到一个品质因数M，所述品质因数M通过下式(十八)得到，

其中，R²表示拟合优度，

9.一种核电厂点源体源组合的复合辐射源强逆推系统，其特征在于，该系统用于执行权利要求1-8所述的核电厂点源体源组合的复合辐射源强逆推方法。

10.根据权利要求9所述的核电厂点源体源组合的复合辐射源强逆推系统，其特征在于，该系统包括探测器、伽玛射线平均能量计算模块和辐射源强度计算模块；

所述探测器有多个，包括预定位置探测器和核电厂辐射值监测探测器，

所述预定位置探测器设置在核电厂辐射区域内与辐射源之间距离确定的预定位置，且在所述预定位置探测器外部任选地包覆有可拆卸的屏蔽层；

所述预定位置探测器用于将探测到的辐射剂量率信息传递至伽玛射线平均能量计算模块，

所述核电厂辐射值监测探测器分布在核电厂的辐射区域中，用于将分别探测到的核电厂中剂量率信息传递至辐射源强度计算模块，

所述伽玛射线平均能量计算模块用于计算伽玛射线的平均能量E，

所述辐射源强度计算模块用于计算核电厂中辐射源强度。