基于统计抽样的敏感性系数计算方法

马续波，唐 辉，杨 乐，朱帅涛，陈义学

(华北电力大学 核科学与工程学院，北京 102206)

不确定性和敏感性分析在各领域都有重要地位[1-3]，如对物理模型输入参数的敏感性分析可知哪个参数的重要性更大，确定出问题的主要因素，以方便构造更好的模型。在核工程领域，为保证核电厂的绝对安全，以往的计算模型通常采用留有很大余量的保守模型进行设计，而进行反应堆工程设计中的关键参数的不确定性和敏感性分析方法可较好理解核电厂设计的保守性以及改进核电厂的经济性。以往的不确定性和敏感性分析方法通常有两种，分别是确定论方法和统计抽样方法，但基于确定论方法的不确定性分析方法通常仅能考虑低阶的不确定性(一般能到二阶，高阶较难考虑)，并且对分析的响应量有一定要求。所以发展了广义微扰理论的敏感性计算方法，以可更好考虑响应量的扩展量，如控制棒价值、堆芯功率等。而基于统计抽样方法的不确定性分析方法，其优点是对响应量无要求，且易计算响应量的不确定度，其缺点是较难计算响应量输入参数的敏感性系数[4-5]以及耗费的计时较长。伴随计算机计算能力的发展，计算时间较长问题也在逐步得到改善。本文推导基于统计抽样方法计算响应量对输入参数的敏感性系数的理论公式，然后利用算例进行验证分析，验证理论方法的可行性。

1 敏感性系数计算理论

假设响应量R为N个输入变量的函数：

R=R(α1,…,αN)

(1)

而N个输入变量的真值是由变量的最佳估计值和相应的不确定度组成的：

(2)

对响应量进行泰勒展开，且在一阶近似下可得：

(3)

通常定义输入变量αi相对于响应量R的敏感性系数为：

(4)

如果计算得到了响应量R对每个参数的敏感性系数，则可根据误差传递公式(式(5))计算响应量R的协方差，进而利用式(6)求得响应量R的标准偏差ΔR，即为R的不确定度[3]。

var(R)=SVαST

(5)

ΔR=[var(R)]1/2

(6)

其中：S为敏感性系数向量，是N个参数组成的行向量；ST为向量S的转置向量；Vα为N×N的输入变量的协方差矩阵。

通过上面的分析可知：求ΔR需先求得响应量R对每个输入参数αi的敏感性系数矩阵S，根据输入变量之间的协方差矩阵可求得ΔR。而基于统计抽样方法的不确定性分析方法的思路与上面有所不同，其基本思想为：先根据输入变量α中每个参数的分布规律以及每个参数之间的协方差矩阵Vα进行抽样，然后利用式(1)可得到响应量R的样本，对样本进行统计分析即可得到ΔR。基于统计抽样方法的优点在于：1) 不仅可考虑输入参数对R的一阶扰动，所有阶的扰动均可考虑；2) 对响应量无限制，方法的通用性更好。但其缺点也很明显，即不易求出输入变量αi相对于响应量R的敏感性系数。在核工程中，由于离散能群的不等宽效应，可能产生因截面所属能量宽度较大引起的较大敏感性系数[6]，为此，本文引入平均勒能量相对敏感性系数：

(7)

(8)

式中：i为样本编号；j为输入变量即能群截面的编号。上式可变为：

(9)

(10)

上式可写为：

(11)

利用上式，可求得响应量对每个能群截面的敏感性系数：

(12)

2 测试和验证

2.1 双群临界公式验证

为验证基于统计抽样的敏感性系数计算方法的可行性，本文提出了利用裸堆双群近似下的临界公式进行验证。裸堆双群近似的临界公式如式(13)所示。

(13)

其中：(υΣf)1、(υΣf)2分别为第1和第2能群的产生截面；D1、D2分别为第1和第2能群的扩散系数；Σ1→2为第1能群到第2能群的散射截面；Σa,1、Σa,2分别为第1、2能群的吸收截面；Σr=Σa,1+Σ1→2，为输运截面；B2为几何曲率常数，假定几何曲率常数B2=2.9×10-4。

选取了式中的7个参数进行研究，7个参数的编号及初始值列于表1。假定每个参数的相对误差为1%，且服从正态分布。假定7个参数之间是相互独立的，不考虑参数之间的相关性情况下，表2列出了采用不同的样本总数的敏感性系数的结果，表2中的基准值是指直接利用式(4)的计算结果。由表2的结果可见：基于抽样方法的敏感性系数计算方法是正确的，且对于敏感性系数较大的参数，随着抽样样本数的增加，计算结果与基准值吻合越好，但对于敏感性系数很小的参数，如D2和Σr，由于统计性原因导致误差较大。图1a示出了敏感性系数的相对误差随样本总数的变化。由图1a可见，样本数达到2 000时，除了两个敏感性系数较小的参数，其他各主要的误差基本可控制在5%以内。这也说明了该方法的可行性。

表1 裸堆双群近似的临界公式参数编号及初始值Table 1 Parameter number and initial value of two groups multiplication factor of bare homogeneous ractor

基于抽样的敏感性系数计算方法的1个优点是可方便考虑研究参数之间的相关性。本文研究了前3个参数的相关性对敏感性系数计算的影响。假定前3个参数的相关性系数矩阵为Σ，有：

表2 不同抽样样本数对结果影响Table 2 Result effects with respect to sampling number

a——敏感性系数的相对误差随抽样样本数的变化关系；b——(υΣf)1、(υΣf)2和Σ1→2之间的相关性系数对敏感性系数计算的影响图1 抽样样本数和参数之间相关性对敏感性系数的影响Fig.1 Sensitivity coefficients variety with total sampling sizes and correlation between parameters

(14)

其中，a为两个参数的相关性的大小。图1b示出了设定的相关性系数与敏感性系数的关系(样本数为10 000)，由图1b可见，在相关性系数较小的情况下，相关性系数对3个参数的敏感性系数的计算结果影响不大，但在相关性系数较大的情况下，有可能导致相关性系数出现跳跃现象，且在后面的压水堆单栅元的敏感性分析中也出现了类似现象。

2.2 压水堆单栅元问题研究

核工程设计中，通过对复杂系统的不确定性分析可减小系统的近似程度，从而提高核设计的经济性。通过对系统的敏感性系数分析，可进一步明确输入变量的重要性，进而为提高输入参数精度指明方向。近年来，国内外对此也进行了大量研究[7-17]。为研究核工程计算过程中的不确定性，OECD/NEA专门提出了基准题，国际上所有研究单位可利用此基准题对比校核各自程序的适用性[9]。本文采用基准题中的三哩岛(TMI)燃料栅元为研究对象，研究利用统计抽样方法计算235U裂变截面的敏感性系数。为解决式中的协方差求逆问题，采用两种方法解决此问题。方法1是去掉协方差矩阵中的相关项，只保留主对角线部分进行求逆；方法2是假定每个能群的截面的误差均为0.1%或1%的协方差矩阵。

(15)

式中，V′σ,σ为简化后的协方差矩阵。

235U裂变截面的协方差矩阵采用69群能群结构，基于ENDF/B-Ⅶ.1核评价数据库，经NJOY程序加工而成[18]。235U裂变截面的相关性系数矩阵如图2所示。由图2可知，部分能群之间存在很强的相关性。

为验证不同计算方法的正确性，开发不确定性与敏感性分析程序SUACL，该程序目前主要具备压水堆开展截面的不确定性和敏感性系数分析功能。计算中对TMI的输运计算采用了加拿大蒙特利尔大学开发的反应堆计算模拟软件DRAGON4.0。采用方法1或方法2在计算敏感性系数时均采用式(15)，不再采用式(12)，因为产生样本时的协方差矩阵已用V′σ,σ替换了Vσ,σ，由于V′σ,σ很易求出其逆矩阵，因此可很方便求出。图3示出了方法1计算结果与直接微扰计算结果对比，图中的Full matrix表示采用图2给出的协方差矩阵，抽样样本为500时计算得到的相对敏感性系数；Diagonal matrix表示对图2给出的协方差矩阵去掉了相关性系数，仅保留了主对角线部分，即方法1，计算得到的相对敏感性系数；Direct perturbation表示直接利用式(4)计算得到的结果。由图3可见，采用方法1，即主对角元矩阵，而不考虑相关性，计算结果与直接微扰的结果吻合很好。这种结果很易理解，因为利用式(4)的直接微扰方法中并未考虑能群之间的相关性。而考虑能群之间的相关性系数后会出现震荡，通过之前的双群近似下的临界公式分析可知，震荡主要是由于考虑了能群之间的相关性导致。

图2 235U 69群裂变截面的相关性系数矩阵Fig.2 Correlation coefficient matrix of 69 groups fission cross section of 235U

图3 利用方法1计算的235U裂变截面的相对敏感性系数Fig.3 Relative sensitivity coefficient of 235U fission cross section with method 1

方法1中采用主对角线元素就可得到较好的相对敏感性系数，由此推断：如果每个能群均采用相同的扰动量，如均为1%或0.1%，那么也应能计算出较好的敏感性系数，为此采用方法2进行验证。图4示出了方法2中各能群截面的微扰量为1%和0.1%的计算结果与直接微扰法的比较，样本总数为500。由图4可见，在本例子中，微扰量1%和0.1%均可计算出235U裂变截面的相对敏感性系数，但微扰量为0.1%的计算结果较1%的微扰量的计算结果要好，即与直接微扰法相比，误差更好。

图4 利用方法2计算的235U裂变截面的相对敏感性系数Fig.4 Relative sensitivity cofficient of 235U fission cross section with method 2

根据敏感性系数可利用式(5)计算得到对应响应量的不确定度。表3列出了采用不同的敏感性系数计算方法下235U裂变截面对栅元kinf的相对不确定度及其误差。由表3可见，微扰量为0.1%时较微扰量为1%情况好，此结果也可由图4看出，采用Full matrix的情况，虽然敏感性系数跳跃性很强，但很多能群正负抵消，导致最后对栅元的kinf的相对不确定度的误差仅1.36%。

表3 不同敏感性系数计算方法下kinf的不确定度Table 3 Uncertainty of kinf with different sensitivity cofficient calculation methods

3 结论

本文推导了基于抽样方法进行敏感性系数计算的理论公式，首先利用简单的例子进行了验证分析，验证了方法的可行性，同时也发现了该方法在计算敏感性系数较小的参数时结果可能会有较大误差。针对核工程中的单栅元基准题问题，利用此方法的难点为协方差矩阵的求逆，本文提出了两种替代的方法进行敏感性系数计算，从而避免了协方差矩阵求逆困难的问题，进而对所提出的方法进行了验证，验证了方法的可行性。本文还发现微扰量的不同会导致计算的敏感性系数的精度不同，微扰量的选取方法是下一步要研究的问题。