CN113221955B - 一种针对反应堆物理分析中高维输入参数的不确定性传播方法 - Google Patents

Abstract

本发明属于核反应堆堆芯设计和安全评估领域，涉及高维输入参数的不确定性传播方法，包括步骤如下：根据给定输入参数个数，读取输入参数及相应维度相对协方差矩阵；对该矩阵进行主成分分析，确定降维后的特征维度；在总体分布空间上进行等区间划分，确定每个区间上随机样本的累计概率密度；通过求解正态分布的反函数，计算相应累积概率密度对应的随机样本；计算随机样本协方差矩阵；对该矩阵进行奇异值分解，获得左奇异矩阵U，并构造对角矩阵D；利用左奇异矩阵U及对角矩阵D，构建新的随机样本空间；利用线性关系获得最终样本。本发明同时考虑统计相关性和高维两个关键问题，使得更适用于高维输入下的反应堆物理分析中的不确定性传播用途。

Description

一种针对反应堆物理分析中高维输入参数的不确定性传播 方法

技术领域

本发明涉及核反应堆堆芯设计和安全评估领域，具体涉及针对反应堆物理分析中高维输入参数的不确定性传播方法。

背景技术

随着高性能计算机技术的逐步发展，建模与仿真技术逐渐成为研究者模拟客观物理现象并开展进一步分析、设计的主要手段。然而，由于人们对实际物理现象认知不明确以及对模拟系统输入参数测量不可避免地存在误差，因此不确定性天然存在于建模与仿真过程中，这会使得最终预测结果也带有一定的不确定性，从而影响模拟结果的可信度，因此预测结果的不确定性必须被量化出来。要准确量化不确定性，必须保证输入参数的不确定性在建模与仿真过程中被合理传播，即既能保证输入参数的不确定性信息被完整保留，也能确保参数间的相关性不被破坏，同时为了尽量减少计算代价，处理具有高维输入参数的分析过程也是亟待解决的问题。

目前可以用于不确定性传播的技术有简单随机采样方法(SRS)、重要性采样方法(IS)及拉丁超立方体采样方法(LHS)；为了解决高维问题而产生的降维技术可分为线性降维和非线性降维，线性降维主要有：主成分分析(PCA)、独立成分分析(ICA)等；非线性降维主要有：核主成分分析法(KPCA)、等度量隐射(ISOMAP)等。上述不确定性传播技术均属于抽样方法，其原理为从N个总体中逐个不放回的提取n个样本个体(n≤N)，且保证每个个体被提取到的机会相等。提取过程中不考虑总体的分布特性而仅提取个体的方法为简单随机采样方法(SRS)；从与原分布不同的另一个分布中提取个体而对原先分布的性质进行估计的方法为重要性采样方法(IS)；将总体分成互不相交的层，然后按照等数分配原则，从各层独立地提取个体的方法为拉丁超立方体采样方法(LHS)。

对于不确定性传播方法来说，简单随机采样方法(SRS)、重要性采样方法(IS)及拉丁超立方体采样方法(LHS)三者共同的缺点为：在获取样本个体的过程中不可避免地会引入统计相关性，使得最终提取到的样本无法真实描述输入参数之间的相关性，使得不确定性传播出现偏差。对于简单随机采样方法(SRS)而言，通过此方法提取到的样本不能很好地描述总体分布特征，因此必然会使不确定性的传播出现信息遗漏。上述各方法中的降维方法分别介绍如下。主成分分析(PCA)通过计算具有n维特征的数据协方差矩阵，得到协方差矩阵的特征值和特征向量，选择特征值最大(即方差最大)的k个特征所对应的特征向量组成的矩阵从而实现数据特征的降维；独立成分分析(ICA)是在子分量为非高斯分布且相互独立的假设下，将多元信号分离为加性子分量的以实现降维；核主成分分析法(KPCA)是通过构造核函数，将原空间上线性不可分的样本变为线性可分，再使用主成分分分析从而实现降维；等度量隐射(ISOMAP)的原理为将计算测地线距离的问题转化为计算近邻连接的两点间最短路径问题，从而实现样本点在低维空间中的坐标以实现降维。

其中，对于降维技术而言：独立成分分析(ICA)处理的数据为非高斯分布，然而对于反应堆物理分析而言，其输入参数往往被认为是高斯分布；同时，ICA是从总体中找到最大独立化的分量，这也与反应堆物理计算对输入参数信息进行降维的目的相违背；核主成分分析法(KPCA)虽然可以最大限度的提取样本的信息，但是其提取过程的实际意义并不明确，导致物理意义无法解释，同时由于该方法需要进行核函数的计算，因此计算量较大；等度量隐射(ISOMAP)常用在图像识别和运动追踪中，对于针对核反应堆物理分析中的不确定性传播不具有适用性。

对于以上不确定性传播方法会不可避免地引入统计相关性的问题，这是由于在实现样本提取的过程中会用到随机数，而计算机产生的随机数并不能保证不同序列间完全独立，从而使得引入的统计相关性传递至最终样本。对样本不能很好地描述总体分布特征问题，这是因为在采用简单随机采样方法(SRS)提取个体时对于个体的选择是完全随机的，因此在个体数较少时很容易出现在总体某一区域内集中的问题，从而导致分布特性失真。

对于降维技术而言，其主要缺陷表现为在核反应堆物理分析中不确定性传播领域的不适用性，这是由于其方法本身的限制造成的。

对于因获取样本个体的过程中引入统计相关性并使得不确定性传播出现偏差的问题，从业人员往往会通过多次采样以量化这样的偏差，并在分析中采用定性分析的方法来考量此问题；对于提取到的样本不能很好地描述总体分布特征的问题，从业人员会通过采用对总体分层并从各层中提取样本个体的方法解决此问题，如采用拉丁超立方体抽样(LHS)方法，但此方法同样面临不确定性传播出现偏差的问题。对于高维问题，从业者往往会提取高于输入参数维度的样本掩盖此问题，但会使得计算代价显著增加。

对于上述两个问题，从业者一般仅针对单个问题进行处理，因此二者共同作用产生的效应往往被忽略。同时，由于使用采样方法传播不确定性本身就具有一定的不确定性，因此若非对传播过程进行细致的统计学分析，问题本身就很难被发现。对于高维问题带来的计算量大的问题，从业者一般采用投入大量计算量的方法，因此该问题很容易被掩盖。

发明内容

对于因获取样本个体的过程中引入统计相关性的问题，本发明采取在过程中对随机样本空间的协方差矩阵进行奇异值分解，利用左奇异矩阵及构造矩阵得到新的随机样本空间的方法，这样即可在保证样本总体分布类型的同时，也能使得样本间的相关性被完全保留；对于高维问题，本发明采取对输入参数的协方差矩阵进行主成分分析的方法，获得能够描述原协方差矩阵信息的维度更低的协方差矩阵，从而达到减少样本个体数目的目的。本发明同时考虑以上两个关键问题，使得更适用于高维输入下的反应堆物理分析中的不确定性传播用途。

本发明提供一种针对反应堆物理分析中高维输入参数的不确定性传播方法，包括步骤如下：

步骤1：根据给定输入参数个数，读取输入参数及相应维度相对协方差矩阵；

步骤2：对读取的相对协方差矩阵进行主成分分析，确定降维后的特征维度；

步骤3：在总体分布空间上进行等区间划分，确定每个区间上随机样本的累计概率密度；

步骤4：通过求解正态分布的反函数，计算相应累积概率密度对应的随机样本；

步骤5：计算随机样本协方差矩阵；

步骤6：对随机样本空间协方差矩阵进行奇异值分解，获得左奇异矩阵U，并构造对角矩阵D；

步骤7：利用左奇异矩阵U及对角矩阵D，构建新的随机样本空间；

步骤8：利用线性关系获得最终样本。

根据上述方法可以通过编制计算机程序来实现。因此本发明还涉及一种针对反应堆物理分析中高维输入参数的不确定性传播方法的计算流程及其应系统，其包括对应各个步骤的模块，即参数及相应维度相对协方差矩阵输入和读取模块；相对协方差矩阵进行主成分分析模块；随机样本的累计概率密度计算模块；累积概率密度对应的随机样本计算模块；随机样本协方差矩阵模块计算模块；随机样本空间协方差矩阵奇异值分解模块；构建新的随机样本空间模块；最终样本获取模块。上述可以通过计算程序存储在存储器中，由于处理器运行该程序来实现所述的方法。

本发明通过主成分分析，可以确定特定问题的最小样本数目，从而减少采样次数和样本个体数目，达到减少计算量的目的；通过采用奇异值分解技术，可以消除随机抽样过程中引入的统计相关性，确保不确定性合理传播。以上两点均为直接效果，在二者的共同作用下会产生以下深层次效果：使得采样过程变得更加快速高效；在使得采样次数减少的情况下，由于减少了采样过程中引入的统计相关性，因此使得不确定性传播更加准确，显著降低后续分析的误差。

本发明的上述方法的研发获得过程中，发明人也进行了更多的尝试。列举一二如下：

一是为了解决高维抽样问题，曾经提出了结合进化算法的拉丁超立方体抽样(LHS-EA)方法，其原理为：由于样本相关系数矩阵C为对称非正定矩阵，目标函数如下所示：

式(3)中c_ij为相关系数矩阵C的第i行、第j列的相关系数。结合进化算法的拉丁超立方体抽样方法的原理是：使用LHS抽取样本空间后，若目标函数c_s的值越接近0，则参数间的相关性也就越小，最后将随机排序后能收敛到最小c_s的序列保留。但该方法存在显著缺陷，即：进化次数的选取需要人为设置，并且当进化次数设置较大时，程序运算量巨大；使用该方法虽然降低了参数间相关性，但是参数间的相关性仍较大，不满足设定的相关性要求。

二是为了解决采样过程中引入统计相关性的问题是曾提出了结合楚列斯基分解的拉丁超立方体(LHS-CDC)采样技术，其原理为：

对LHS方法得到的随机样本空间Z，对其相关系数矩阵C进行如下处理：

C＝QQ^T (21)

对于随机样本空间Z来说(此时只考虑nS>nX的抽样情况)，相关系数矩阵C是对称正定矩阵，满足使用Cholesky分解的前提条件，对上三角矩阵Q求得逆矩阵Q^-1，利用逆矩阵Q^-1对随机样本空间Z进行一一变换，如式(22)所示。

Z^*＝Q^-1Z (22)

经过上述处理得到新的随机样本空间Z^*，通过统计分析得到其相关系数矩阵C^*，对数值结果分析可知，相关系数矩阵C^*比C更接近单位对角矩阵。该方法虽然针对采样过程中引入统计相关性的问题有所改善，但事实上并未完全解决该问题。

附图说明

图1本发明方法过程示意图。

图2输入参数原始相关系数矩阵。

图3传统技术与本发明改进的对比。

图4不确定性传播验证。

具体实施方式

下面通过具体实施方式对本发明进行详细的说明，以期更好的理解本发明，但不构成对本发明的限制。

实施例一

本发明的实现步骤如下：

步骤1：根据给定输入参数个数，读取输入参数本征值及相应维度相对协方差矩阵，相对协方差矩阵用于描述输入参数自身不确定性及参数间相关性，输入参数本征值用于最后步骤同相对扰动量作用得到最终样本；

步骤2：对读取的相对协方差矩阵进行主成分分析，确定降维后的特征维度；

对于一个样本空间Z_nS×nX，其协方差矩阵为∑_nX×nX，根据奇异值分解理论可以得到协方差矩阵可以表示为：

∑_nX×nX＝USU^T (1)

结合主成分分析原理，对于矩阵∑，矩阵S_nX×nX为奇异值矩阵，对奇异值进行分析处理，选取主要的t(t<nx)个奇异值构造新的奇异值矩阵S_1t×t来表示协方差矩阵∑_t×t，如下式所示：

∑_t×t＝U₁S₁U₁ ^T (2)

此时即完成了对协方差矩阵为∑_nX×nX的降维。

相比而言，进化算法也是实现高维抽样的一种方案，其原理为：选择目标函数如下所示：

在抽取样本空间后，若目标函数c_s的值越接近0，则参数间的相关性也就越小，每次进化将随机排序后能收敛到最小c_s的序列保留，通过设置进化次数，获取最后保留的序列。但该方法存在显著缺陷：实施过程运算量巨大；相关性降低效果并不显著。因此该方案并未被本发明最终采纳。

步骤3：在总体分布空间上进行等区间划分，在该区间内随机抽取(0,1)的随机数，确定每个区间上随机样本的累计概率密度，方法为：

p_i＝(i-1)/nS+r_u/nS (4)

步骤4：通过求解正态分布的反函数，计算相应累积概率密度对应的随机样本，方法为：

为第j个输入参数分布函数的反函数；

步骤5：计算随机样本协方差矩阵，协方差计算方法为：

Cov(x_i,x_j)＝E(x_ix_j)-E(x_i)E(x_j) (6)

其中x_i，x_j为随机样本序列，将计算得到协方差按照对应位置排列即可得到协方差矩阵；

步骤6：对随机样本空间协方差矩阵进行奇异值分解，获得左奇异矩阵U，并构造对角矩阵D，方法为：

基于步骤4得到的样本空间Z_s＝[Z_s1,Z_s2,....,Z_snX]^T，该样本空间的协方差矩阵∑_s的格式是：

参数Z_si和Z_sj由抽样样本得到的协方差估计值如式(8)所示：

将参数间的协方差表达式(8)代入到协方差矩阵∑_s中，得到协方差矩阵∑_s格式如下：

其中，nS为样本数量。H是一个n×n满矩阵，其元素均为1。对协方差矩阵∑_s进行奇异值分解，

∑_s＝USV^T (10)

由于协方差矩阵∑_s为对称正定矩阵，因此有U＝V，则式(10)可转换为如下形式：

U^T∑_sU＝S (11)

其中，S为对角阵且其元素为矩阵∑_s的奇异值，构造对角矩阵E和对角矩阵D。其格式分别如式(12)-(14)所示。

其中，s₁,s₂,....,s_nX为矩阵∑_s的奇异值，并且按照从大到小的顺序排列。

步骤7：利用左奇异矩阵U及对角矩阵D，构建新的随机样本空间；

进而可得：

(D^-1)^TSD^-1＝E (15)

结合公式(11)可得：

(D^-1)^TU^T∑_sUD^-1＝E (16)

将式(9)带入公式(16)可得：

由公式(17)可知，对于样本空间

其协方差矩阵

为单位对角矩阵，表明样本空间

中的各参数间完全相互独立，

即为构建的新样本空间。

通过对随机样本的相关系数矩阵进行楚列斯基分解，并进一步构建新的随机样本也是实现相关性控制的可行方案，其原理为：

对步骤4得到的随机样本空间Z，对其相关系数矩阵C进行楚列斯基分解如下：

C＝QQ^T (18)

对于随机样本空间Z而言，由于此时已经进行了降维处理，因此相关系数矩阵C是对称正定矩阵，满足使用Cholesky分解的前提条件，对上三角矩阵Q求得逆矩阵Q^-1，利用逆矩阵Q^-1对随机样本空间Z进行变换，如式(19)所示。

经过上述处理得到新的随机样本空间

通过统计分析得到其相关系数矩阵C^*，统计可知相关系数矩阵C^*比C更接近单位对角矩阵。该方法虽然针对采样过程中引入统计相关性的问题有所改善，但事实上并未完全解决该问题。

步骤8：利用线性关系获得最终样本，方法为：

其中μ为输入参数本征值，A为对输入参数协方差矩阵采用奇异值分解获得的左奇异矩阵。

应用实施例

以一个44×44的核截面不确定性传播问题为例。对于一个反应堆物理计算的44能群问题，其输入参数为44能群的核截面，因此表征其不确定性的协方差矩阵为44×44的矩阵，即44×44维问题。该问题的输入参数相关系数矩阵如图2所示，该图反映了每一能群对应截面与其它能群截面的皮尔逊相关系数。按照实施例一的步骤进行验证：

步骤1：输入参数个数为44个，读取输入参数及相应44×44维度相对协方差矩阵；其中，44群能群结构下，每一能群对应的²³⁵U辐射俘获反应截面，从核截面库和协方差数据库中可以得到它的本征值和相对协方差并作为本发明的计算输入。

步骤2：采用本发明方案对读取的相对协方差矩阵进行主成分分析并确定了降维后的特征维度为18，因此本方案样本空间个数大于18即可，此处选择样本空间个数为30个，作为对比方案的LHS方法样本空间设置为100个；

步骤3：在总体分布空间上进行等区间划分，确定每个区间上随机样本的累计概率密度；

步骤4：通过求解正态分布的反函数，计算相应累积概率密度对应的随机样本；

步骤5：计算随机样本协方差矩阵；

步骤6：采用本发明的方案，对随机样本空间协方差矩阵进行奇异值分解，获得左奇异矩阵U，并构造对角矩阵D；

步骤7：利用左奇异矩阵U及对角矩阵D，构建新的随机样本空间；

步骤8：利用线性关系获得最终样本后，根据样本进行相关性保留验证及不确定性传播验证。

其中，作为对比方案的LHS方法无上述步骤4-步骤7。

上述验证是为了表明本发明在实现降维的同时能保证输入参数的相关性被完整保留，采用LHS方案作为对比，验证内容包括：(1)本发明能否实现降维；(2)本发明在降维条件下，与LHS方案相比是否保证了参数的相关性保留；(3)输入参数不确定性信息是否完整保留。

为了展示本发明的方法，对以上提到的不确定性传播方法具有的缺陷所进行的改善，以拉丁超立方体抽样(LHS)方法进行采样并作为计算对比验证，同时采用本技术(下文称改进技术)进行采样。正常情况下，样本个体的数目应当大于输入参数个数，对于本不确定性传播问题即为44，由于本发明可以通过降维技术降低样本个体的数目，所以在计算时只提取了30个样本个体。图2为输入参数初始相关系数矩阵，作为参数相关性保留的验证标准，通过拉丁超立方体抽样(LHS)方法及改进技术获得的参数相关系数矩阵对比，如图3所示。从图3表示的相关系数柱状图形状即可发现，采用LHS技术获得的参数样本相关性系数与初始相关性矩阵存在较大偏差，而采用本发明获得的参数样本相关性与初始相关性矩阵几乎完全一致，这表明采用本发明获得的样本相比完整地保存了输入参数的相关性，即采用本技术传播的不确定性更加合理。

另一方面，为了确保采用本技术的“小样本”方案能使不确定性信息完整地传播，需要同时对输入参数的不确定性进行验证，验证结果如图4所示。图4分别展示了44能群核截面初始不确定性信息(基准)、采用对比方案LHS方法传播的不确定性信息、采用本发明即改进方案传播的不确定性信息。可以看到，采用本发明得到的描述不确定性的曲线与初始不确定性曲线几乎完全重合，因此可以判断本发明对于不确定性的传播效果明显好于作为对比方案的LHS方案，进一步证明了本发明在高维输入参数下能保证不确定性的合理传播。

Claims (11)

1.一种针对反应堆物理分析中高维核截面输入参数的不确定性传播方法，其特征在于，包括步骤如下：

步骤1：根据核反应堆物理计算不确定性分析需求输入核截面参数个数，读取输入核截面参数及相应维度的核截面相对协方差矩阵；

步骤2：对读取的核截面相对协方差矩阵进行主成分分析，确定降维后的特征维度，根据特征维度建立准确表征原始不确定性信息的核截面相对协方差矩阵；

步骤3：在选取的核截面总体分布空间上进行等区间划分，确定每个区间上随机样本的累计概率密度；

步骤4：通过求解正态分布的反函数，获取相应累积概率密度对应的随机样本；

步骤5：建立随机样本的协方差矩阵；

步骤6：对随机样本空间协方差矩阵进行奇异值分解，获得左奇异矩阵U，并构造对角矩阵D；

步骤7：利用左奇异矩阵U及对角矩阵D，构建新的随机样本空间；

步骤8：利用线性关系获得表征原始核截面不确定性信息的核截面最终样本数据后，根据核截面样本数据进行相关性验证及不确定性验证，准确传播高维核截面的不确定性。

2.如权利要求1所述的不确定性传播方法，其特征在于，所述步骤1具体如下：

对于一个样本空间Z_nS×nX，其协方差矩阵为∑_nX×nX，根据奇异值分解理论得到协方差矩阵可以表示为：

其中，nS为样本数量，nX为参数个数，矩阵S_nX×nX为奇异值矩阵，矩阵U_nX×nX为左奇异矩阵。

3.如权利要求2所述的不确定性传播方法，其特征在于，所述步骤2具体如下：

对奇异值进行分析处理，选取能表征矩阵∑_nX×nX大部分信息的t个奇异值构造新的奇异值矩阵S_1t×t来表示协方差矩阵∑_t×t，如下式所示：

∑_t×t＝U₁S₁U₁ ^T (2)

此时即完成了对协方差矩阵为∑_nX×nX的降维，其中，t<nX。

4.如权利要求3所述的不确定性传播方法，其特征在于，所述步骤3具体为：在总体分布空间上进行等区间划分，在该区间内随机抽取(0,1)的随机数，确定每个区间上随机样本的累计概率密度，方法为：

p_i＝(i-1)/nS+r_u/nS (4)。

5.如权利要求4所述的不确定性传播方法，其特征在于，所述步骤4具体为：通过求解正态分布的反函数，计算相应累积概率密度对应的随机样本，方法为：

为第j个输入参数分布函数的反函数。

6.如权利要求5所述的不确定性传播方法，其特征在于，

步骤5具体为：计算随机样本协方差矩阵，其中协方差计算方法为：

Cov(x_i,x_j)＝E(x_ix_j)-E(x_i)E(x_j) (6)

其中x_i，x_j为随机样本序列，将计算得到协方差按照对应位置排列即可得到协方差矩阵。

7.如权利要求6所述的不确定性传播方法，其特征在于，所述步骤6具体方法为：基于步骤4得到的样本空间Z_s＝[Z_s1,Z_s2,....,Z_snX]^T，该样本空间的协方差矩阵∑_s的格式是：

参数Z_si和Z_sj由抽样样本得到的协方差估计值如式(8)所示：

将参数间的协方差表达式(8)代入到协方差矩阵∑_s中，得到协方差矩阵∑_s格式如下：

其中，nS为样本数量，H是一个n×n满矩阵，其元素均为1；

对协方差矩阵∑_s进行奇异值分解，如式(10)：

∑_s＝USV^T (10)

由于协方差矩阵∑_s为对称正定矩阵，因此有U＝V，则式(10)可转换为如下形式：

U^T∑_sU＝S (11)

其中，S为对角阵且其元素为矩阵∑_s的奇异值，构造对角矩阵E和对角矩阵D，其格式分别如式(12)-(14)所示：

其中，s₁,s₂,....,s_nX为矩阵∑_s的奇异值，并且按照从大到小的顺序排列；

8.如权利要求7所述的不确定性传播方法，其特征在于，所述步骤7具体方法为：

利用左奇异矩阵U及对角矩阵D，构建新的随机样本空间；

进而可得：

(D^-1)^TSD^-1＝E (15)

结合式(11)可得：

(D^-1)^TU^T∑_sUD^-1＝E (16)

将式(9)带入式(16)可得：

对于样本空间

其协方差矩阵

为单位对角矩阵，表明样本空间

中的各参数间完全相互独立，

即为构建的新样本空间。

9.如权利要求5所述的不确定性传播方法，其特征在于，所述步骤4中，还包括通过对随机样本的相关系数矩阵进行楚列斯基分解，并进一步构建新的随机样本，具体的其方法为：

对步骤4得到的随机样本空间Z，对其相关系数矩阵C进行楚列斯基分解如下：

C＝QQ^T (18)

对上三角矩阵Q求得逆矩阵Q^-1，利用逆矩阵Q^-1对随机样本空间Z进行变换，如式(19)所示：

10.如权利要求8所述的不确定性传播方法，其特征在于，所述步骤8的具体方法为：

其中X为最终样本数据，μ为输入参数本征值，A为对输入参数协方差矩阵采用奇异值分解获得的左奇异矩阵。

11.如权利要求1至10任一项所述的不确定性传播方法，其特征在于，所述高维核截面输入参数的不确定性传播是44×44的核截面不确定性传播。

Images