基于氪和氙的同位素的核爆炸识别方法

摘要

本发明涉及核物理领域，可用于基于在自然发生的放射性气体（NORG）中测得的核活性来识别核爆炸的系统。其技术效果是提高了确定效率和从对于每种氪和氙的同位素总活性中的各种裂变类型进行沉积物的精确预估的可靠性。

说明书

发明领域

本发明涉及核物理，可用于确定大气中的放射性物质的系统中。

背景技术

在为了监视全面禁止核试验条约（CTBT）的遵守情况而监测放射性状况 的过程中，要根据放射性惰性气体（RNG）的同位素进行核爆炸的识别。

一种远距离探测核材料的方法已众所周知，该方法包括通过使用电子束 照射激发被拦截物体，然后记录从该物体发射的中子数（美国专利第4320298 号,1982）。

上述方法的缺点在于其不适用于大气条件的检测和对不包含氘元素的对 象的检测。

多种远距离探测核材料的方法已为人所知，这些方法包括通过使用中子 束激发被拦截物体，然后记录从该物体发射的中子数或伽玛-射线（SU  Inventor's Certificate  #439740（SU发明人证书#439740）;SU Inventor's  Certificate#1349478（SU发明人证书#1349478）;美国专利第4483817号）。

这些方法的缺点在于它们不适用于由于某些或其他原因而禁止外部辐射 作用于拦截的对象上的情况。

一种远距离探测核材料的方法已众所周知，该方法包括在拦截物体周围 检测能量范围在从0.1到2.0MeV内的γ射线的强度（Sagdeyev,R.Z.et al. Problems of Monitoring sea-based cruise missiles with nuclear warheads. Preprinted by IKI of the USSR AS,Pr-1373.-M.,1988（Sagdeyev,R.Z等人所著的 “在监测带有核弹头的海基巡航导弹中的问题”，USSR AS的IKI预印， Pr-1373.-M.，1988年））。

这种方法的一个缺点是存在误检测核材料的可能性，因为带有这种能量 的辐射也可能是含有放射性物质的非爆炸性对象产生的。

通过放射性氪和氙同位素确定核爆炸的方法已众所周知。

已知有一种通过放射性氪和氙同位素确定核爆炸的方法（SU Inventor's  Certificate#366771（SU发明人证书#366771））。在这种方法中，测量RNG 活性的方式如下：将样品带入大气（在一项事件后）中，进行一段时间的研 究。在此期间，利用单晶闪烁伽马谱仪NaJ（Tl）测定氪和氙的同位素活性。

因此，所测得的氪和氙的同位素的活性被用于形成一种与RNG源对氪和 氙的同位素的总活性的未知贡献相关的线性代数方程组（SLAE）。

该方程组通过使用最小二乘法（LSM）来求解。

在技术本质上最接近的是一种可以消除上面方法的一些缺点的方法，即 考虑待解方程组中的矩阵元素的错误，并通过使用A.N.Tikhonov的正则化方 法确保获得一个稳定的解（Greshilov A.A.,Tetjukhin A.A."An Algorithm for  Identifying Sources of Radioactive Noble Gases".Bulletin of MGTU Named After  N.E Bauman.Natural Sciences Series,2003.No.2,p.3-19.（“用于识别放射性惰 性气体源的算法”，Greshilov A.A.和Tetjukhin A.A，自然科学系列，在N.E. Bauman后以MGTU命名的公报，2003年，第2期，第3-19页））

已知的通过放射性氪和氙同位素识别核爆炸的方法包括：

1.测量分离的放射性氪和氙同位素在大气中的活性（t为测量时 间），其中，n是测量同位素的数量。

2.通过氪和氙同位素的活性关系为不同类型的可裂变材料（一种裂变类型 可理解为采用裂变光谱的中子或者带有14MeV的能量的中子产生的U²³⁵,U²³⁸, PU²³⁹的重核的裂变体之一）确定分离发生的时间间隔[t_H，t_K]，它们是在充分 考虑和未充分考虑分离的情况下建立的；

3.在时间间隔[t_H，t_K]上设置间隔为Δt的时间网格；

4.为每个网格节点tq形成和记忆氪和氙同位素的比活度的二维信号 {a_ij(t_q，t)}，这些节点所在的行对应特定的同位素，列对应特 定的裂变类型；N为被研究的同位素的数量；m为被研究的裂变类型的数 量；t为样品测量时间；tq为假定的分离时间；

5．将测量的活性值和二维信号的成分{a_ij(t_q，t)}分别作为统计上不相 关的量，这些量按照带有均方根偏差σ(A_i(t))和σ(a_ij(t_q，t))的正态法则（normal  law）分布，在设定的时间tq形成规则的信号：

$F_l=\frac{1}{2}\underset{i=1}{\overset{n}{\Sigma }}(\frac{{({\tilde{A}}_i\left(t\right)-\underset{j=1}{\overset{m}{\Sigma }}{a}_{\mathrm{ij}}^{\mathrm{true}}(t_q,t)·\left(\rho N_j\right))}^2}{{\sigma }^2\left({\tilde{A}}_i\left(t\right)\right)}+\underset{j=1}{\overset{m}{\Sigma }}\frac{{(a_{\mathrm{ij}}(t_q,t)-a_{\mathrm{ij}}^{\mathrm{true}}(t_q,t))}^2}{{\sigma }^2\left(a_{\mathrm{ij}}(t_q,t)\right)}),$

其中：

(ρN_j)，为放射源对活性总量所做的贡献；

为比活度的未知真值；

a_ij(t_q，t)为通过放射性同质异位核转化链（RTIC）的独立的和累积的成分 计算的特定活性，其中有误差；

l为用于和估计的搜索的迭代数。

6.设置γ1,γ2数量，它们表示和估计的精度特征。

7.通过使用A.N.Tikhonov的正则化，反复查找和中的一 个最小的信号F1，使其符合 $|a_{\mathrm{ij}}(t_q,t)-{\hat{a}}_{\mathrm{ij}}^{\mathrm{true}}(t_q,t)|\le 3\sigma \left(a_{\mathrm{ij}}(t_q,t)\right),i=\mathrm{1,2},...,n,j=\mathrm{1,2},...,m,$的限制，直到满 足以下条件：

$\left|\right|\frac{{\left(\widehat{{\mathrm{\rho N}}_j}\right)}_l-{\left(\widehat{{\mathrm{\rho N}}_j}\right)}_{l-1}}{{\left({\widehat{\mathrm{\rho N}}}_j\right)}_l}\left|\right|<{\gamma }_1,$

$\left|\frac{F_l-F_{l-1}}{F_l}\right|<{\gamma }_2,$

其中指数（l-1）是在前面的迭代中获得的值。

8.通过关系式$D_{\mathrm{ij}}={(-\frac{{\partial }^2{F}_l}{\partial \left(\widehat{{\mathrm{\rho N}}_i}\right)\partial \left(\widehat{{\mathrm{\rho N}}_j}\right)})}^{-1},i,j=\mathrm{1,2},...,m.$对估计协方差矩阵进行确 定。

但是，这个已知的方法存在以下缺点：1）只使用了Tikhonov正则法，它 需要另外确定正则化参数，而尚没有发现找到它的独特方法；Tikhonov的正 则法使解“平滑（smooth）”，这可能在识别中造成很大的错误；

2）所考虑的方程组是过度定义的（overdefined），即一个两维的信号 {a_ij(t_ij，t)}的行数大于它的列数—方程的数目大于未知数的个数。当氪和氙的同 位素的放射性在事件后几天内被测量，且未知数ρN_j（各种来源的贡献）的数 量大于被测量的同位素数目即量的数目时，在这样的实际情况中运用这 种方法的可能性尚未得到研究。

因此，当测得2-4种同位素的活性，并且是在相应的事件后5-6天内取样 的，已知的解决方法并不能保证在最有可能的情况下识别核爆炸。

显然，由于实际应用的低效性，上述缺陷使得这种方法难以付诸实际应 用。

发明内容

所要求保护的方法的技术效果是对核爆炸事实的确定可靠性的提高，所 测量的同位素的数目小于所考虑的未知数（裂变类型）的数目。

所要求保护的方法的功效通过以下方式得以确保实现：

1）同时考虑氪和氙同位素活性源和核爆炸的各种机制的各种假想组合；

2）开发多准则的数学规划（programming）机制并将其纳入核爆炸的识 别方法中，以便能够考虑一个解的估计应该符合的所有可能类型的附加条件 （解的非负性，解的边界性）、正则化参数（由A.N.Tiknonov提出的）以及 没有定义的清楚地格式化的程序；

3）通过将RTIC元素的独立和累积输出平均，由裂变频谱中子和在单裂变 类型中带有14MeV能量的中子进行两种类型的U²³⁵裂变的组合，以及由裂变 频谱中子和单裂变类型中的带有14MeV能量的中子进行两种类型的Pu²³⁹裂变 的组合。

要求保护的发明的技术效果是通过开发一种通过放射性氪和氙同位素识 别核爆炸的方法实现的，其特征是：在事件后的t时刻对描述在测量站附近的 大气中每种同位素的总活性的变化的信号Ai(t)进行测量；建立同位素活性和 时间之间的关系而无需适当考虑同位素活性之间的分离与关系，这些同位素 活性是为所有考虑的裂变类型逆时序（reverse time）从测量点提取的；确定 分离间隔[t_H，t_K]；在该时间间隔[t_H，t_K]内设置时间网格；形成裂变类型的组 合；对于该分离间隔内的每一个网格节点tq和每一个组合，根据分离时间tq 和测量时间t计算描述每一种同位素的“具体”活性值的二维信号A和（某种可裂 变材料和中子能量的）潜在源；设定测量信号Ai(t)误差的均方值σ(Ai(t))和二 维信号成分{a_i；j(t_q,t)}的误差的均方根值σ(a_ij(t_q,t))；设定定义了信号估计pNj 和J2信号计算的精度的量γ1,γ2；

通过利用信号J1的形成和信号J2识别核爆炸，信号J1确定从信号Ai(t)与 二维信号{a_i；j(t_q,t)}的行乘上信号pNj的值所得的积之间的差的平方和所得到的 解的精度，信号J2确定了pNj信号的类型；从所设置的组合中的信号J1和J2 形成限制信号和目标信号（目标函数）；以及通过使用迭代程序并改进（refine） 每次迭代，寻找pNj信号中成分的二维信号{a_i；j(t_q,t)}的值；根据确定每种裂变 方式可能对总活性作出的活性贡献的点估计，通过差异平方的和（sum of  discrepancy squares）选择最佳的裂变类型组合。

该方法的特征还在于将来自不同能量组的中子的U²³⁵和Pu²³⁹的两种裂 变类型组合成单一裂变形式和这是通过与不同裂变类型相应的同 重元素链中的元素（isobaric chain elements）的独立输出以及累积输出的加权 求和来实现的；根据权重c₁和c₂建立二维网格，为二维信号元素{a_i；j(t_q,t)}cc 的每对权重值c₁,c₂计算比活度；根据第1项计算信号估计(pN_j)_cc；选择一个估 计(pN_j)_cc使

$\underset{i=1}{\overset{n}{\Sigma }}{({\tilde{A}}_i\left(t\right)-\underset{j=1}{\overset{m}{\Sigma }}{\left(a_{\mathrm{ij}}^{\mathrm{true}}(t_q,t)\right)}_{c_1{c}_2}{\left(\widehat{{\mathrm{\rho N}}_j}\right)}_{c_1{c}_2})}^2$

的值最小。

Kr_83m,Kr_85m,Kr₈₅,Kr₈₈,Xe_131m,Xe_133m,Xe₁₃₃,Xe₁₃₅作为单独的氪和氙同位 素。

为了通过少量的测量同位素来确立(establish)核爆炸的事实，通过对RTIC 元素的对应于上述裂变类型的独立输出和累积输出取平均，将U²³⁵裂变两种类 型（通过裂变光谱的中子和带有14MeV的能量的中子）组合成单一裂变类型， 将Pu²³⁹裂变的两种类型（通过裂变光谱的中子或者带有14MeV的能量的中子） 组合成单一裂变形式。

来自个别源的活性贡献通过形成几个目标函数和使用多准则数学规 划的方法来估计，把一个多准则数学规划任务减少到有限制的单准则任务， 通过迭代计算过程得到对上述有限制的单准则任务的解，根据值（j型裂 变的数量）确定各个源对总活性的贡献，也就是，识别一个核爆炸的参数。 一个赋予了Fl信号一个最小值的RNG同位素来源的组合在求解处理中被认为 是真值。

通过氪和氙的放射性同位素识别核爆炸的方法的最佳实施方式。

所要求保护的发明的实质将在下文中通过其文字说明和图形资料来解 释，其中：

图1显示了在一个没有分离的U^f₂₃₅和Pu^f₂₃₉（分别是实线和虚线）的裂变 的情况中同位素A(Xe133m)/A(Xe133)的活性的相对变化，并适当考虑与前面 的同位素分离的情况（有标记的线），点1是在时间t=12小时时测得的同位素 活性的关系。

图2显示了一个用于获得分离时间和解的估计的算法的总体框图。

图3显示了通过使用融合分析（confluent analysis）找到条件不佳的线性代 数方程组的解的框图（相当于图2所示的总体框图的方框1）。

图4显示了同位素Xe¹³³和Xe¹³⁵的累积输出与裂变光谱的中子或者带有 14MeV的能量的中子的相关贡献的依赖关系的图形。

为了执行这个方法来通过放射性氪和氙的同位素识别核爆炸，需要根据 大气中RNG来确定核爆炸来源的下列参数：

每种裂变类型（可裂变的材料和中子能量对大气中测量的总 RNG活性的贡献。

作为本要求保护的方法的一个实际上合理的假设，（同位素活性的）注 册信号被认为具有确定性(deterministic)，受到额外的干扰，其参数估计应该被 定义。

在裂变瞬间，i-次同位素表现为不同类型的裂变的结果，其测得的活性 可表示如下[3]：

$\underset{j=1}{\overset{m}{\Sigma }}{a}_{\mathrm{ij}}(\theta ,\eta ,\lambda ,t,t_q)\rho N_j={\tilde{A}}_i\left(t\right),i=\mathrm{1,2},...,n,---\left(1\right)$

其中：a_ij(θ，η，λ，t，t_q)为i型同位素在j型裂变中一次衰变的活性，这是通过 在时间t>tq的分离计算得出的，即比活度；

θ为参数矢量，表征测量到的同位素与在前的(preceding)同位素的分离；

η为同位素（j型裂变）的独立输出的矢量；

λ为衰变常数矢量；

t为观测时间；

tq为在放射衰变链（chains of radioactive transformations）中氪和氙同位素 与在前的同位素分离的假定分离时间；

ρ是i型同位素在样品中的的比率（ρ值通常是未知的）；

Nj为j型的裂变的次数。

直到分离时间tq，比活度由以下公式决定：

$a_{\mathrm{nj}}^{l_q}(\theta ,\eta ,\lambda ,t_q)=\{{\eta }_n{\lambda }_n\exp (-{\lambda }_n{t}_q)+\underset{p=1}{\overset{p_{\max }}{\Sigma }}\underset{i_p=1}{\overset{n_p-1}{\Sigma }}{\eta }_{i_p}{\lambda }_n[\underset{r_p=i_p}{\overset{n_p-1}{\Pi }}{\gamma }_{r_p}{\lambda }_{r_p}\underset{s_p=i_p}{\overset{n_p}{\Sigma }}\frac{\exp (-{\lambda }_{s_p}{t}_q)}{\underset{q_p\ne s_p}{\overset{n_p}{\underset{q_p=i_p}{\Pi }}}({\lambda }_{q_p}-{\lambda }_{s_p})}\left]\right\},---\left(2\right)$

其中：

ηi是i型同位素的独立输出；

n_p是研究同位素沿p-分支的数量；

n是{n_p}中的最大值；

p_max是支链的数目；

(n_p–1)是沿衰变的p-分支中在所研究的同位素之前的同位素的数量；

是链中r成员的比率，这是沿p-链的（r-1）处获得的；

是同位素的衰变常数，这些同位素沿p-分支分别有数目 i_p，r_p，s_p，q_p，其中的i_p≤r_p≤n_p-1；i_p≤s_p≤n_p；i_p≤q_p≤n_p且q_p≠s_p；

tq是所研究的同位素与在先的同位素发生瞬间分离的时间，该时刻后同位 素按衰变常数λn进行指数衰变。

在按照衰变常数λi的同位素衰变的分离时间后，λi： $a_{\mathrm{ij}}(\theta ,\eta ,\lambda ,t,t_q)=a_{\mathrm{ij}}^{t_q}(\theta ,\eta ,\lambda ,t_q)\exp (-{\lambda }_i(t-t_q)),i=\mathrm{1,2},...,n,j=\mathrm{1,2},...,m,---\left(3\right)$其中， 是根据公式（2）在tq时刻计算得到的比活度。

类型（1）方程式用公式表示了每个测定的氪和氙的同位素，在结果中， 形成了SLAE

其中，定义了放射源ρNj对氪和氙同位素的总活性的未知贡献。

解决识别RNG任务的第一步是定义氪和氙的同位素分离时间tq。

通过完成“构建(building)”从测量时间起“逆时间”中不同裂变类型中相关 的同位素活性而不考虑在前的同位素的影响，以及通过确定从多个实验点绘 出的多条线的交叉点和适当考虑沿同位素衰变链上的在前同位素的影响而建 立的相关活性，可以寻找分离时间所属的时间间隔区。

图1显示了两个氙同位素（Xe133m,Xe135）相对活性的图线，其中实验点 1对应在事件后t=12小时测量活性的时间。为了简化图，只显示了“边界”线， 对应U^f₂₃₅和Pu^f₂₃₉（而不是6个可能的裂变类型：和

如图1所示，分离时间属于事件后t_H=3到t_K=4小时的时间间隔内。

当在时间间隔[t_H,t_K]内设置一个时间网格，和求解矩阵式（system）（4） 以得到对应于网格节点的t_q时，时间因为它的矩阵式（4）的差的平方和 （sum of squares of discrepancies）是最小的，作为分离时间。

图2显示了可以找到估计的总体算法的框图。

解决识别核爆炸任务的第二个步骤是，为每一个固定分离时间tq确定解估 计 

在一个给定的tq，矩阵式（4）与未知的呈线性相关。

由于二维信号A＝{a_ij(t_q，t)}的成分可能无法准确计算（已知独立输出有误 差）且同位素的活性的测量也有误差，我们认为二维信号A的元素和测量活性 是按照正态分布（normal law）的独立的随机数量，数学期望等于 和方差等于${\sigma }^2\left(a_{\mathrm{ij}}(t_q,t)\right)\mathrm{and}{\sigma }^2\left({\tilde{A}}_i\left(t\right)\right),$分别为：

$a_{\mathrm{ij}}(t_q,t)=a_{\mathrm{ij}}^{\mathrm{true}}(t_q,t)\pm {ϵ}_{\mathrm{ij}},$(5)

${\tilde{A}}_i\left(t\right)={\tilde{A}}_i^{\mathrm{true}}\left(t\right)\pm {\delta }_i,$

其中：

是同位素的具体和测量活性的真值（这是我们未知的）；

ε_ij是在确定比活度a_ij(t_q，t)时的误差；

δ_i是测量大气中RNG的活性的误差。

为了考虑在测得的活性和二维信号{a_ij(t_q，t)}的成分中的误差，使用了 正交回归[4]的定义，以及由于随机变量a_ij(t_q，t)和信号的独立性，信号可 写为：

$F_l=\frac{1}{2}\underset{i=1}{\overset{n}{\Sigma }}(\frac{{({\tilde{A}}_i\left(t\right)-\underset{j=1}{\overset{m}{\Sigma }}{a}_{\mathrm{ij}}^{\mathrm{true}}(t_q,t)·\left(\rho N_j\right))}^2}{{\sigma }^2\left({\tilde{A}}_i\left(t\right)\right)}+\underset{j=1}{\overset{m}{\Sigma }}\frac{{(a_{\mathrm{ij}}(t_q,t)-a_{\mathrm{ij}}^{\mathrm{true}}(t_q,t))}^2}{{\sigma }^2\left(a_{\mathrm{ij}}(t_q,t)\right)}),---\left(6\right)$

其中：

(ρN_j)，放射源对总活性的贡献；

是未知的比活度的准确值，在寻找ρN_j的过程中它的估计变得更 精确，a_ij(t_q，t)是按照由独立和累积输出根据公式（2）-（3）计算得出的放射 性衰变的同重元素链的元素的比活度，其具有误差；

是样本中测量的RNG。

除了一个未知放射性源贡献矢量ρN_j，j＝1，2，...，m，信号（6）还包括计 算的活性的未知真值f＝1，2，...，n，j＝1，2，...，m,通过使用融合性分析 [3，4]找出该信号的估计值。

信号（6）的最低点应符合以下条件符合：

$\frac{\partial F_l}{\partial \left(\rho N_j\right)}{|}_{\rho N_j=\widehat{{\mathrm{\rho N}}_j}}=0,j=\mathrm{1,2},...,m;---\left(7\right)$

$\frac{\partial F_l}{\partial a_{\mathrm{ij}}^{\mathrm{true}}(t_q,t)}{|}_{a_{\mathrm{ij}}^{\mathrm{true}}(t_q,t)={\hat{a}}_{\mathrm{ij}}^{\mathrm{true}}(t_q,t)}=0,i=\mathrm{1,2},...,n;j=\mathrm{1,2},...,m---\left(8\right)$

最小信号（6）的一个结构图如图3所示，它相当于在图2所示的结构图的 方框1。

尽管方程组（7）-（8）在固定点tq是线性的，由于方程组（7）缺少条件， 从计算的角度来说该任务是不正确的。在方程组（7）的矩阵的最大和最小的 适当的数之间的关系可以是1026数量级（order of 1026）。因此，为了求解， 需要采用特定的技术，在要求专利保护的方法中是多准则的数学规划 （multi-criterion mathematical programming），该方法不需要像其他解决不正 确的任务的方法那样确定正则化参数值。

作为第一步，在SLAE(7)通过多准则的数学规划（它 是一种压缩允许值的区域的方法，即目标规划）求解，估计的第一接近 值被找到了。为了获得在给定时刻tq下的真值的估计值，在获得估 计值j＝1，2，...，m的每一个步骤中使用了条件（8）[3]，其结果是能够另外 求解出以下形式的带有未知数m的n个线性方程组：

$\left(\begin{array}{c}\underset{r=1}{\overset{m}{\Sigma }}\frac{\left(\widehat{{\mathrm{\rho N}}_r}\right)\left(\widehat{{\mathrm{\rho N}}_v}\right)}{{\sigma }^2\left({\tilde{A}}_i\left(t\right)\right)}{\hat{a}}_{\mathrm{ir}}^{\mathrm{true}}(t_q,t)+\frac{{\hat{a}}_{\mathrm{iv}}^{\mathrm{true}}(t_q,t)}{{\sigma }^2\left(a_{\mathrm{iv}}\right)}=\frac{a_{\mathrm{iv}}(t_q,t)}{{\sigma }^2\left(a_{\mathrm{iv}}(t_q,t)\right)}+\frac{\left(\widehat{\rho N_v}\right){\tilde{A}}_i\left(t\right)}{{\sigma }^2\left({\tilde{A}}_i\left(t\right)\right)},\\ i=\mathrm{1,2},...,n,v=\mathrm{1,2},...,m.\end{array}\right)$

由此获得的值估计应符合自然条件，即属于测量值a_ij(t_q，t)的不确定范围 内。

$|a_{\mathrm{ij}}(t_q,t)-{\hat{a}}_{\mathrm{ij}}^{\mathrm{true}}(t_q,t)|\le 3\sigma \left(a_{\mathrm{ij}}(t_q,t)\right),i=\mathrm{1,2},...,n,j=\mathrm{1,2},...,m.$

如果不符合这个条件，那么，这等式不符合这种不 等式，应该由最近的边界点的值替换。由于这个原因，与那些在迭代过程中 先前的数值相比，信号Fl的值可能会增加到新的常数值，这将导致在降低迭代 过程的收敛速度，或出现波动。

为了在再次计算估计${\hat{a}}_{\mathrm{ij}}^{\mathrm{true}}(t_q,t),j=\mathrm{1,2},...,m$后，函数值不增加，估计值 集合（Fl的相应分量与他们在先前迭代中的值相比增加 了）应该被前面步骤的相应值代替。

在调整了真实值的估计后，通过多准则的数学规划的方法找出 ρN_j，j＝1，2，...，m的下一个近似解而不是像在类似的解法中那样使用 A.N.Tikhonov提出的正则方法。

停止该算法的标准是信号Fl的值和相邻迭代的向量内的ρNj矢量的成分之 间差异很小，即满足以下不等式：

$\left|\right|\frac{{\left(\widehat{{\mathrm{\rho N}}_j}\right)}_l-{\left(\widehat{{\mathrm{\rho N}}_j}\right)}_{l-1}}{{\left(\widehat{{\mathrm{\rho N}}_j}\right)}_l}\left|\right|<{\gamma }_1,---\left(9\right)$

$\left|\frac{F_l-F_{l-1}}{F_l}\right|<{\gamma }_2,$

其中：

是下一个接近于第一次迭代的解；

γ₁，γ₂是一些数字（十进制的小数部分，如0.001），定义了的估计的 计算值的精度。

当通过多准则的数学规划方法求解时：

形成一个数学规划的二准则的任务：

$J_1=\underset{i=1}{\overset{n}{\Sigma }}{({\tilde{A}}_i\left(t\right)-\underset{j=1}{\overset{m}{\Sigma }}{a}_{\mathrm{ij}}^{\mathrm{true}}(t_q,t)\left(\rho N_j\right))}^2\to \underset{\rho N_j}{\min },$(10)

$J_2=\underset{j=1}{\overset{m}{\Sigma }}\rho N_j\to \underset{\rho N_j}{\min },$

限制条件为$\rho N_j\ge 0,j=\overline{1,m}.$

在这里，J1是一个调节类型（4）方程组的差异平方和的信号，以确保为 氪和氙同位素测量活性匹配一个估计解J2是形成一种解类型的信 号。

2）利用阈值优化方法或目标规划，通过将除了一以外的上述所有的函数 转换（transfer）至极限，该算法从数学规划（10）的二准则任务转变为单准 则任务。

阈值的优化方法（或e-极限（e-limitation）的方法）产生目标函数和限制 条件的各种可能的组合。该算法适用于以下类型：

$\underset{\rho N_j}{\min }\underset{i=1}{\overset{n}{\Sigma }}{({\tilde{A}}_i\left(t\right)-\underset{j=1}{\overset{m}{\Sigma }}{a}_{\mathrm{ij}}^{\mathrm{true}}(t_q,t)·\left(\rho N_j\right))}^2$当$\underset{j=1}{\overset{m}{\Sigma }}\rho N_j\le \delta ;\rho N_j\ge 0,j=\overline{1,m};---\left(11\right)$

$\underset{\rho N_j}{\min }\underset{j=1}{\overset{m}{\Sigma }}\rho N_j$当$\underset{i=1}{\overset{n}{\Sigma }}{({\tilde{A}}_i\left(t\right)-\underset{j=1}{\overset{m}{\Sigma }}{a}_{\mathrm{ij}}^{\mathrm{true}}(t_q,t)·\left(\rho N_j\right))}^2\le \beta ;\rho N_j\ge 0,j=\overline{1,m}.---\left(12\right)$

任务（11）是用于二次规划，任务（12）是用于非线性规划。

右侧的限制条件δ和β的估计可以通过带有极限的函数J1 和J2的独立的最小化中获得。在这里，可以使用任何数学规划方法。

目标规划包括两种解模型---阿基米德模型和有优先顺序的模型（the model with priorities）。

如果使用阿基米德模型，所有的目标函数被转换至极限内，它们与极限 的偏差的加权总和被最小化：

$\underset{d_1,d_2}{\min }\{-(w_1{d}_1+w_2{d}_2)\}$于J₁(ρN_j)+d₁≤β，J₂(ρN_j)+d₂≤δ，       (13)

其中：

w_i是加权系数，d_i是与极限的偏差。

如果使用有优先顺序的模型，目标函数继续被转换至极限，目标函数值 与极限的偏差被最小化。这里，在这一步被发现的偏差值di在第i+1步用做一 个最佳偏差：

第1步：$\underset{d_1}{\min }\left({-d}_1\right)$于J₁(ρN_j)+d₁≤β；

第2步：$\underset{d_2}{\min }(-d_2)$于$J_1\left(\rho N_j\right)+d_{\mathrm{lonm}}{|}_{d_{\mathrm{lonm}}=d_1}\le \beta ,J_2\left(\rho N_j\right)+d_2\le \delta .---\left(14\right)$

为了通过使用少量氙同位素识别核爆炸，利用了这些同位素的输出不充 分地依赖于中子能量的事实。

图4中说明了该情况，该图示出了Xe133和Xe135同位素的平均累积的输 出的绘图(plot)，其取决于与裂变频谱中子相应的累积输出的部分和与带有14 MeV的能量的中子相应的累积输出的部分的比例。已知累积输出的误差可达 到5％。

从图4中可以看到，平均输出值（作为例子，示出了对应于相同比例的裂 变光谱中子和带有14MeV的能量的中子的值）在这个误差内。

两种U²³⁵裂变类型（由裂变光谱的中子和带有14MeV的能量的中子产生） 组合成一种裂变类型，而两种Pu²³⁹裂变类型（由裂变光谱的中子和带有14MeV 的能量的中子产生）组合成一种裂变类型，使用这样的组合来通过2-4个测量 的同位素进行识别，其结果是减少了可识别的裂变类型的数量（2种代替了4 种）。

在这里，比活度的计算公式如下：

1）对于可裂变的U²³⁵-根据公式(2)并利用独立的输出向量 的，其中和是当U²³⁵分别被裂变光谱的中子和带 有14MeV能量的中子裂变后同重元素链中元素的独立输出，c₁是考虑到独立 输出和占其总和c₁∈[0，1]的比例的参数。

2）对于可裂变的Pu²³⁹，根据公式（2）并利用独立的输出向量 其中，和是当Pu²³⁹分别被裂变光谱的中子和 带有14MeV能量的中子裂变后同重元素链中元素的独立输出，c₁是用于解释 独立输出和占他们的和c₂∈[0，1]的比例的参数。

通过c₁和c₂设置一个分别有间距Δc₁和Δc₂的二维网格，并为不同的c₁和c₂找到信号的最小值（6），在公式（2）的差异平方和是最小的时候来源 和为氪和氙同位素的总活性所贡献的量被认为是真实的。

因此，所建议的识别核爆炸参数的方法，可以如下进行：

1、测量大气中的单独的氪和氙的同位素（Kr_83m,Kr_85m,Kr₈₅, Kr₈₈,Xe_131m,Xe_133m,Xe₁₃₃,/Xe₁₃₅）的活性（t为测量时间），其 中，n是测得的同位素的数目。

2、建立同位素活性关系对时间的依赖关系，但是不考虑从事 件（和爆炸）出现的时间到测得氪和氙的同位素（例如Kr_85m/Xe₁₃₅） 活性的时间内不同类型裂变材料的分离。

3、在不考虑分离的情况下，在“逆时序”内确定同位素（例 如，Kr_85m/Xe₁₃₅）的关系，该“逆时序”是从测得的实验点开始直到 它们匹配相同同位素的相对活性值。

4、根据匹配值确定一段时间间隔[t_H，t_K]，在该时间间隔中氪和 氙的同位素与其在放射性衰变的同重元素链上先前的同位素分离。

5、在间隔[t_H，t_K]内设置具有间距Δt的时间网格。

6、对网格的每个节点，形成并记录氪和氙同位素的比活度的 二维信号

其中，行对应特定的同位素，列对应特定的裂变类型；m是所 考虑的裂变类型的数目。这里的a_ij(t)对应于j型裂变的i型同位素的比 活度，根据公式（2），计算直到分离时间tq的该活性，根据公式 （3），计算在分离时间tq后的t时刻的

7、考虑到已测量的活动值和二维信号{a_ij(t_q，t)}的成分在是统计学上 的独立事件，均符合数学期望为和均方差为的正态分布，根据公式(6)形成一维信号F_l。

8、设置表征估计和的精度的数γ₁,γ₂。

9、在首次迭代为数学规划形成一个二准则的任务

$J_1=\underset{i=1}{\overset{n}{\Sigma }}{({\tilde{A}}_i\left(t\right)-\underset{j=1}{\overset{m}{\Sigma }}{a}_{\mathrm{ij}}^{\mathrm{true}}\left(t_{q,}t\right)·\left(\rho N_j\right))}^2\to \underset{\rho N_j}{\min },$

$J_2=\underset{j=1}{\overset{m}{\Sigma }}\rho N_j\to \underset{\rho N_j}{\min },$

取$a_{\mathrm{ij}}^{\mathrm{true}}(t_q,t)=a_{\mathrm{ij}}(t_q,t)$（随后的迭代中取$a_{\mathrm{ij}}^{\mathrm{true}}(t_q,t)={\hat{a}}_{\mathrm{ij}}^{\mathrm{true}}(t_q,t)),$限制（limitation） 是$\rho N_j\ge 0,j=\overline{1,m}.$

10、使用压缩允许的值域或目标规划的方法，该算法通过将除一之外的 所有函数上述所有的函数转换到极限条件，从符合公式(11)-(14)的数学规划二 准则的任务转变为单准则的任务。

11、使用二次规划、非线性规划和目标规划（阿基米德模型和有优先顺 序的模型），求取单独来源的活性对总活性的贡献评估的第一个近似值

12、在获得第一个近似值后，二维信号的成分{a_ij(t_q，t)}被更精确地规 定。为此，求解n个方程m个未知数的线性方程组：

$\left(\begin{array}{c}\underset{r=1}{\overset{m}{\Sigma }}\frac{\left(\rho N_r\right)\left(\rho N_v\right)}{{\sigma }^2\left({\tilde{A}}_i\left(t\right)\right)}{\hat{a}}_{\mathrm{ir}}^{\mathrm{true}}(t_q,t)+\frac{{\hat{a}}_{\mathrm{iv}}^{\mathrm{true}}(t_q,t)}{{\sigma }^2\left(a_{\mathrm{iv}}(t_q,t)\right)}=\frac{a_{\mathrm{iv}}(t_q,t)}{{\sigma }^2\left(a_{\mathrm{iv}}(t_q,t)\right)}+\frac{\left(\rho N_v\right){\tilde{A}}_i\left(t\right)}{{\sigma }^2\left({\tilde{A}}_i\left(t\right)\right)},\\ i=\mathrm{1,2},...,n,v=\mathrm{1,2},...,m,\end{array}\right)$

在第一次迭代中找到的近似估计被替换成(ρN_j)。

13、然后，该算法检查新值是否满足成分a_ij(t_q，t)的不确定性的 自然区域：

$|a_{\mathrm{ij}}(t_q,t)-{\hat{a}}_{\mathrm{ij}}^{\mathrm{true}}(t_q,t)|\le 3\sigma \left(a_{\mathrm{ij}}(t_q,t)\right),i=\mathrm{1,2},...,n,j=\mathrm{1,2},...,m.$

如果不符合这个条件，那么，不满足该不等式的由 最近的边界点的值取代。

由于这个原因，与迭代过程先前的步骤相比，信号Fl的值可能会在变量 的新值处更大，从而导致迭代过程中的收敛速度降低并且出现振荡。

为了在重新计算估计值${\hat{a}}_{\mathrm{ij}}^{\mathrm{true}}(t_q,t),j=\mathrm{1,2},...,m$后使信号F_l不增加，对于与 先前的迭代相比函数F_l的相应成分增大的估计值集合应当由先前步骤中的相应值取代。

14、重复上述9-13项中提到的操作，直到满足以下条件：

$\left|\right|\frac{{\left(\widehat{{\mathrm{\rho N}}_j}\right)}_l-{\left(\widehat{{\mathrm{\rho N}}_j}\right)}_{l-1}}{{\left(\widehat{{\mathrm{\rho N}}_j}\right)}_l}\left|\right|<{\gamma }_1,$

$\left|\frac{F_l-F_{l-1}}{F_l}\right|<{\gamma }_2,.$

15、为了通过少量的同位素（2-4个同位素）来识别核爆炸，对c₁和c₂上 设置节距为Δc₁和Δc₂的二维网格，c₁和c₂是对RTIC元素的独立和累积输出求 和的权重。

16、为c₁和c₂的每个值计算独立的输出向量和 ${{\eta }^{\mathrm{Pu}}}^{239}=c_2{\eta }^{{\mathrm{Pu}}_f^{239}}+(1-c_2){\eta }^{{\mathrm{Pu}}_{14}^{239}}.$

17、为每对向量和计算并保存比活度的二维信号 ${\left\{a_{\mathrm{ij}}(t_q,t)\right\}}_{c_1{c}_2},i=\overline{1,n},j=\mathrm{1,2}.$

18、为每个二维信号执行上述6-13项中提到的操作，直到满足第14项的条 件。

19、在解决了识别核爆炸的任务后在相应的二维信号找到的信 号估计之中，选择$\underset{i=1}{\overset{n}{\Sigma }}{({\tilde{A}}_i\left(t\right)-\underset{j=1}{\overset{m}{\Sigma }}{\left(a_{\mathrm{ij}}^{\mathrm{true}}(t_q,t)\right)}_{c_1{c}_2}{\left(\widehat{{\mathrm{\rho N}}_j}\right)}_{c_1{c}_2})}^2$的值为最小的估计。

20、由于取关系式$\underset{i=1}{\overset{n}{\Sigma }}{({\tilde{A}}_i\left(t\right)-\underset{j=1}{\overset{m}{\Sigma }}{\hat{a}}_{\mathrm{ij}}^{\mathrm{true}}(t_q,t)\left(\widehat{{\mathrm{\rho N}}_j}\right))}^2/\underset{j=1}{\overset{m}{\Sigma }}{\left(\widehat{{\mathrm{\rho N}}_j}\right)}^2$为最小的分离时间tq，如算 法总框图（图2）所示，方框“在此间隔内引入分离间隔[t_H，t_K]和节距Δt”对应 于第4，5项；方框“通过Tq使信号最小化”与第20项相对应。

该算法的主要部分置于方框1“搜索的估计中，其示意图如图3所 示。

方框“输入初始数据σ(a_ij(t_q，t))，γ₁，r_2”与本方法的第1，7， 8项相对应。

方框“瞬间裂变的情况下二维信号成分{a_ij(t_q，t)}的计算”对应于本方法的第 6，17项；

方框“通过c₁和c₂形成网格。形成独立的和累积的输出的平均量”对应于 第15和16项；

方框“初始化对应于本方法的第9项。

在本方法的第9-11项中提到的操作都是在方框“通过多准则的数学规划方 法搜索近似值”中执行；

方框“线性方程组求解，以获得下一个近似值”实现了第12项，方 框“用边界点值替换估计“由先 前的迭代的值替换估计值和将它们绑定的条件对应于本方法的第 13项。

方框“符合输出条件”实现了本方法的第14项，方框“选择平均输出的最优 解”实现了第19项。

本方法的产业应用

本方法在2.4GHz英特尔赛扬处理器，内存为768MB和MATLAB 7.0数学 包上的个人计算机上模拟。

模拟爆炸后6天的取样的情况，以测量5个同位素（Kr_85m，Xe_131m，Xe_133m， Xe₁₃₃，Xe₁₃₅）的核活性。在假设分离时间已知并且等于事件发生后3个小时 的条件下计算结果。

测得的核活性的值用高斯噪声叠加来“制造噪音”，根均方差等于它们 的“真”值的5％。

以下是裂变类型的组合（ρN_j变量的可能的集合）：

1）$U_{\mathrm{th}}^{235}+\mathrm{Xe}133$背景（两个不明来源）;

2）$U_{\mathrm{th}}^{235}+U_f^{235}+U_{14}^{235}$（三个不明来源）;

3）$U_{\mathrm{th}}^{235}+{\mathrm{Pu}}_f^{239}+{\mathrm{Pu}}_{14}^{239}$（三个不明来源）;

4）$U_f^{235}+U_{14}^{235}$（两个不明来源）;

5）$U_f^{235}+U_{14}^{235}+\mathrm{Xe}133$背景（三个不明来源）;

6）${\mathrm{Pu}}_f^{239}+{\mathrm{Pu}}_{14}^{239}$（不明来源）;

7）$P{u}_f^{239}+{\mathrm{Pu}}_{14}^{239}+\mathrm{Xe}133$背景（三个不明来源）;

8）$U_f^{235}+U_{14}^{235}+{\mathrm{Pu}}_f^{239}+{\mathrm{Pu}}_{14}^{239}$（四个不明来源）;

9）$U_f^{235}+U_{14}^{235}+{\mathrm{Pu}}_f^{239}+{\mathrm{Pu}}_{14}^{239}+{\mathrm{Xe}}^{133}$背景（五个不明来源），

其中是反应堆发射（反应堆上的数据是从文献中获得）。

真解是组合4，的相对贡献是100，的相对贡献是100。模拟结果 如表1所示。

“解法”行描述了如在本方法中提出的四种方法用于解决识别核爆炸的 任务（二次规划，非线性规划，阿基米德模型，优先顺序的模型），它们与 在类似的算法(Tiknonov正则法)中使用的解决方法比较。

“裂变类型组合数”行表示了裂变类型的组合的数量，该数量保证在从 总共的9组中相应的解决方法（Tiknonov正则法，多准则的数学规划方法）的 方程组（4）中差异的平方和最小。

“矩阵方程的限制性数量的数量级（Order of conditionality number for  system matrix）”行表示了方程组（4）的矩阵的限制性数量的数量级，对应 于第二行表示的裂变类型的组合。

“解估计”行表示的是在第二行表示的裂变类型组合中存在的裂变类型 贡献的估计。例如，对于Tiknonov正则法，从偏差平方和的角度来说，组合2 是最好的。三种裂变类型对应于这个组合，计算的相对贡献是35.99;计算的相对贡献是42.97,而计算的相对贡献是110.27。这也同样适用于其 他解决方法。

“偏差平方和”行表示方程组（4）的偏差平方和的值，计算该值是为了 确定裂变类型组合和它们对氪和氙的同位素的总活性的贡献。

“使用算法的时间，分钟”行表示的是通过各方法获得贡献估计所需的 时间，用分钟表示。

表1

通过各种方法解决识别任务的结果（放射性源为和真解是100 和100）

由表1可以看出，Tikhonov正则化方法给出了否定的结果，该解包括包含 了一个核反应堆的显著的相对贡献（真解中缺乏这一点）。

通过多准则规划（二次，非线性，目标规划（阿基米德模型和有优先顺 序的模型））获得了肯定的结果，使用变量的非负性的附加条件。

最佳的结果对应的是第九个裂变类型的组合。这并不与真解矛盾，因为 在真解中未显示的那些裂变类型的贡献是不显著的（其中大部分是零）。

因此，该被提议的方法的优点是：

通过以下方式提高了核爆炸参数识别的效率和准确度：利用多准则数学 规划的方法同时检查各种裂变类型组合，不需要正则化参数的计算；利用对 解的附加限制并对应于任务的物理特性（physical statement）（解的非负性， 有界性）；同时使用不同的目标函数，并通过两种能量组的中子将一种材料 中的两种裂变类型组合成为一种单一的裂变类型，也就是，减少方程组中的 未知数数量。