三维辐射场是反映一个核设施内真实外照射分布的数据库^[1]。通过分析三维辐射场，不仅可以确定辐射热点的位置，还可以估算出设备或管道中放射性物质的等效活度和分布情况，进而建立有效的屏蔽措施、确定合理的作业方案和降低人员的剂量，实现现场精细化的辐射防护最优化分析^[2]。然而在实际的辐射场重构中，一方面由于源项的活度往往是未知的或存在较大的不确定度，直接将其作为输入用于辐射场模拟计算时会产生较大的误差；另一方面通过测量获得的辐射场数据十分有限，未测量位置处的剂量率通常会用其附近实测的剂量率数值进行插值和外推，而插值和外推往往会引入较大的计算误差。这些都在一定程度上制约了三维辐射场数据以及后续仿真优化技术的实际应用^[3-5]。

事实上，核设施现场的辐射防护人员对所关心区域的热点或源项的位置、分布及数量是很清楚的，同时现场沉积源项所含有的核素种类基本变化不大，可通过探测器测量其特征能量峰及计数获得源项的核素种类与比例，这就使得利用部分辐射场剂量率的测量值来反向计算源项活度成为了可能^[6]。利用计算出的源项活度再进一步计算并重构出未测量位置处的剂量率数据，要比利用简单的插值和外推法得到的数据更加准确可靠。在源活度反演计算与辐射场重构方面，国外的相关研究开展较早，比较有代表性的如比利时核能研究中心(Belgian Nuclear Research Centre, SCK·CEN)和美国电力科学研究院(Electric Power Research Institute, EPRI)。SCK·CEN以点核积分理论为基础，采用了最小二乘法与迭代求解的方式，实现了均匀分布下源活度的反演计算，并开发出了带有源活度反演计算功能的点核积分软件VISIPLAN^[7-8]。EPRI同样采用点核积分理论，通过最小二乘法和线性回归分析实现了均匀分布下源活度的反演计算^[9]。

点核积分方法是辐射场模拟与屏蔽计算常用的方法之一^[10]，其原理简单且易于编程实现，可以处理复杂场所的剂量分布问题，同时具有较高的计算效率，与解析方法和Monte Carlo方法并列为辐射防护剂量模拟计算的三大基本方法。本文基于点核积分理论，采用了最小二乘法和Gauss-Seidel迭代算法，建立了基于源活度反演计算的γ辐射场数据重构方法，并结合某核电现场的实地测量对该重构方法进行了验证与分析，结果表明基于源活度反演计算重构出的辐射场数据与相对应的测量值之间的相对偏差约为10%，能够满足核设施现场辐射防护最优化分析计算的需要。

1 材料和方法

假设某空间内存在n个辐射源项，每个源项的活度均匀分布且第i个源项对应的活度值设为s_i，同时该空间内设置有m个测量位置，设第j个位置处的剂量率测量值为d_j，根据点核积分方法可得第j个位置处的剂量率数据，理论上其与d_j相等，从而构建方程，

$ \begin{array}{*{20}{c}} {\sum\limits_{i = 1}^n {{s_i}} \cdot \int_{{V_i}} {\int_0^{{E_{i\max }}} F } \left( {{E_i}} \right) \cdot }\\ {{K_{i,j}}\left( {{E_i},{r_i} \to {r_j}} \right) \cdot {\rm{d}}E{\rm{d}}r = {d_j}.} \end{array} $

(1)

其中：i的取值为1~n，j的取值为1~m; K_{i, j}(E_i, r_i→r_j)称为点核函数，E_i为γ射线的能量，r_i为第i个源项所在位置，r_j是第j个计算点的位置; F(E_i)是转换因子，用于辐射物理量之间的转换。式(1)实际是一个拥有n个未知数、m个方程的方程组，并且要求m>n，这样便形成了一个超定方程组。

若令h_{j, i}代替式(1)中点核积分部分，如式(2)所示，并同时将式(1)写成矩阵的形式，如式(3)所示。

$ {h_{j,i}} = \int_{{V_i}} {\int_0^{{E_{i\max }}} F } \left( {{E_i}} \right) \cdot {K_{i,j}}\left( {{E_i},\;{r_i} \to {r_j}} \right) \cdot {\rm{d}}E{\rm{d}}r. $

(2)

$ \begin{array}{*{20}{l}} {\mathit{\boldsymbol{H}} \cdot \mathit{\boldsymbol{S}} = \mathit{\boldsymbol{D}}},\\ {\mathit{\boldsymbol{H}} = {{\left( {{h_{j,\;i}}} \right)}_{m \times n}}},\\ {\mathit{\boldsymbol{S}} = {{\left( {{s_1},{s_2}, \cdots ,{s_n}} \right)}^{\rm{T}}}},\\ {\mathit{\boldsymbol{D}} = {{\left( {{d_1},{d_2}, \cdots ,{d_m}} \right)}^{\rm{T}}}.} \end{array} $

(3)

式(3)即为求解源活度的超定方程组，其中将矩阵H称为响应矩阵，其元素h_{j, i}的物理含义是单位活度下第i个源在测量位置j处所造成的剂量率贡献。由于方程数大于未知数的个数，一般超定方程组不存在解向量，只能求出一个最接近的近似解。常用的方法为最小二乘法，即若存在向量S，使得式(4)的值达到最小，则称向量S为式(3)所示的超定方程组的最小二乘解。

$ \sum\limits_{j = 1}^m {{{\left( {{d_j} - \sum\limits_{i = 1}^n {{s_i}} \cdot {h_{j,i}}} \right)}^2}} . $

(4)

根据定理可知，最小二乘解必然存在且唯一，向量S是超定方程组的最小二乘解的充要条件是S是式(5)方程组的解，

$ {\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{H}} \cdot \mathit{\boldsymbol{S}} = {\mathit{\boldsymbol{H}}^{\rm{T}}} \cdot \mathit{\boldsymbol{D}}\mathit{\boldsymbol{.}} $

(5)

式(5)称为式(4)对应的法方程组，其中H^TH是正定矩阵，因此式(5)的解向量存在且唯一。通过求解式(5)方程组的解向量，即可得出满足条件的各个源项的活度值。

基于上述理论，可以将源活度的反演计算转化成求解线性方程组的问题。然而对于不同的计算模型及源项分布，建立起来的方程组的复杂程度也各有不同，同时考虑到方程求解的效率，本文采用Gauss-Seidel迭代算法进行求解，其迭代公式如下：

$ x_i^{(k + 1)} = \frac{1}{{{a_{ii}}}}\left( {{b_i} - \sum\limits_{j = 1}^{i - 1} {{a_{ij}}} x_j^{(k + 1)} - \sum\limits_{j = i + 1}^n {{a_{ij}}} x_j^{(k)}} \right). $

(6)

其中：向量序列{x^(k)}为线性方程组的解向量，{a_ij}是线性方程组的系数矩阵，{b_i}是线性方程组的常数项。若当k趋于无穷时，x^(k)也趋向于稳定值x^*，则x^*就是该线性方程组的解，此时称迭代是收敛的，否则迭代发散。Gauss-Seidel迭代法具有较快的收敛速度，同时若系数{a_ij}为对称正定矩阵，则对于任意迭代初始值x⁽⁰⁾，Gauss-Seidel迭代法均收敛。

根据上述理论，反演计算得到源活度后，将其作为输入并使用点核积分法即可重构出区域的三维辐射场数据。值得指出的是，基于点核积分理论反演的是源的等效活度值，由于点核积分理论的保守性，往往会比真实的活度值偏低，它的目的是用于重构未测量位置处的γ辐射场数据。同时在源活度反演计算过程中，剂量率实测值的位置选取会对源活度反演结果影响较大，需要对其进行深入研究。本文计算所用的点核积分软件由中国辐射防护研究院自主研发，采用了较新的累积因子数据^[11-12]，具有较高的计算精度与效率。

2 结果与讨论 2.1 测量点位选取的影响

在源活度反演计算时，所选取的剂量率测量位置非常重要，尤其对于多个源项的活度计算，需要兼顾每个源项才能使计算结果更加准确。不合适的测量位置的选取会对源活度反演结果产生较大的影响，有时甚至会与真实活度值差几个数量级。本文以算例的形式详细计算并阐述源活度反演计算时所选取的剂量率测量位置对结果的影响。

该算例中存在2个γ射线能量均为0.8 MeV的点源，分别位于(20, 0, 20)和(60, 0, 40)位置处，如图 1所示。屏蔽体为一立方体，材料为水，其在X轴方向上长100 cm、Y轴方向上宽10 cm、Z轴方向高100 cm，其中心坐标为(50, 25, 50)，该算例中坐标的单位均为cm。

利用Monte Carlo软件在测量位置处的剂量率计算结果作为其假定的测量值，在模拟计算时2个点源的活度均设为7.20×10⁹ Bq。同时根据测量点的位置不同，将其分为源1附近的点、源2附近的点、两源中间的点和边缘位置的点，其布局见图 1所示。结合不同测量位置的剂量率，利用点核积分软件对双点源的活度进行反演计算，计算结果如表 1所示。在本算例中通过比较Monte Carlo法与点核积分法在正向计算同一位置处的剂量率数据，发现点核积分法平均要比Monte Carlo法偏大约20%，因此在利用点核积分反向计算2个点源的活度值时，理论上会比真实活度偏小约20%，即约为5.75×10⁹ Bq。

从表 1数据可以看出，源活度反演计算的结果会随着测量点位置的不同产生较大的波动。通过对表 1中的数据进行细致地比较分析，可以得出如下结论。

Case1中所选取的测量点位置均距离源1较近，其剂量率大部分由源1贡献，用该组测量值反演出的源1的活度较为准确可靠，而反演出的源2的活度相对波动较大且不可靠，同时当该组测量值存在误差时，对源1活度的影响较小，对源2活度的影响很大；Case2的情况刚好相反，所选取的剂量率大多由源2贡献，反演计算出的源2的活度值较为准确可靠；Case3中所选取的测量点位，既有源1贡献较大的测量点，又有源2贡献较大的测量点，用该组测量值反演出的源1和源2的活度均比较准确可靠；Case4中测量点处的剂量率数值都比较小，且距离2个源的位置均较远，用该组测量值反演出的源1和源2的活度均不可靠，且当该组测量值即使存在较小的误差也会使得源活度的计算结果波动很大；源活度反演的结果均比理论值7.20×10⁹ Bq要小，是因为相比Monte Carlo法，点核积分法通过引入累积因子来考虑粒子散射的方式要偏保守，所以在反向计算源活度时会使计算结果偏小。

结合表 1给出的不同测量位置处计算的源活度，利用点核积分软件计算同一组未测量点位置处的剂量率，同时利用Monte Carlo软件和预设的源活度计算同一位置处的剂量率，对比结果如表 2所示。可以看出，各组计算结果差异较大，利用Case3给出的活度值计算得到的剂量率结果与同一位置处的剂量率测量值之间差异最小，这也验证了基于该组源活度计算出的剂量率数据的准确性。虽然点核积分软件由于其偏保守性会在反向计算源活度时引入偏小量，但将反演出的源活度作为输入，再次利用点核积分软件正向计算剂量率数据时又会引入偏大量，这样一正一反2次计算会大大降低计算方法本身引入的误差。

2.2 γ辐射场数据重构的现场验证

为了验证本文提出的基于源活度反演计算的辐射场重构方法，在某核电站反应堆厂房中选取了某房间进行现场测量实验, 在该房间内选取了一段管道进行了测量。如图 2所示，现场用铅屏蔽材料盖在辐射热点上，热点主要集中在管道下面的细长圆柱体内。

通过高纯锗探测器分别对热点源项及环境本底进行了测量，得到了热点源项的γ能谱与本底谱。同时用卷尺和周围剂量当量率仪测量了空间的辐射场数据，得到了位置坐标及该位置上的H^*(10)数据，测量时每个位置处取5个读数。在测量过程中发现源项主要集中在U形管的下端，下端表面接触剂量率为5.0 mSv/h，U形管的横管表面接触剂量率为0.11 mSv/h，U形管的竖管表面接触剂量率为0.028 mSv/h。在测量U形管源项时，由于受到来自下端热点的干扰较大，无法用现有屏蔽措施去除，因此在实地测量时，仅测量了下端热点源项的γ能谱，即本实验中热点源项的活度通过测量得到，而U形管源项的活度未知，需要进行反演计算。

通过对热点源项及环境本底的γ能谱进行谱分析处理^[13-14]，得到热点源项的核素种类为⁶⁰Co，该核素主要产生2种能量的γ射线，分别是1.173和1.332 MeV，同时利用Monte Carlo方法计算得到了热点源项的总活度值为2.48×10⁷ Bq。

U形管的管道外径为16.8 cm，管壁厚0.34 cm。将U形管下方热点源项的圆柱分为3段，分别为上段、中段和下段，如图 2所示。它们的管道壁厚均为0.34 cm，上段与下段的外径均为2.4 cm，上段长3.5 cm，下段长3 cm，中段的外径为1.67 cm，长7.5 cm，下段底面距离地面的距离为14 cm。U形管两侧的弯管对应的圆心角均为90°，其左侧竖管长14.2 cm，右侧竖管长88.7 cm。整个U形管道的表面离后墙面的距离为23.5 cm。在利用点核积分软件建模计算时，为了提高计算效率忽略了最上面管道中的阀门的影响，如图 3所示。

该阀门的结构相对比较复杂，不利于三维建模，并且对计算管道正前方的辐射场影响非常小。同时在设置U形管源项时，将其用1根横管和2根竖管代替，简化了计算模型，源项设置均为均匀分布在管道内壁的面源，如图 3中虚线所示。

在图 3坐标系下，现场测量的辐射场原始H^*(10)数据如表 3所示，其中坐标单位为cm，H^*(10)数据单位为μSv/h，同时将测量数据分成2组，一组用于反演U形管的源项活度，另一组作为辐射场计算的验证数据。

由于U形管源项未进行测量，若忽略其影响，在利用点核积分软件正向计算辐射场时，距离热点较近的点剂量率的计算值与测量值相差不大，但随着计算点与热点之间距离的增大，U形管源项的剂量率贡献比重会变大，若只考虑热点源项会使得计算值偏小较多，如表 4所示，表中给出了仅考虑热点源项条件下的第1组测量位置处的剂量率测量值与计算值之间的对比。

为了得到更加准确的剂量率计算结果，需要根据部分剂量率的实测值，反演计算出U形管源项的活度值，然后将其与热点源项同时作为输入，利用点核积分软件再进行正向计算。在进行U形管源活度的反演计算时，由于管道是相通的，假设其源项核素仍为⁶⁰Co，同时在设置源项时进行了简化，用1根横管和2根竖管代替U形管。在反演计算活度时，至少需要4个测量位置处的坐标及H^*(10)数据。为此，使用剂量率测量值中的第1组数据进行U形管源活度的反演计算，计算结果如表 5所示。

结合计算出的U形管源项及测量得到的热点源项，使用剂量率测量值的第2组数据进行对比验证。即利用点核积分软件对第2组测量位置处的H^*(10)进行计算，并将计算结果与对应的测量数据进行对比，结果如表 6所示。从表 6数据可以看出，基于源活度反演计算重构出的辐射场数据，无论是参与源活度反演计算的第1组数据还是未参与计算的第2组数据，均与实际的测量数据符合较好，最大的相对偏差在10%左右，能够较好地解决由于源活度缺失或剂量率测量数据少等造成的辐射场难以准确重构的问题。

导致表 6中剂量率的计算值与测量值之间的差异的主要原因有以下5点：

1) 热点源项活度引入的偏差。热点源项活度是通过高纯锗探测器测量能谱并使用谱分析技术计算得到的，这部分会引入偏差。2)剂量率测量值引入的偏差。剂量率的位置使用卷尺测量，H^*(10)数据则通过周围剂量当量率仪测量，这两部分也会引入测量的偏差。3) U形管源项简化引入的偏差。4)点核积分软件计算引入的误差，即体源的离散、对散射的估计以及迭代计算而引入的误差。5)计算时忽略了墙面与地面的反散射影响。

3 结论

在实际计算中，辐射源项的不确定性会对γ辐射场数据的重构产生非常大的影响，同时通过测量获取的辐射场数据也往往较为有限。基于点核积分理论，本文使用最小二乘法与Gauss-Seidel迭代算法建立了基于源活度反演计算的γ辐射场数据重构方法，并研究了所选取不同测量点位置对源活度反演结果的影响，给出了相应测量点位选取的建议，同时结合某核电现场的实地测量对该重构方法进行了计算验证与分析讨论。结果表明，基于该方法重构出的辐射场数据与实测数据之间的相对偏差约为10%，可以较好地解决由于源活度存在较大不确定度或剂量率测量数据少等造成的γ辐射场难以准确重构的问题，能够满足核设施现场辐射防护计算的需求。

本文的研究结果对于核设施内的三维辐射场可视化计算与虚拟仿真具有较大的意义，能够提升三维γ辐射场模拟计算数据的准确性，使其更加切合真实照射的情况，从而更好地用于可视化显示及作业剂量优化的相关工作，为核设施现场的辐射防护最优化提供先进的技术支持手段。