基于蒙特卡罗减方差方法的图像传感器辐射屏蔽结构多目标优化设计

Radiation Shielding Structure Multi-objective Optimization of Imagine Sensors Based on Monte-Carlo Variance Reduction Method

摘要

针对核辐射环境下图像传感器辐射屏蔽问题，提出基于蒙特卡罗（MC）减方差的辐射仿真方法和屏蔽结构多目标优化设计方法。根据中子和光子耦合输运原理构建图像传感器的多辐射源复合屏蔽模型，并针对其厚屏蔽微探测结构特征提出特征插值权窗MC减方差的屏蔽仿真方法，有效提高了其屏蔽仿真的计算效率与精度；通过将遗传算法与辐射屏蔽结构参数化几何建模结合，构建图像传感器的多目标优化设计方法，以获得体积、质量和剂量率参数下最优解并得到Pareto前沿的非劣解组合。数值实验验证了该图像传感器辐射屏蔽结构优化设计方法的有效性。

Abstract

The radiation simulation method and shielding structure multi-objective optimization design method were proposed based on MC variance reduction for the radiation shielding problems of image sensors in nuclear radiation environments. According to the principles of coupled neutron and photon transport， a multi radiation source composite shielding model was constructed for image sensors， and a shielding simulation method of characteristic interpolation weight window MC variance reduction was proposed for the thick-shielding micro-detection structure characteristics， effectively improving the computational efficiency and accuracy of the shielding simulation. By combining the genetic algorithms with parameterized geometric modeling of radiation shielding structures， a multi-objective optimization design method for image sensors was constructed to obtain the optimal solution under volume， mass， and dose rate parameters and obtain a non-dominated solution combination of Pareto front. Numerical experiments verified the effectiveness of the imagine sensor radiation shielding structure optimization design method.

Graphical abstract

关键词

图像传感器 / 辐射屏蔽 / 多目标优化设计 / 蒙特卡罗粒子输运 / 减方差方法

Key words

imagine sensor / radiation shielding / multi-objective optimization / Monte-Carlo（MC） particle transport / variance reduction method

0 引言

视频监控是辐射环境下运行的核设施及核退役作业场所中获取环境信息最常见、最直观的手段，且通常需应对恶劣的工作环境，其中的电子设备尤其是图像传感器易受到高能中子和γ射线辐射的损伤，导致图像传输效率低，使用寿命缩短，影响整体装备的使用寿命和可靠性^［1］。随着核工业应用场景的扩展，对图像传感器的辐射屏蔽提出更高的要求，不仅要保证屏蔽性能（剂量率低），还要追求结构紧凑（体积小）与轻量化（质量小）。图像传感器在复杂辐射环境下的屏蔽性能分析与仿真模拟以及屏蔽结构的多目标优化，一直是各国研究者的研究热点。核工业中互补金属氧化物半导体图像传感器（complementary metal-oxide-semiconductor image sensor， CIS）因其低功耗、抗辐射性强等优点得到广泛应用，成为智能核设备视频监控前端系统的核心部件^［2］。

近年来，国内外的许多研究者已针对CIS在不同辐射条件下的性能展开实验或蒙特卡罗（MC）仿真研究，验证经不同类型的辐射源（如60Co、质子、中子和电子等）照射后CIS的阈值电压漂移、暗电流和光谱响应的变化情况，分析辐射对CIS产生的负面影响与辐射的总剂量、剂量率、暴露时间和辐射类型等因素的关系^［3-5］。在屏蔽方面，NÁFRÁDI等^［6］研究聚变装置中工作的EDICAM相机CIS的暗图像，并使用MC方法对使用一年的伽马能谱和剂量进行估算，以评估相机机构的辐射屏蔽效果。XUE等^［7］使用GEANT4模拟CIS因中子辐射而导致的暗信号退化，并考虑像素规格、核相互作用和损伤效应与实验交叉验证。徐守龙等^［8］对CMOS APS数字模块的辐射容限进行研究，结合MC设计制造屏蔽加固结构并进行实验对比，以研究其屏蔽性能、工作寿命和辐射损伤模式。在核工业中辐射屏蔽的优化是一个多参数、多目标、非线性问题，为追求低剂量率、小体积和轻质量，SONG等^［9］将智能多目标优化方法应用于核反应堆辐射屏蔽。CHEN等^［10］将遗传算法（GA）与MC方法结合而提出一种辐射屏蔽设计的多目标优化策略。上述研究大多只考虑单一辐射类型，而核工业环境下中子和光子并存且相互影响；此外，基于MC仿真的优化方案因在厚屏蔽-微探测问题上计算效率低下，从而导致优化设计成本高。

针对上述不足，本研究面向核工业环境下中子-光子辐射并存的实际需求，首先设计适配的CIS中-光子复合屏蔽结构，并构建对应的辐射屏蔽MC仿真方案；在此基础上，采用GA对屏蔽结构的多目标参数进行优化设计。此外，针对该结构因厚屏蔽-微探测特性而导致的MC仿真效率低下问题，本文提出特征插值权窗方法，以提高优化迭代过程中的MC计算效率。

1 中-光子辐射机制与屏蔽模型

核工业中辐射源主要为中子与伽马光子。中子与介质主要发生弹性散射、非弹性散射、俘获吸收等反应^［11］，其中弹性散射导致的位移损伤为CIS中子辐射损伤的主要来源^［5］，而发生非弹性散射和俘获吸收会产生光子，光子在介质中会发生康普顿散射、光电效应、对产生、三产生、吸收等作用^［12］，光子作用进而会产生电离损伤。CIS的辐射屏蔽需在光学模块外设计屏蔽壳以抵御中子和伽马射线，因材料对不同能量粒子的屏蔽效果受反应截面影响，且屏蔽中子时需考虑非弹性散射产生的二次伽马射线，故本研究采用直筒复合屏蔽结构，形成一种直筒式设计，如图1所示，前部为双层透光层，外部是双层屏蔽包壳和支撑结构。该设计结构简单、成本可控，通过不同材料组合使该结构具备了对光子和中子的良好综合屏蔽性能^［13］。

MC程序运行依赖构造实体几何（CSG）模型^［14］，需进行相应简化，如图2所示，保留前端双层透光玻璃层，外围双层屏蔽包壳，光学模块中将CIS的芯片基底、基底硅、掺杂层和氧化层等各自简化为单一介质层。

多目标优化的目的是平衡减少辐射剂量与小型化、轻量化目标，为此需明确目标函数、决策变量及约束条件，建立多目标优化数学模型以寻求最佳方案，该模型可表示为

\left.\begin{array}{l}\mathrm{m}\mathrm{i}\mathrm{n}F(\mathit{X})=(F_R(\mathit{X}),F_V(\mathit{X}),F_M{(\mathit{X}))}^{\mathrm{T}}\\ \mathrm{s}.\mathrm{t}.F_R(\mathit{X})=R_{\mathrm{N}}(\mathit{X})+R_{\mathrm{P}}(\mathit{X})\le R_{\mathrm{0}}\\ F_V(\mathit{X})={\displaystyle \sum _{h=\mathrm{1}}^H}{V}_h\le V_{\mathrm{0}}\\ F_M(\mathit{X})={\displaystyle \sum _{h=\mathrm{1}}^H}{V}_h{\rho }_h\le M_{\mathrm{0}}\end{array}\right\}（1）

\mathit{X}=(x_{\mathrm{1}},x_{\mathrm{2}},\dots ,x_{n-\mathrm{1}},x_n)\mathit{X}\in \mathbf{R}

L_i\le x_i\le U_ii=\mathrm{1,2},\dots ,n

h_j(x_i)\le \mathrm{0}j=\mathrm{1,2},\dots ,p

g_k(i)=\mathrm{0}k=\mathrm{1,2},\dots ,q

其中，\mathit{X}为属于可行参数空间\mathbf{R}的n维设计变量向量；x_i为第i个设计变量，受到下限L_i和上限U_i的约束，且存在p个不等式约束h_j和q个g_k等式约束；辐射剂量率F_R(\mathit{X})为中子剂量率R_{\mathrm{N}}(\mathit{X})和光子剂量率R_{\mathrm{P}}(\mathit{X})之和，不大于剂量率约束参数R_{\mathrm{0}}；屏蔽结构总体积F_V(\mathit{X})为H个组成屏蔽层体积V_h之和，不大于体积约束V_{\mathrm{0}}；总质量F_M(\mathit{X})为各屏蔽层体积V_h与对应密度{\rho }_h的乘积之和，不大于质量约束M_{\mathrm{0}}。

2 蒙特卡罗粒子输运过程与源的构建

2.1 蒙特卡罗粒子输运过程

在辐射屏蔽领域，MC方法已成为应对复杂挑战的有力工具，它植根于概率论和统计原理，可模拟辐射粒子在介质中输运的实际物理过程。典型的MC软件具有高效的数据处理和可视化能力，如MCNP、Geant4和OpenMC。其中OpenMC是一个开源的MC代码，可以模拟包括中子、光子、电子与耦合输运等各类问题^［15］。本研究采用OpenMC作为辐射剂量的计算工具。MC方法是对粒子在介质中随机游走过程的直接模拟，根据MC输运原理^［11］，粒子的状态由状态参数描述为

\mathit{S}=(\mathit{r},E,\mathit{\Omega },t,W)（2）

其中，\mathit{r}为粒子的坐标，E为粒子能量，\mathit{\Omega }为粒子的方向，t为时间，W为粒子权重。粒子从第m-1个状态序列{\mathit{S}}_{m-\mathrm{1}}到第m个状态序列{\mathit{S}}_m需确定碰撞位置的坐标{\mathit{r}}_m，且满足

\begin{array}{c}T({\mathit{r}}_{m-\mathrm{1}}\to {\mathit{r}}_m|E_{m-\mathrm{1}},{\mathit{\Omega }}_{m-\mathrm{1}})=\\ {\mathrm{\Sigma }}_t({\mathit{r}}_m,E_{m-\mathrm{1}})\mathrm{e}\mathrm{x}\mathrm{p}(-{\int }_{\mathrm{0}}^l{\mathrm{\Sigma }}_t({\mathit{r}}_{m-\mathrm{1}}+l{\mathit{\Omega }}_{m-\mathrm{1}},E_{m-\mathrm{1}})\mathrm{d}l)·\end{array}

\frac{\delta ({\mathit{\Omega }}_{m-\mathrm{1}}-\frac{{\mathit{r}}_m-{\mathit{r}}_{m-\mathrm{1}}}{\left|{\mathit{r}}_m-{\mathit{r}}_{m-\mathrm{1}}\right|})}{{\left|{\mathit{r}}_m-{\mathit{r}}_{m-\mathrm{1}}\right|}^{\mathrm{2}}}（3）

其中，T({\mathit{r}}_{m-\mathrm{1}}\to {\mathit{r}}_m\left|E_{m-\mathrm{1}},{\mathit{\Omega }}_{m-\mathrm{1}})\right.表示能量和方向分别为E_{m-\mathrm{1}}和{\mathit{\Omega }}_{m-\mathrm{1}}时，由位置{\mathit{r}}_{m-\mathrm{1}}到{\mathit{r}}_m的条件概率密度；迁徙长度l=\left|{\mathit{r}}_m-{\mathit{r}}_{m-\mathrm{1}}\right|；\delta (\cdot )表示狄拉克函数；{\mathrm{\Sigma }}_t(\mathit{r},E)表示粒子能量为E时，\mathit{r}位置的宏观总截面。{\mathit{r}}_m由抽样得到：

{\mathit{r}}_m={\mathit{r}}_{m-\mathrm{1}}-\frac{\mathrm{l}\mathrm{n}\xi }{{\mathrm{\Sigma }}_t({\mathit{r}}_{m-\mathrm{1}},E_{m-\mathrm{1}})}{\mathit{\Omega }}_{m-\mathrm{1}}（4）

其中，\xi为（0，1）之间的随机数。随后通过俄罗斯轮盘赌确定碰撞介质中核的种类以及反应类型。最后，基于前述粒子与介质反应的过程，角度和能量的值通过其概率密度分布函数进行连续性随机变量抽样获得。利用上述过程对粒子的状态序列完成更新，直至粒子被吸收。过程中部分粒子能量沉积在碰撞点（即为吸收剂量率），可定义为

H_{\mathrm{t}}(E)=\varphi (E)\sum _i{\rho }_i\sum _r{k}_{i,r}(E)（5）

其中，H_{\mathrm{t}}(E)为能量沉积，\varphi (E)为通量率，{\rho }_{\mathrm{N},i}为第i种核素的密度，k_{i,r}为第i种核素发生第r种反应的动能释放系数。能量沉积的单位为\mathrm{e}\mathrm{V}/(\mathrm{k}\mathrm{g}·\mathrm{s})，评估辐射强度时常换算成吸收剂量率R（E），单位Gy/h，定义为

R(E)=\mathrm{1.6}\times {\mathrm{10}}^{-\mathrm{19}}\times \mathrm{3600}{H}_{\mathrm{t}}(E)（6）

2.2 辐射源的构建与仿真方案

实验中通常采用Co-60点源或圆柱源作为辐射源，通过调整试样的摆放位置与朝向来测试屏蔽效果的空间均匀性^［16］。仿真实验中可采用薄球壳源抽样构建各向同性的空间辐射场，常被用于航空器辐射MC仿真屏蔽计算^［17］。本研究也采用该方法，并对其适用性展开进一步讨论。

当粒子在薄球壳内均匀分布，球坐标系下的位置为\mathit{r}=(R,\theta ,\phi )，R的概率密度分布函数为^［11］

f(R)=\left\{\begin{array}{ll}\frac{\mathrm{3}{R}^{\mathrm{2}}}{R_{\mathrm{2}}^{\mathrm{2}}-R_{\mathrm{1}}^{\mathrm{2}}}& R_{\mathrm{1}}\le R\le R_{\mathrm{2}}\\ \mathrm{0}& \mathrm{其他}\end{array}\right.（7）

其中，R_{\mathrm{1}}\mathrm{、}{R}_{\mathrm{2}}分别为球壳源内外径。采样点的球坐标中r_{\mathrm{0}}为球半径，{\theta }_{\mathrm{0}}为方位角，{\phi }_{\mathrm{0}}为俯仰角。图3展示了r_{\mathrm{0}}的抽样过程^［11］，其中X为抽样结果。

采样点在笛卡儿坐标系下的坐标(x_{\mathrm{0}},y_{\mathrm{0}},z_{\mathrm{0}})可表示为

\left.\begin{array}{l}{x}_{\mathrm{0}}=r_{\mathrm{0}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{0}}\\ y_{\mathrm{0}}=r_{\mathrm{0}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{0}}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{0}}\\ z_{\mathrm{0}}=r_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{0}}\end{array}\right\}（8）

\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{0}}=\mathrm{2}{\xi }_{\mathrm{4}}-\mathrm{1}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{0}}=\frac{{\xi }_{\mathrm{5}}^{\mathrm{2}}-{\xi }_{\mathrm{6}}^{\mathrm{2}}}{{\xi }_{\mathrm{5}}^{\mathrm{2}}+{\xi }_{\mathrm{6}}^{\mathrm{2}}}

\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{0}}=\frac{\mathrm{2}{\xi }_{\mathrm{5}}{\xi }_{\mathrm{6}}}{{\xi }_{\mathrm{5}}^{\mathrm{2}}+{\xi }_{\mathrm{6}}^{\mathrm{2}}}{\xi }_{\mathrm{5}}+{\xi }_{\mathrm{6}}^{\mathrm{2}}\le \mathrm{1}

其中，{\xi }_{\mathrm{1}}~{\xi }_{\mathrm{6}}为在区间（0，1）内均匀分布的独立随机数。源点的各向同性使得球壳源具有各向同性，但其内部的通量率分布不均，中心点通量率为

\varphi =\frac{I}{\mathrm{4}\mathrm{\pi }{R}_{\mathrm{1}}^{\mathrm{2}}}（9）

其单位为\mathrm{c}{\mathrm{m}}^{-\mathrm{2}}\cdot {\mathrm{s}}^{-\mathrm{1}}；I为源发射强度，单位为{\mathrm{s}}^{-\mathrm{1}}。球壳源通量率分布与设定值的比值R_D如图4所示，通量率值会从中心向四周增大，在内径一半的空间内通量率相较于设定值的误差在10%以下，该误差可以接受，因此仿真设定中球壳源内径大小应当比屏蔽体外接半径大一倍为佳。

3 特征插值减方差方法

3.1 方差和权窗

MC方法以大量粒子模拟的样本均值为计算结果，组内方差可反映近似值与真实值的误差，因此统计区域的粒子密度会直接影响结果准确性，且屏蔽效果越好，统计方差越大，结果越难以收敛。近似值与真实值之间的误差可由近似方程表示为^［11］

（10）

式（10）表示在\mathrm{1}-\alpha的置信水平下，近似值\widehat{Z}与真实值E(Z)的误差为x\sigma /\sqrt[]{N}。增大粒子数N可以提高精度，已知收敛速度为O(N^{-\mathrm{1}/\mathrm{2}})，若精度提高一个数量级，则N需要增大100倍，总时长相应延长100倍^［11］。到达一定精度后增大N的作用将有限，难以解决精度问题。

权窗法通过控制粒子权重来增加重要空间区域和能量范围内的粒子数量以此减少方差，是MC方法中解决效率问题的常用方法^［18］。该方法原理是：为充分利用粒子的历史，即当粒子抽样发生吸收作用时不直接结束其进程，而是减小粒子权重后继续模拟，仅在权重降至一定值时结束。吸收反应后的权重可表示为

W\mathrm{\text{'}}=W(\mathrm{1}-\frac{{\sigma }_{\mathrm{a}}(E)}{{\sigma }_{\mathrm{t}}(E)})（11）

其中，{\sigma }_{\mathrm{a}}为吸收微观截面，表示单个靶核与入射粒子发生相互作用的概率；{\sigma }_{\mathrm{t}}为微观总截面，表示所有可能反应的总概率。即反应后权重W\mathrm{\text{'}}等于反应前权重W与存活概率的乘积。

通过设置空间位置和能量区间的权窗上下限，可选择对粒子进行扩增或提前终止进程，从而分别提高计算精度与缩短计算时间。当粒子权重W大于权窗上限W_{\mathrm{u}}时，则会分裂成M个粒子，可定义为

M=\mathrm{m}\mathrm{i}\mathrm{n}(M_{\mathrm{m}\mathrm{a}\mathrm{x}},[\frac{W}{W_{\mathrm{u}}}])（12）

其中，M_{\mathrm{m}\mathrm{a}\mathrm{x}}为设定的最大分裂粒子数，分裂后粒子权重为W/M。在重要区域设置较低权窗上限，可以增大粒子密度以提高计算精度。当粒子权重W小于权窗下限W_{\mathrm{l}}时，可通过轮盘赌决定是否结束进程，生存权重可定义为

W_{\mathrm{s}}=\mathrm{m}\mathrm{i}\mathrm{n}(M_{\mathrm{m}\mathrm{a}\mathrm{x}}W,f_{\mathrm{s}}{W}_{\mathrm{l}})（13）

其中，f_{\mathrm{s}}为大于1的主观生存率。当粒子权重W小于权窗下限W_{\mathrm{l}}时，抽取区间（0，1）内均匀分布的随机数\xi，若\xi<\mathrm{1}-W/W_{\mathrm{s}}，则结束进程；若\xi≥\mathrm{1}-W/W_{\mathrm{s}}则粒子存活且将权重值设为W_{\mathrm{s}}。在非重要区域设置较大的权窗下限可提前结束低权重粒子计算。

3.2 特征插值权重减方差方法

根据权窗原理，粒子会沿着权窗分布的逆梯度方向运动，合理设置权窗的上限和下限并引导粒子向重要统计区域移动能够节约计算资源，从而达到减小方差、加速计算的目的。合理的权窗分布并不影响结果的准确性，而通过粒子调配影响全局方差分布，从而将有限的计算资源集中在更重要的空间位置，则权窗上下限的分布会直接影响权窗加速效果。

本研究中屏蔽结构仿真的权窗设置具有如下特点：①选稍大球壳源以保证辐射通量率均匀，但会导致许多粒子未抵达重要统计区域；②重要统计区域CIS几何尺寸过小，造成厚屏蔽微探测问题；③依赖粒子密度先验信息，优化迭代过程中几何尺寸与材料的不断更新，需要全局计算且耗时较长。基于此，本研究提出了一种特征插值权窗设置方法，以提高迭代计算过程中的计算精度和效率。传统先验权窗是计算得到全局通量率分布并令权窗上限W_{\mathrm{u}}和下限W_{\mathrm{l}}为

W_{\mathrm{u}}=\frac{\varphi }{{\varphi }_{\mathrm{m}\mathrm{a}\mathrm{x}}}W_{\mathrm{l}}=\frac{W_{\mathrm{u}}}{N_{\mathrm{m}\mathrm{a}\mathrm{x}}}（14）

其中，\varphi、{\varphi }_{\mathrm{m}\mathrm{a}\mathrm{x}}分别为当前位置的通量率值与全局最大通量率值。在迭代优化时，更新先验数据权窗效率低，采用固定先验权窗无法适配结构参数更新带来的通量分布变化，导致减方差效果会下降。考虑到对称性，在x负半轴上均匀选取n个特征点，将其原权窗值作为特征权重值，其表达式如下：

(u_i,W_i)=(\alpha x_i,W_{\mathrm{c},i}^{\beta })i=\mathrm{0,1},\dots ,n（15）

其中，\alpha、\beta分别为几何放缩因子和权重放缩因子，x_i为采样点坐标值，W_{\mathrm{c},i}为特征点权重值。为得到平滑的权重分布函数，对特征点进行三次样条插值以获取平滑曲线，插值函数可定义为^［19］

\begin{array}{c}A(u)=\frac{A(u_{i-\mathrm{1}})(u_i{-u)}^{\mathrm{3}}}{\mathrm{6}{\lambda }_{i-\mathrm{1}}}+\frac{A(u_i)(u-u_{i-\mathrm{1}}{)}^{\mathrm{3}}}{\mathrm{6}{\lambda }_{i-\mathrm{1}}}+\\ \frac{(W_{i-\mathrm{1}}-\frac{A(u_{i-\mathrm{1}}){\lambda }_{i-\mathrm{1}}^{\mathrm{2}}}{\mathrm{6}})(u_i-u)}{{\lambda }_{i-\mathrm{1}}}+\end{array}

\frac{(W_i-\frac{A(u_i){\lambda }_{i-\mathrm{1}}^{\mathrm{2}}}{\mathrm{6}})(u-u_{i-\mathrm{1}})}{{\lambda }_{i-\mathrm{1}}}（16）

其中，u、{\lambda }_{i-\mathrm{1}}为参数，{\lambda }_{i-\mathrm{1}}=u_i-u_{i-\mathrm{1}}。点u_i处的函数具有连续二阶导数，添加如下自由边界条件：

A^{″}(u_{\mathrm{0}})=W_{\mathrm{0}}^{″}=\mathrm{0}A^{″}(u_n)=W_n^{″}=\mathrm{0}（17）

{\mu }_i{A}^{\text{'}}(u_{i-\mathrm{1}})+\mathrm{2}{A}^{\text{'}}(u_i)+{\lambda }_i{A}^{\text{'}}(u_{i+\mathrm{1}})={\delta }_i（18）

\left.\begin{array}{l}{\mu }_i=\frac{{\lambda }_{i-\mathrm{1}}}{{\lambda }_{i-\mathrm{1}}+{\lambda }_i}\\ {\lambda }_i=\mathrm{1}-{\mu }_i\\ {\delta }_i=\frac{\mathrm{6}(\frac{W_{i+\mathrm{1}}-W_i}{{\lambda }_i}-\frac{W_i-W_{i-\mathrm{1}}}{{\lambda }_{i-\mathrm{1}}})}{{\lambda }_i+{\lambda }_{i+\mathrm{1}}}\end{array}\right\}（19）

其中，{\mu }_{\mathrm{0}}=\mathrm{0},{\delta }_{\mathrm{0}}={\delta }_n=\mathrm{0},{\lambda }_n=\mathrm{0}，A^{\text{'}}\mathrm{、}{A}^{″}分别为插值函数的一阶与二阶导数。设权窗中心到原点的距离为u，将其代入式（16）~式（19）得到插值函数，并计算特征插值权窗权重值分布。

如图5和图6所示，与先验权窗相比，特征插值权窗构建了逆梯度方向指向几何中心的平滑权窗分布，该分布以类漏斗式引导粒子向中心核心统计区域输运并分裂扩增，以减小区域方差。权窗效果常用品质因子V_FOM（figure of merit，FOM）评判计算效率，其定义为^［20］

V_{\mathrm{F}\mathrm{O}\mathrm{M}}=\frac{\mathrm{1}}{D^{\mathrm{2}}t}（20）

其中，D为统计相对误差，t为计算用时。为研究采样点数对FOM值的影响，分别测试 4 种不同采样数权窗下各50组随机模型的性能。

图7展示了数据分布特征，其中灰色箱体代表25%~75%分位数区间即四分位距（interquartile range，IQR），须线为1.5倍IQR内的正常数据范围。当采样数n=4时权窗性能较差，随着采样数的增加，FOM值得到了增大，但n>6后FOM值提升效果有限，n=10时权窗性能最佳。故令n=10，测试不同权窗的计算效率并分析，如图8~图10所示。由图 8可见，特征插值权窗的FOM值提升率为300%~700%，计算效率得到大幅提高，且稳定性优于先验权窗。图9表明两种权窗均有效减小了相对误差且效果相近。图 10显示先验数据权窗的用时不稳定且有激增，特征插值权窗的用时少且稳定。

4 屏蔽结构多目标优化

本研究选用遗传算法^［21］开展多目标优化设计，设计参数如图11所示，尺寸参数x_{\mathrm{1}}~x_{\mathrm{6}}和材料编号m_{\mathrm{1}}~m_{\mathrm{4}}分别采取实数和整数编码策略。表1列出了备选材料，其中m_{\mathrm{1}}\mathrm{、}{m}_{\mathrm{2}}只在前4种透光材料中选择。本文采用适应度函数法^［14］，通过权重设置将式（1）中多目标优化转换为单目标优化，得到目标函数为

F(\mathit{X})=C_R\frac{F_R(\mathit{X})}{R_{\mathrm{0}}}+C_V\frac{F_V(\mathit{X})}{V_{\mathrm{0}}}+C_M\frac{F_M(\mathit{X})}{M_{\mathrm{0}}}（21）

式中：C_R\mathrm{、}{C}_V\mathrm{、}{C}_M分别为总剂量率、总体积和总质量的加权系数，且C_R+C_V+C_W=\mathrm{1}，并施加约束C_R>C_V\ge C_W。

MC仿真中辐射源输入参考自反应堆压力容器外表面辐射通量率数据^［22］，中子通量率为\mathrm{9.23}\times {\mathrm{10}}^{\mathrm{7}}\mathrm{c}{\mathrm{m}}^{-\mathrm{2}}\cdot {\mathrm{s}}^{-\mathrm{1}}，主要为快和慢中子；光子通量率为\mathrm{4.15}\times {\mathrm{10}}^{\mathrm{8}}\mathrm{c}{\mathrm{m}}^{-\mathrm{2}}\cdot {\mathrm{s}}^{-\mathrm{1}}，主要为伽马射线。C_R\mathrm{、}{C}_V\mathrm{、}{C}_M的具体取值见表2。

目标函数优化的结果如图12所示，可得目标函数在90代后均趋于收敛。剂量率、总体积和总质量的优化结果分别如图13~图15所示。图13中，剂量率加权系数高，在迭代过程中收敛于较低值。在图14和图15中，由于受小体积、小质量与低剂量率之间矛盾关系的影响，曲线在前期持续振荡，后期趋于稳定。结合图13和表3可知，第4组中目标函数值最小和剂量率第2小，较无屏蔽时的辐射剂量率0.562 Gy/h降低了91.6%，与预设计的铅与铅玻璃结构（初始值：剂量率0.157 Gy/h，体积5.92×10³ cm³，质量64.6 kg）相比，体积减小11%，质量减小66%。结合表2和图13~图15可知，第1组中剂量率加权系数最高，因此辐射剂量率最低，剂量率下降了93%，但其体积与预设计相当，质量减小了22%。图16展示了各组最优解中各目标相对于预设计初始值的比例。

结合表3分析材料组成共性：透镜材料均选用铅玻璃，表明在综合因素下铅玻璃是中-光子屏蔽的优良材料；石英玻璃使用频率高，其密度小、高光谱透射率及抗辐射性能佳，且长时间暴露于辐射不易变暗^［23］；包壳均选择钨，中-光子综合屏蔽性能优异；PEB和PE的选择比例高，其质地轻，氢原子含量高，对中子慢化屏蔽作用好，被广泛应用于核工业相关辐射屏蔽^［24］。

此外通过计算加权系数C_R\mathrm{、}{C}_V\mathrm{、}{C}_M与优化程度之间的皮尔逊相关系数（表4），以此评估加权系数对各目标结果的影响，进而指导加权系数的设置。由表4可得，剂量率与剂量率加权系数成较强负相关，与体积加权系数成较强正相关，与质量加权系数成较弱正相关。这意味着，若要减小剂量率，可适当增大剂量率加权系数或减小体积加权系数。此外，体积与剂量率加权系数、体积加权系数均成强负相关，故要减小体积，可适当增大剂量率及体积加权系数。

除各组最优解外，GA的全局搜索特性能获取Pareto前沿，是非劣解在多目标之间达到的一个平衡态^［25］，可为屏蔽结构的设计提供更多样化的选择方案。迭代后各组非劣解如图 17~图 22所示，这些非劣解在各维度上呈条带状分布，边界构成Pareto前沿，可为决策提供在不同目标之间进行平衡的丰富选择，对后续指导设计方案的优化、提高产品性能、降低成本等方面具有重要意义。

5 结论

针对核环境下互补金属氧化物半导体图像传感器（CIS）的辐射屏蔽问题，本文提出了一种基于蒙特卡罗（MC）减方差方法的辐射仿真和屏蔽结构多目标优化设计方案，得到如下结论：

1）根据辐射作用原理设计了CIS复合辐射屏蔽结构模型，基于其厚屏蔽微探测特性，提出球壳源与权窗减方差仿真方案并使用MC方法进行辐射屏蔽仿真分析，实现中-光子耦合输运计算，反映屏蔽结构和辐射源之间复杂相互作用。

2）为解决迭代过程中实时更新权窗效率低下的问题，提出一种特征插值权窗MC减方差方法。通过对比若干不同权窗计算结果，发现特征插值权窗高效稳定，可提高优化过程整体计算效率。

3）构建CIS辐射屏蔽结构的参数化模型，并通过遗传算法实现多目标参数优化，得到一系列参数组合下的最优解与Pareto前沿，为决策者提供了不同目标间平衡的丰富选择，为CIS辐射屏蔽结构的多目标优化设计问题提供了一套解决方案，对核工业装备的优化设计有一定参考价值。

参考文献

[1] CHANG Y， ZOU S， LIN G， et al. Research on Radiation Damage and Reinforcement of Control and Sensing Systems in Nuclear Robots［J］. Electronics， 2024， 13（7）：1214.

[2] KHANAM Z， ASLAM B， SAHA S， et al. Gamma-Induced Image Degradation Analysis of Robot Vision Sensor for Autonomous Inspection of Nuclear Sites［J］. IEEE Sensors Journal， 2021， 22（18）： 17378-17390.

[3] JOHNSTON A H， RAX B G， LEE C I. Enhanced Damage in Linear Bipolar Integrated Circuits at Low Dose Rate［J］. IEEE Transactions on Nuclear Science， 1995， 42（6）： 1650-1659.

[4] ANTONSANTI A， VIRMONTOIS C， LAUENSTEIN J M， et al. Probing Dark Current Random Telegraph Signal in a Small Pitch Vertically Pinned Photodiode CMOS Image Sensor after Proton Irradiation［J］. IEEE Transactions on Nuclear Science， 2022， 69（7）： 1506-1514.

[5] GAO T， YIN C， CHEN Y， et al. Simulation of Total Ionizing Dose Effects Technique for CMOS Inverter Circuit［J］. Micromachines， 2023， 14（7）： 1438.

[6] NÁFRÁDI G， CZIFRUS S， KOCSIS G， et al. Analysis of Dark Current Images of a CMOS Camera during Gamma Irradiation［J］. Fusion Engineering and Design， 2013， 88（12）： 3169-3175.

[7] XUE Y， WANG Z， CHEN W， et al. Modeling Dark Signal of CMOS Image Sensors Irradiated by Reactor Neutron Using Monte Carlo Method［J］. Science China Information Sciences， 2018， 61： 1-10.

[8] 徐守龙， 邹树梁， 彭聪. CMOS APS数字模组辐射环境中的屏蔽加固［J］. 核动力工程， 2018， 39（4）： 165-170.

[9] XU Shoulong， ZOU Shuliang， PENG Cong. Shielding Reinforcement Of CMOS APS Digital Modules in Radiation Environment［J］. Nuclear Power Engineering， 2018， 39（4）： 165-170.

[10] SONG Y， ZHANG Z， MAO J， et al. Research on Fast Intelligence Multi-objective Optimization Method of Nuclear Reactor Radiation Shielding［J］. Annals of Nuclear Energy， 2020， 149： 107771.

[11] CHEN Z， ZHANG Z， XIE J， et al. Multi-objective Optimization Strategies for Radiation Shielding Design with Genetic Algorithm［J］. Computer Physics Communications， 2021， 260： 107267.

[12] 许淑艳. 蒙特卡罗方法在实验核物理中的应用［M］. 北京：原子能出版社， 2006： 19-67.

[13] XU Shuyan. Application of Monte Carlo Method in Experimental Nuclear Physics［M］.Beijing： Atomic Energy Press， 2006： 19-67.

[14] BINJOLA A. Interaction of Radiation with Matter［M］∥ Practical Radiation Oncology.Singapore： Springer， 2020： 3-11.

[15] HU G， HU H， YANG Q， et al. Study on the Design and Experimental Verification of Multilayer Radiation Shield against Mixed Neutrons and γ-Rays［J］. Nuclear Engineering and Technology， 2020， 52（1）： 178-184.

[16] BROWN F B， BARRETT R F， BOOTH T E， et al. MCNP Version 5［J］. Transactions of the American Nuclear Society， 2002， 87（273）： 02-3935.

[17] ROMANO P K， JOSEY C J， JOHNSON A E，et al. Depletion Capabilities in the Openmc Monte Carlo Particle Transport Code［J］. Annals of Nuclear Energy， 2021， 152： 107989.

[18] ALSAYED Z， BADAWI M S， AWAD R， et al. Study of Some γ-Ray Attenuation Parameters for New Shielding Materials Composed of Nano ZnO Blended with High Density Polyethylene［J］. Nuclear Technology and Radiation Protection， 2019， 34（4）： 342-352.

[19] BARTHEL J， SARIGUL-KLIJN N. Importance of Spherical Shell Models for Radiation Shielding Designs on Space Missions［J］. Journal of Spacecraft and Rockets， 2019， 56（5）： 1658-1661.

[20] PAN Q， RAO J， HUANG S， et al. Improved Adaptive Variance Reduction Algorithm Based on RMC Code for Deep Penetration Problems［J］. Annals of Nuclear Energy， 2020， 137： 107113.

[21] AKIMA H. A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures［J］. Journal of the Association for Computing Machinery， 1970， 17（4）：589-602.

[22] DAVIS A， TURNER A. Comparison of Global Variance Reduction Techniques for Monte Carlo Radiation Transport Simulations of ITER［J］. Fusion Engineering and Design， 2011， 86（9/11）： 2698-2700.

[23] 姜一啸， 吉卫喜， 何鑫， 等.基于改进非支配排序遗传算法的多目标柔性作业车间低碳调度［J］. 中国机械工程， 2022，33（21）： 2564-2577.

[24] JIANG Yixiao， JI Weixi， HE Xin， et al. Low-carbon Scheduling of Multi-objective Flexible Job Shop Based on Improved Non-dominated Sorting Genetic Algorithm［J］. China Mechanical Engineering， 2022， 33（21）： 2564-2577.

[25] 陈德锋. 基于MCNP程序的AP1000反应堆屏蔽分析研究［D］. 哈尔滨： 哈尔滨工程大学， 2013.

[26] CHEN Defeng. Research on Shielding Analysis of AP1000 Reactor Based on MCNP Program［D］. Harbin： Harbin Engineering University， 2013.

[27] AKHTAR S M J， ASHRAF M， KHAN S H. A Study of Neutron and Gamma Radiation Effects on Transmission of Various Types of Glasses， Optical Coatings， Cemented Optics and Fiber［J］. Optical Materials， 2007， 29（12）： 1595-1603.

[28] ZHANG S， WANG Z， LI C， et al. Study on Neutron Shielding Performance of Hot Cell Shielding Door for Nuclear Power Plant［J］. Annals of Nuclear Energy， 2022， 166： 108752.

[29] 田宗睿， 智鹏鹏， 云国丽， 等. 基于自适应增量Kriging模型的多目标稳健优化设计方法［J］. 中国机械工程， 2023， 34（8）： 931-940.

[30] TIAN Zongrui， ZHI Pengpeng， YUN Guoli， et al. Multi-objective Robust Optimization Design Method Based on Adaptive Incremental Kriging Model［J］. China Mechanical Engineering， 2023， 34（8）： 931-940.

基金资助

国家自然科学基金(12402233)

国家自然科学基金重点项目(52235005)