相较于带电粒子，空间中子的测量更复杂。首先是中子探测器的优化设计，目前国际上的中子探测器设计复杂，资源消耗多。如1998年月球勘探者号(LP)搭载伽马/中子谱仪(GS/NS)使用了多组探测单元对中子能谱进行分段测量[2]；2001年国际空间站上的邦纳多球中子谱仪(BBND)[3]使用了6个不同屏蔽厚度的球体；2012年好奇号火星车上搭载的辐射评估探测器(RAD)使用的闪烁体探测器只能对快中子探测[4]。其次是鉴别中子，排除带电粒子尤其是伽玛射线的影响。最后是中子能谱反演理论，由于中子不带电，中子探测器并不能直接测量入射中子能谱，而是测量中子进入探测器后产生的次级粒子的能量沉积谱。这样直接探测得到的能量损失与入射中子的能量并不是一一对应的，需要通过复杂的反演获得原始的入射中子能谱。由于在中子能谱反演时，存在诸如响应矩阵病态严重、响应矩阵和测量值均存在不确定性等难点，迄今为止还没有一个完全的理论上的解法。目前，在中子反演中主要采用以下几种方法：J-log法、非负最小二乘(NNLS)法、概率迭代法、遗传算法等[5-8]。
中国空间站目前计划将于5年之内建设完成。考虑到中子对宇航员的危害，需要在中国空间站开展空间热中子、超热中子和快中子等探测，服务于航天器与航天员辐射安全保障和空间科学研究。
在中国空间站上计划搭载一个新型的基于Cs2LiY16:Ce(CLYC)闪烁体的中子探测器，用于空间中子能谱的探测。CLYC闪烁体是近年来发展起来的一种新型闪烁体材料，它的密度为3.31 g/cm3, 发射峰波长为373 nm[9]。该闪烁体有很高的光产额，对于伽马射线的荧光产额高达20 000光子/MeV，其能量分辨率为4%@662 keV, 优于应用较为广泛的NaI(Tl)。与其他空间中子探测器相比，CLYC闪烁体探测中子的优势主要体现在如下3个方面：
1) ?探测效率高。CLYC闪烁体对热中子的探测效率接近100%，对能量在100 MeV以下的中子探测效率高于20%。
2) ?可探测的中子能量范围广。CLYC闪烁体对中子的探测基于核反应法，不仅可以测量快中子，也可以测量热中子以及超热中子。
3) ?中子伽马鉴别能力强。中子和伽马射线在CLYC闪烁体中引起的脉冲形状不同，利用2种粒子脉冲形状的不同可以进行中子和伽马射线的鉴别。
现阶段CLYC闪烁体尚未应用于空间探测中，在地面上常用于中子伽马射线的甄别，对中子能谱反演的研究较少。在对该中子探测器能谱反演时，本文给出了中子反演常用的NNLS法和概率迭代法[6, 8]，并针对这2种算法中存在的相对误差较大以及反演运算效率低的问题，提出了一种基于增广矩阵的非负最小二乘(AM-NNLS)法[10]，通过数值实验验证了该方法的有效性。
1 中国空间站中子探测器的探测原理 中国空间站上中子探测器是能量粒子探测器的一部分，安装在空间站试验舱I的舱外暴露平台上，对0.025 eV~100 MeV能量范围的空间中子进行探测。该中子探测器的基本设计结构如图 1所示，图中D1为主探测器(CLYC闪烁体)，用于探测中子；包围在主探测器外面的部分D2为反符合探测器(塑料闪烁体)，作为反符合探测器去除带电粒子对中子探测的影响[11]。
 |
| 图 1 反符合结构 Fig. 1 Anti-coincidence structure |
| 图选项 |
为了避免空间存在的带电粒子、伽马射线对中子探测的影响，中子探测器探测中子的逻辑工作方式如图 2所示。入射粒子进入探测器后若只在CLYC闪烁体中产生响应(触发)而在塑料闪烁体中无响应则认为该粒子为中性粒子，否则认为该粒子为无效粒子。对中性粒子，利用脉冲形状甄别(PSD)的方法进一步分辨，筛选出中子并记录其能量沉积信息。PSD算法基于的理论基础是：中子伽马射线与CLYC闪烁体相互作用的物理机制不同，引起的脉冲形状不同。粒子与CLYC闪烁体作用的物理机制有4种：自陷激子(STE)、Ce3+直接的电子空穴俘获发光、二元Vk-电子扩散以及芯价发光(CVL)[12]。不同种类粒子与CLYC闪烁体相互作用时有不同的物理机制。不同物理机制产生闪烁脉冲的上升沿和下降沿存在差别[12]。中子与CLYC闪烁体相互作用只涉及前3种物理机制，由于伽马射线与CLYC闪烁体相互作用时存在CVL，因此中子和伽马射线会在CLYC闪烁体中引起不同的脉冲形状[13]。该中子探测器所用的PSD算法在王琦标等的论文中有详细的介绍[14]。
 |
| 图 2 中子探测器的逻辑工作方式 Fig. 2 Logic schematic diagram of neutron detector |
| 图选项 |
下面将对CLYC闪烁体探测中子的原理进行分析。CLYC闪烁体中含有丰富的、对中子敏感的6Li和35Cl，根据中子在CLYC闪烁体中核反应产生的次级带电粒子的能量沉积，可以间接获取入射中子的信息。热中子(0.025 eV)到快中子能量区间内主要的2个核反应如表 1所示[15-16]。
表 1 中子与CLYC闪烁体反应过程[15-16] Table 1 Reaction of neutron with CLYC scintillator[15-16]
| 反应 | Q/MeV | 能量区间 |
| 6Li+n→3H+ɑ | +4.78 | 热中子及快中子 |
| 35Cl+n→35S+P | +0.615 | 快中子 |
表选项
中子能量较低时，以6Li(n, ɑ)3H反应为主, 与热中子反应截面高达939 barns，随着中子能量的上升，该反应截面呈减小趋势，35Cl(n, 35S)P逐渐占主导地位。利用GEANT4软件分别对0.025 eV和5 MeV的单能中子仿真，沉积能谱如图 3所示[16-17]。仿真结果显示，热中子在CLYC闪烁体中能量沉积的峰值在4.78 MeV附近，为6Li(n, ɑ)3H特征峰。当入射中子能量为5 MeV时，沉积能谱中出现较为明显的双峰，能量较高的峰来自核反应6Li(n, ɑ)3H；较低能量的峰来自核反应35Cl(n, 35S)P。
 |
| 图 3 热中子以及5 MeV的单能中子在CLYC闪烁体中的沉积能谱 Fig. 3 Deposited energy spectrum of thermal neutron and 5 MeV monoenergetic neutron in CLYC scintillator |
| 图选项 |
热中子、超热中子以及快中子均可在CLYC闪烁体中产生高于探测器的触发阈值20 keV的能量沉积，因此通过CLYC闪烁体可实现热中子到快中子的探测。当入射中子能量一定时，在探测器中产生的能量沉积非固定值，即两者不存在一一对应的关系，需要通过反演的方法获得原始的入射中子能谱。
2 中子能谱的反演 假设探测器对中子的响应函数为u，入射中子的能谱为f，测量能谱为z，那么这三者之间满足一定的映射关系：
 | (1) |
式中：i和j分别为沉积能谱和入射能谱的能道数，该方程也称为观测方程。因此，中子反演问题的本质是已知观测数据z，求方程组中未知参数f的问题。
2.1 响应函数的计算 响应函数是反演能谱的关键因素，是入射能量在某个范围内的粒子，在探测器中产生的能量沉积落在各探测能道情况的表征，与粒子入射角度、能量、探测器类型以及逻辑工作方式等有关[18]。响应函数的计算公式为[19]
 | (2) |
式中：Ed为入射粒子的能量；Eb为入射粒子在探测器中的沉积能; K用于描述入射能量为Ed、沉积能为Eb的粒子在仪器中敏感程度，K的表达式为
 | (3) |
式中:ε为一个能量在Eb和Eb+ΔEb之间的粒子从点(x, y)以天顶角θ，方位角φ入射时，在探测器上产生的沉积能在Ed和Ed+ΔEd之间的概率[18]；dF和dΩ分别为面积微元和立体角。
在中国空间站轨道，认为空间中子各向同性入射，将式(2)简化后得到该探测器的响应函数公式[18]为
 | (4) |
式中:r为入射源半径; nj为能道j的入射中子数; s为入射源面积; Nij为入射能量在第j个能道，沉积能量在第i个能道的粒子数。
利用GEANT4对该中子探测器进行蒙特卡罗仿真，仿真中使用的物理列表为QGSP_BERT_HP，模型中加入了全热中子截面库，包括20 MeV以下中子的高精度模型，可实现对热中子到快中子的仿真[17]。仿真时设置的入射源-中子的能量范围为0.025 eV~100 MeV，入射源的形状为球面源，入射方向为球面源内各向同性，按照图 2的逻辑方式，筛选出符合D1D2的中子，并记录其入射能量和在CLYC闪烁体(D1)中的沉积能量。将入射能量(0.025 eV~100 MeV)和沉积能量(20 keV~100 MeV)分别按照指数线性均分为24个能道，测量能道(沉积能道)边界为
 | (5) |
入射能道边界为
 | (6) |
得到探测器对0.025 eV~100 MeV能量范围中子的响应函数如图 4所示，不同的线型表示沉积能量对应的能道(测量能道)。可以看出，中子能量在0.025 eV~100 MeV时均能在测量能道中产生响应，这也论证了该中子探测器进行热中子、超热中子以及快中子测量的可行性。
 |
| 图 4 CLYC闪烁体中对热中子和快中子(0.025 eV~100 MeV)的响应函数 Fig. 4 Response functions of thermal neutron and fast neutron (0.025 eV-100 MeV) in CLYC scintillator |
| 图选项 |
以上仿真得到的响应函数近似于真实的响应函数，但是受泊松涨落的影响，两者有偏差。因此，在蒙特卡罗仿真过程中需要考虑泊松不确定性，并将这种不确定性添加到能谱的反演中[20]。根据式(4)的响应函数计算公式可知响应函数的误差由2部分组成：入射粒子分布的涨落和触发粒子的涨落。响应函数的误差为
 | (7) |
式中:σNij和σnj分别为触发粒子数和入射粒子数的标准差; Δij为响应函数的误差。响应函数的不确定性Δij/uij如图 5所示，不同的线型为探测能道。
 |
| 图 5 响应函数在每个能道的不确定度 Fig. 5 Uncertainty of response function with different channels |
| 图选项 |
2.2 测量值的生成 由于中国空间站上的中子探测器还未工作，无实测数据，因此需要通过正演的方法生成虚拟测量值。
空间中子能谱一般遵循幂律谱分布，假设中子入射谱符合f(E)=aE-1.2，a为幂律谱的一个无量纲参数，无具体实际含义，E为连续的能量值，则每个能道j的入射粒子强度可近似认为
 | (8) |
式中:Wj为能道j的宽度；Ej, mid为能道j上下边界能量值乘积的开方(
 | (9) |
将假定的入射谱f代入式(9)中可以确定假设的入射谱中a的值：
 | (10) |
由式(10)可以得到入射能道j中的粒子强度为
 | (11) |
将式(11)代入观测方程中即可得到理论的测量计数zi。
实际上，仪器每个能道的测量计数存在涨落，一般认为符合以λi=zi为期望值的泊松分布[21]。为了模拟反演方法对实测数据的处理效果，需要加入泊松误差后再做反演。根据计数涨落特性生成测量值Mi，每个能道的测量值Mi出现的概率为
 | (12) |
根据式(12)，以每个能道的理论测量值zi为期望生成一系列测量值，利用Bootstrap有放回随机抽样的方法产生虚拟测量值zi[22]。
2.3 反演方法比较 2.1节和2.2节在考虑不确定性的情况下，计算了响应函数和虚拟测量值，下面将基于中子在中子探测器中的响应情况，开展数值实验，对概率迭代法和NNLS法分析比较，同时提出一种AM-NNLS法的反演方法。
2.3.1 概率迭代法 概率迭代法是基于式(1)观测方程以及响应函数的物理含义得到的一种迭代算法，该算法的推导如下[8]：
沉积能道i的粒子中来自入射能道j的概率为
 | (13) |
粒子个数为
 | (14) |
 | (15) |
 | (16) |
得到的迭代关系为
 | (17) |
在不加约束的情况下，此迭代方法易出现数值溢出。为避免此情况，需对方程解的范围作约束。其计算步骤如图 6所示。反演得到的中子能谱如图 7中蓝色实线所示。纵坐标为入射能谱的微分强度，该结果与原始谱的符合性较好，但反演运算效率低，耗时156 s。
 |
| 图 6 概率迭代法的反演步骤 Fig. 6 Inversion steps of probabilistic iterative method |
| 图选项 |
 |
| 图 7 基于概率迭代法和NNLS法的中子微分谱 Fig. 7 Neutron differential spectrum based on probabilisticiterative method and NNLS |
| 图选项 |
2.3.2 NNLS法 最小二乘法的本质是求解使残差向量的平方和极小化的解。
 | (18) |
最小二乘法求解方程组(尤其系数矩阵为病态时)的解常包含负数，对于能谱，这样的解不具有物理意义，因此需要约束解的范围(fj>0)，这种方法称为NNLS法[23-24]。利用L-BFGS-B迭代[25]方法对式(18)求解，耗时9.9 s，反演运算效率高于概率迭代法。反演结果如图 7中紫色实线所示。对比图 7中2种反演结果，可知NNLS法运算效率高于概率迭代法，但是反演结果的可信度低于概率迭代法，这是由于响应函数病态性对NNLS法的影响更大。
2.3.3 AM-NNLS法 为了有效地解决响应矩阵的病态性导致的NNLS解相对误差偏大的问题，提出了AM-NNLS法，该算法的思想是在NNLS法中引入胡圣荣和戴纳新的增广方程组的思想[10]，将式(1)改写为
 | (19) |
式中:I为单位矩阵；令
 | (20) |
此时，中子能谱f只是式(20)的局部解，将式(20)作为目标函数，得到的中子微分谱如图 8中紫色实线所示，真实入射能谱与反演能谱的相对误差如图 9中虚线所示。其中蓝色实线为由NNLS反演得到的中子能谱与真实入射谱之间的相对误差，绿色虚线为基于AM-NNLS法反演得到的中子能谱的相对误差。从图中可以看出，AM-NNLS法反演得到的能谱与真实能谱的符合性优于NNLS法，相对误差明显减小。此外AM-NNLS法的反演所需的时间为5 s，比NNLS法的运算效率提高了50%左右。
 |
| 图 8 NNLS法与AM-NNLS法的能谱对比 Fig. 8 Comparative energy spectrum of non-negative least squares method and AM-NNLS |
| 图选项 |
 |
| 图 9 相对误差 Fig. 9 Relative error |
| 图选项 |
2.4 讨论 由3种方法的反演结果可知，概率迭代法和AM-NNLS法的反演误差小，结果可信度高，适合基于中国空间站中子探测器的中子能谱的反演。为了检验测量值改变时，这2种反演方法的有效性，在不同测量值条件下对2种反演方法进行数值实验。
当测量值N=102，103，104，105时，用2种方法得到的反演谱如图 10所示。反演得到的能谱与真实入射能谱之间的相对误差比较如图 11所示，结果显示随着测量值的增大，反演相对误差明显减小。在1 keV以下能量点和真实能谱符合度高；在1 keV以上时反演相对误差明显增大。入射能量在1 eV~1 MeV时，AM-NNLS法解谱得到的相对误差小于概率迭代法，但当入射能量低于1 eV或高于1 MeV时，这种优势不再明显。
 |
| 图 10 不同测量总计数时反演得到的能谱 Fig. 10 Inverted energy spectra for different measurements |
| 图选项 |
 |
| 图 11 测量总计数不同时2种反演方法得到的能谱与真实入射谱之间的相对误差 Fig. 11 Relative error between energy spectrum obtained from two inversion methods and actual input spectrum when measured values are different |
| 图选项 |
2种方法反演运算时间如表 2所示。可知，AM-NNLS法反演运算效率高于概率迭代法，通过AM-NNLS法可以用更短的运算时间获得较为精确的中子。
表 2 概率迭代法与AM-NNLS法反演能谱所需时间对比 Table 2 Time comparison between probabilistic iterative method and AM-NNLS for inversion of energy spectrum
| 测量总计数 | 耗时/s | |
| 概率迭代法 | AM-NNLS | |
| 102 | 138 | 5 |
| 103 | 143 | 4 |
| 104 | 156 | 5 |
| 105 | 150 | 5 |
表选项
3 结论 CLYC闪烁体是目前国际上最新的中子探测器材料，其在仪器资源消耗、中子探测效率、中子/伽玛射线鉴别等方面都具有一定的优势。在未来的空间中子测量应用中将具有广大的应用前景。
本文介绍了中国空间站上基于CLYC闪烁体的中子探测器的探测原理，通过GEANT4仿真证明了该空间中子探测器测量热中子、超热中子以及快中子的可行性。针对中子能谱反演，提出了AM-NNLS法。在考虑响应函数和测量值不确定度的情况下，利用概率迭代法、NNLS法以及AM-NNLS法进行数值实验。3种方法对中子反演结果进行了对比，结果表明：
1) ?反演误差。AM-NNLS法得到的反演误差比NNLS法少93%，比概率迭代法降低11%。
2) ?反演运算效率。AM-NNLS法比NNLS法的解谱时间缩短了50%左右，比概率迭代法缩短了90%以上，这对空间站实时/准实时的数据反演提供了保障。
3) ?在测量计数增大时, AM-NNLS法和概率迭代法的反演相对误差均减小，因此为了得到更加精确的反演结果，在实际应用中，可以适当地增加探测结果的累积时间获得更多的探测计数。
中子探测器反演时，使用的能道划分方法并不唯一，实际应用中针对探测结果以及探测器的工作状态，可以灵活划分能道，如确定哪一个测量范围作为反演对象、能道数划分为多少最合适等都是需要进一步研究的内容。
致谢 感谢中国科学院国家空间科学中心空间站能量粒子探测器研制团队对本文工作的贡献。
参考文献
| [1] | MOUNTFORD P J, TEMPERTON D H. Recommendations of the international commission on radiological protection(ICRP)1990[J]. European Journal of Nuclear Medicine, 1992, 19(2): 77-79. |
| [2] | MAZUR J E, CRAIN W R, LOOPER M D, et al. New measurements of total ionizing dose in the lunar environment[J]. Space Weather-The International Journal of Research & Applications, 2011, 9(7): 1-12. |
| [3] | KOSHIISHI H, MATSUMOTO H, CHISHIKI A, et al. Evaluation of the neutron radiation environment inside the International Space Station based on the bonner ball neutron detector experiment[J]. Radiation Measurements, 2007, 42(9): 1510-1520. DOI:10.1016/j.radmeas.2007.02.072 |
| [4] | HASSLER D M, ZEITLIN C, WIMMER-SCHWEINGRUBER R F, et al. The radiation assessment detector (RAD) investigation[J]. Space Science Reviews, 2012, 170(1-4): 503-558. DOI:10.1007/s11214-012-9913-1 |
| [5] | TANIGUCHI T, UEDA N, NAKAZAWA M, et al. Systematic study on spectral effects in the adjustment calculations using the NEUPAC-83 code[M]. Berlin: Springer, 1985: 685-691. |
| [6] | KOHLER J, EHRESMANN B, MARTIN C, et al. Inversion of neutron/gamma spectra from scintillator measurements[J]. Nuclear Instruments and Methods in Physics Research B, 2011, 269: 2641-2648. DOI:10.1016/j.nimb.2011.07.021 |
| [7] | 王冬, 何彬, 张全虎. 用遗传算法求解中子能谱[J]. 原子能科学技术, 2010, 44(10): 1270-1275. WANG D, HE B, ZHANG Q H. Unfolding neutron spectrum using genertic algorithm[J]. Atomic Energy Science and Technology, 2010, 44(10): 1270-1275. (in Chinese) |
| [8] | 杨鑫, 李润东, 刘汉刚, 等. 基于概率迭代的NDP反演方法[J]. 计算物理, 2012, 29(6): 891-900. YANG X, LI R D, LIU H G, et al. An unfolding method of NDP based on probability iteration[J]. Chinese Journal of Computational Physics, 2012, 29(6): 891-900. DOI:10.3969/j.issn.1001-246X.2012.06.014 (in Chinese) |
| [9] | D'OLYMPIA N W, NATHAN W.Development of novel neutron and gamma-ray scintillators: Cs2LiYCl6: Ce and CeBr3[D].Boston: University of Massachusetts, 2013: 8-40. |
| [10] | 胡圣荣, 戴纳新. 病态线性方程组的新解法:增广方程组法[J]. 华南农业大学学报, 2009, 30(1): 119-121. HU S R, DAI N X. A novel method for solving Ⅲ-conditioned liner system:Augmented system method[J]. Journal of South China Agricultural University, 2009, 30(1): 119-121. (in Chinese) |
| [11] | 李肖.反符合探测杯在空间粒子探测中的应用[D].北京: 中国科学院, 2015: 5-23. LI X.The application of anti-coincidence detective cup technology in space particle detection[D].Beijing: University of Chinese Academy of Science, 2015: 5-23(in Chinese). |
| [12] | WHITNEY C M, SOUNDARA-PANDIAN L, JOHNSON E B, et al. Gamma-neutron imaging system utilizing pulse shape discrimination with CLYC[J]. Nuclear Instruments and Methods in Physics Research Section A:Accelerators, Spectrometers, Detectors and Associated Equipment, 2015, 784: 346-351. DOI:10.1016/j.nima.2014.09.022 |
| [13] | LEE D W, STONEHILL L C, KLIMENKO A, et al. Pulse-shape analysis of Cs2LiYCl6:Ce scintillator for neutron and gamma-ray discrimination[J]. Nuclear Instruments and Methods in Physics Research Section A:Accelerators, Spectrometers, Detectors and Associated Equipment, 2012, 664(1): 1-5. |
| [14] | WANG Q B, TUO X G, DENG C, et al. Characterization of a Cs2LiYCl6:Ce3+ scintillator coupled with two silicon photomultiplier arrays of different sizes[J]. Nuclear Instruments and Methods in Physics Research Section A:Accelerators, Spectrometers, Detectors and Associated Equipment, 2019, 942: 1-5. |
| [15] | D'OLYMPIA N W, CHOWDHURY P, LISTER C J, et al. Pulse-shape analysis of CLYC for thermal neutrons, fast neutrons, and gamma-rays[J]. Nuclear Instruments and Methods in Physics Research Section A:Accereators, Spectrometers, Detectors and Associated Equipment, 2013, 714: 121-127. |
| [16] | BROWN D, CHADWICK M B, CAPOTE R, et al. ENDF/B-Ⅷ.0:The 8th major release of the nuclear reaction data library with CIELO-project cross sections, new stantards and thermal scattering data[J]. Nuclear Data Sheets, 2018, 148(2): 1-142. |
| [17] | AGOSTINELLI S, ALLISON J, AMAKO K, et al. GEANT4-A simulation toolkit[J]. Nuclear Instruments and Methods in Physics Research Section A:Accelerators, Spectrometers, Detectors and Associated Equipment, 2003, 506(3): 250-303. |
| [18] | ZHANG S Y, ZHANG X G, WANG C Q, et al. The geometric factor of high energy protons detector on FY-3 satellite[J]. Science China Earth Sciences, 2014, 57(10): 2558-2566. DOI:10.1007/s11430-014-4853-0 |
| [19] | KHARYTONOV A, B?HM E, WIMMER-SCHWEINGRUBER R F. Regularization methods used in error analysis of solar particle spectra measured on SOHO/EPHIN[J]. Astronomy & Astrophysics, 2009, 495(2): 663-675. |
| [20] | GUO J N, BANJAC S, ROSTEL L, et al. Implementation and validation of the GEANT4/AtRIS code to model the radiation environment at Mars[J]. Journal of Space Weather and Space Climate, 2019, 9(3): 1-30. |
| [21] | NOWAK R D, KOLACZYK E D. A statistical multiscale framework for Poisson inverse problems[J]. IEEE Transactions on Information Theory, 2000, 46(5): 1811-1825. DOI:10.1109/18.857793 |
| [22] | WILLIAM H P, BRIAN P F, SAUL A, et al. Numerical recipes:The art of scientific computing[M]. Cambridge: Cambridge University Press, 2007: 809-816. |
| [23] | LAWSON C L, HANSON R J. Solving least squares problems[M]. Philadelphia: Society for lndustrial and Applied Mathematics, 1987: 1-5. |
| [24] | B?HM E, KHARYTONOV A, WIMMER-SCHWEINGRUBER R F. Solar energetic particle spectra from the SOHO-EPHIN sensor by application of regularization methods[J]. Astronomy & Astrophysics, 2007, 473(2): 673-682. |
| [25] | MORALE J L, NOCEDAL J. Remark on 'algorithm 778:L-BFGS-B:Fortran subroutines for large-scale bound constrained optimization'[J]. ACM Transactions on Mathematical Software, 2011, 38(1): 1-4. |
