最優化方法應用于裂變核反應理論計算的研究

田 源，續瑞瑞，張 玥，陶 曦，王記民，金永利，張 智，孫小東，葛智剛

(中國原子能科學研究院 核數據重點實驗室，中國核數據中心，北京 102413)

核數據是核能、核工程、核技術應用和核科學研究的基礎數據，雖然實驗測量是獲得核數據的基礎，然而，受各種因素的限制(如實驗成本高，更重要的是在實驗上所要求的入射粒子種類和能量、靶核制備以及某些測量條件等方面并不總能實現)，因此，在核裝置的設計與測試過程中，實驗常被用來精確定量少量關鍵核素的重要反應與結構核數據，而絕大部分核工程所需數據須通過理論計算完成。要使用核反應模型計算核數據，一方面需盡可能準確描述相應核反應的反應機制，另一方面則是快速有效確定核反應模型中的參數。本項目圍繞裂變核反應模型，利用先進的最優化方法確定模型中的參數，并應用于錒系核中子誘發裂變反應核數據的研究。

目前國際上最著名的核反應模型程序系統有TALYS[1]和EMPIRE[2]等。我國采用的FUNF[3]是由中國核數據中心張競上研究員開發研制，具有我國自主知識產權的裂變核反應程序系統，可用于計算和分析能量在20 MeV以下中子誘發易裂變核反應相關的數據。多年來，一直在我國核數據理論研究過程中發揮著重要的作用。FUNF的理論框架包括光學模型、Hauser-Feshbach模型以及激子模型等。裂變核反應模型總的可調參數多達50個。因此，需通過一些先進的最優化方法，根據已知的核反應實驗結果，確定裂變核反應模型的參數，并最終用于計算裂變核反應數據。

最優化方法是一種數學方法，它是研究在給定約束條件下如何尋求某些因素(或量)，以使某一(或某些)指標達到最優的一些方法的總稱。本工作要在裂變核反應模型的約束條件下，通過不同的調參方法調整模型中的參數，使得模型計算的裂變核反應數據，如反應道的截面、出射中子的角分布和能譜、各種出射粒子的雙微分截面等，與已知的實驗結果的符合程度達到最優。目前較常用的的調參方法有梯度下降法[4]、牛頓法[5]、共軛梯度法[6]以及啟發式優化算法等。這些方法各有優缺點，分別適用不同的調參環境，適當的調參方法能大幅提高參數調優的計算成本。20世紀70年代，歐洲核子研究中心(CERN)的科學家Fred James綜合了多種最優化方法建立了尋找方程最小值的物理分析工具MINUIT[7]。迄今為止，MINUIT方法已經被廣泛用于分析實驗數據和建立理論模型。美國的EMPIRE程序將MINUIT方法用于參數調節。

對于簡單的核反應模型，參數一般較少，計算也相對簡單，一般的調參方法如格點法、隨機法通過不斷的嘗試即可找到最優的參數組。FUNF裂變核反應模型不僅約有50個參數，計算也相對復雜，單次計算時間隨著入射粒子能量的不同，可能是幾秒也可能多達數分鐘。因此本工作將FUNF裂變核反應模型與MINUIT方法相結合，并引入并行計算等方法，研究如何將最優化方法用于調節FUNF裂變核反應模型中的參數，并改進裂變核反應數據。

1 裂變核反應理論模型

當中子入射能量在20 MeV范圍內時，核子與核發生反應的全過程可被劃分為3個階段：獨立粒子階段、復合系統階段以及發射粒子后形成剩余核階段。獨立粒子階段中，入射粒子在靶核的勢場作用范圍內發生散射與吸收，可采用光學模型進行描述；進入復合核階段中，入射粒子與靶核相互作用，隨著入射粒子能量的變化，入射粒子與靶核中核子發生能量交換的可能性非常豐富多樣，該階段可發生直接反應、復合核粒子發射反應以及預平衡粒子發射反應，因此，發展了包括直接反應、蒸發模型、Hauser-Feshbach模型、激子模型等一系列核反應理論來描述該階段，此外，中子與錒系核相互作用的過程中，需特別考慮核裂變的貢獻對核反應數據的影響；最后一個階段，在上一個反應階段的基礎上，復合核系統發射粒子后變成剩余核。張競上研究員研制的中子與錒系核反應計算程序FUNF，包含了快中子能區內錒系核的主要核反應過程，考慮3次粒子發射過程，多年來已成功為核工程用戶解決多項錒系核反應數據問題，本工作中錒系核反應中子數據模型計算工作由該程序完成。

1.1 光學模型

光學模型是核反應過程的基礎與起點，控制著之后的多項核反應物理過程，該模型的引入成功描述了入射粒子在靶核勢場中進行的散射與吸收的過程。光學模型所建立的復數形式勢場，結合求解薛定諤方程可用來描述散射過程中的多個物理量，包括全截面、反應(去彈)截面、彈性散射截面及角分布和穿透系數，并為復合核計算提供重要的穿透系數。

在FUNF程序中，中子與錒系核全套核反應數據計算中，光學模型勢采用對實驗數據描述準確度較高的唯象光學勢[8]，對應的一般形式如下：

U(r)=Vr(r,E)+iWs(r,E)+iWv(r,E)+

[Vso(r)+iWso(r)](s·l)+Vc(r)

(1)

其中：Vr(r,E)為實部中心勢；Ws(r,E)為虛部勢的面吸收項；Wv(r,E)為虛部勢的體吸收項；Vso(r)為自旋軌道耦合勢的實部；Wso(r)為自旋軌道耦合勢的虛部；s、l分別為入射粒子自旋角動量與軌道角動量；Vc(r)為庫侖勢；E為入射中子在實驗室系下的能量，MeV。

式(1)中，實部中心勢可表示為：

(2)

式(2)中V(E)可表示為能量依賴關系：

V(E)=V0+V1E+V2E2+

(3)

其中：A為靶核質量數；N為靶核中子數；Z為靶核電荷數。

對應虛部勢的面吸收形式為：

(4)

(5)

對應虛部勢的體吸收形式為：

(6)

Wv(E)=max{0,U0+U1E}

(7)

對應自旋軌道耦合勢實部Vso(r)和虛部Wso(r)分別表示為：

(8)

(9)

其中：ar、as、av和aso為光學模型勢的對應項中的彌散寬度參數；Rr、Rs、Rv和Rso為光學模型勢的對應項中的半徑參數，單位為fm，對應計算公式為：

(10)

其中，i=r,s,v,so，對應半徑參數。

此外，庫侖勢對應形式為：

(11)

式中：z、Z為入射粒子和靶核的電荷數；Rc為庫侖半徑。因此，可調的光學勢參數為12個深度參數、5個半徑參數、4個彌散寬度參數，其中有12個參數對截面較為敏感。

1.2 復合核模型

當入射粒子被靶核吸收后，核反應進入第2個階段，即復合核階段。在該過程中，入射粒子與靶核發生多次相互作用，不斷將能量傳遞給靶核，最終達到統計平衡的狀態，與靶核融合形成新的原子核，被稱為復合核過程，其作用時間最長，約10-15s。進入平衡態的復合核通過發射粒子而衰變，剩余核處于分立能級的激發態，實驗表明，截面與原子核的初態和末態的激發能、自旋、宇稱等性質直接關聯，考慮了角動量守恒與宇稱守恒的Hauser-Feshbach理論被發展用來描述該部分貢獻，結合寬度漲落修正可更好描述真實的實驗測量。

當入射粒子能量增大，二次核反應過程開放，預平衡核反應過程變得不可忽略。該過程中粒子與靶核作用時間介于直接過程與平衡過程之間(約10-22～10-16s)，粒子與靶核中少量核子發生能量交換，在體系還未達到統計平衡狀態時就發射粒子退激。因此，Griffin等在Feshbach的中間結構理論基礎上提出了激子模型用以描述不同激子態條件下的粒子發射過程。

在上述兩種理論體系基礎上，張競上等發展的統一的Hauser-Feshbach與激子模型在保證了角動量守恒與宇稱守恒的前提下，可考慮預平衡粒子道分立能級態的貢獻，對于準確計算中子能譜等物理量有重要價值，并成功應用在裂變核反應程序FUNF中。

一般，復合核過程中考慮的核反應包括中子輻射俘獲(n,γ)和中子非彈(n,inl)、(n,p)、(n,d)、(n,t)、(n,α)、(n,2n)、(n,3n)等，這些核反應互相競爭。在實際工作中，為考慮各類出射粒子反應道的發射概率，除輸入各復合核低激發態分立能級、自旋、宇稱等信息外，復合核相關模型中還引入能級密度、對修正等參數考慮高激發態的核結構與核內核子運動形態，從而確定各反應道的分配比例。另一方面，與中重核反應不同，中子與錒系核相互作用的核反應還需考慮核裂變截面(n,f)的貢獻，該反應同樣與其他核反應道競爭。FUNF程序中，通過唯象引入3個彼此獨立的裂變位壘來表征整個核反應中的裂變系統，基于Bohr-Wheeler公式，通過調整各裂變位壘的高度、曲率形狀、壘上能級密度、對修正等參數，確定裂變位壘形狀，計算得到裂變概率，在與其他反應道競爭的條件下，共同確定(n,f)、(n,nf)、(n,2nf)3個裂變道截面的大小。

因此，復合核模型中涉及到的可調參數包括11個反應道的能級密度、對修正參數，其中7個反應道的參數靈敏度較高，共計14個；8個反應道的巨共振參數，其中2個反應道的參數較敏感，共計12個；3個裂變道對應的位壘高度、曲率形狀、壘上能級密度、對修正等參數，共計12個。合計至少需要調整38個參數用于計算各反應道的截面數據。

1.3 直接核反應

形成復合核并非是所有核反應的必經階段。當入射粒子能量較高時，入射粒子與靶核通過很短的相互作用時間(約10-19～10-21s)發生反應，具體過程包括：1) 入射粒子將能量傳遞給靶核表面或者內部的1個或多個核子，并逃出該復合系統；2) 入射粒子將能量傳遞給靶核后飛出該系統，靶核得到部分能量而產生集體激發，并發生集體轉動或振動；3) 入射粒子與靶核發生多次碰撞，最終逃出該復合系統。

在直接反應模型中，通常采用扭曲波玻恩近似(DWBA)與耦合道方法對直接反應的貢獻進行計算，其中DWBA用來描述球形核或近球形核，耦合道方法用來描述變形核。在中子與錒系核反應計算中，直接核反應的貢獻在中子非彈截面中貢獻較大，由于錒系核受到激發后，普遍存在集體轉動能級，因此，一般采用耦合道方法來計算該部分工作，本工作將采用ECIS程序完成直接非彈耦合計算的任務[9]。

2 最優化方法

2.1 MINUIT方法

MINUIT是由CERN的物理學家Fred James在20世紀70年代建立的用于尋找方程最小值的物理分析工具。MINUIT中已包含了多種先進的最優化方法，如蒙特卡羅隨機法、單純形算法以及變尺度法等。經過多年的發展，這套工具已有Fortran、C++、JAVA以及Python等多個版本，在多個領域被廣泛用于調參并取得了巨大的成功。

MINUIT方法中不同的最優化方法可分別針對不同調參環境做出選擇。

1) 蒙特卡羅方法[10]，這種方法通過蒙特卡羅隨機抽樣的方法尋找參數。主要用于調參初期參數初始值不確定的情況下，通過隨機掃點找到可能存在的多個極小點。

2) 可變單純形算法[11]，這種方法又名Nelder-Mead方法，是常見的直接搜索型非線性優化方法，由Nelder和Mead在1965年提出。此方法利用單純形進行搜索(單純形是N維空間的廣義三角形)?？勺儐渭冃嗡惴ㄔ诿總€頂點對用戶提供厄函數進行求值，然后發現更好的點來迭代收縮單純形。當到達預期的邊界或其他終止條件時，方法終止。由于未利用任何求導運算，算法較簡單，但收斂速度較慢，適合變元數不是很多的方程求極值。

3) MIGRAD方法[12]，MIGRAD方法是基于Fletcher提出的變尺度法，它克服了梯度法收斂速度慢和牛頓法計算量、存儲量大的特點，被公認為是求解無約束優化問題最有效的算法之一。

2.2 MINUIT用于FUNF參數優化

圖1為MINUIT方法用于調整FUNF裂變核反應程序參數的流程圖。第1步，將經驗參數作為初始參數代入到FUNF裂變核反應模型中得到1組裂變核反應數據，通過比較得到與已知的實驗結果偏差。第2步，通過下山單純形或MIGRAD變尺度法找到1組參數，再代入到模型中，計算出1個新的實驗數據的偏差。若偏差有所減小，則接受這套新的參數；若偏差增大，則利用下山單純形或MIGRAD變尺度法重新找到1組參數。重復第1、2步，直到滿足FUNF裂變核反應模型計算的結果與實驗數據的偏差小于預設的條件，這時得到參數組則是需1個極小解。利用蒙特卡羅法隨機生成1組新的參數組，重復第1、2步，得到1組新的極小解。反復找到幾組極小值后，比較極小值最終確定FUNF裂變核反應的最優解。

圖1 MINUIT用于調整FUNF裂變核反應程序參數的流程圖Fig.1 MINUIT flowchart for adjusting parameters of FUNF fission nuclear reaction program

2.3 參數優化準則

首先，在調參過程中，判斷參數是否合格的1個重要結果就是最小均方根：

(12)

其中：σExp為需要符合的實驗結果；σth為理論計算的結果；σerr為誤差權重。對于裂變核反應，首先可通過調研已有實驗數據和評價數據，確定各反應道的截面數據σExp。不同的實驗測試之間可能存在分歧，且不同實驗給出的誤差范圍也不盡相同。σerr的選取決定了理論計算趨向于哪個實驗結果，因此這個量的選取尤為重要。需要結合實驗數據的復雜程度、實驗的誤差范圍以及其他評價數據庫的評價結果和評價誤差最終確定調參過程中采用的誤差權重。

其次，在FUNF裂變核反應程序中，光學模型是所有計算的基礎，直接反應、復合核反應和預平衡反應均依賴光學模型的計算結果。因此在調參順序上，首先根據總截面、彈性截面、彈性角分布等數據確定光學模型的參數；接著固定光學模型參數，再根據各反應道的數據信息同時自動調整能級密度、對修正、巨共振以及裂變位壘等參數；最后在所有截面與實驗數據基本一致的情況下，再對個別參數通過提高權重，改變參數范圍、步長的方法進一步精調。最終得到最優化參數組和理論計算的數據結果。

判斷參數優化要遵循以下3個判據。1) 參數是否具有物理意義，FUNF裂變核反應程序中的所有參數有各自物理意義，因此每個參數均只能在合理的范圍內變化，超出這個范圍會導致非物理的結果出現。2) 是否得到χ2值最優。由于涉及參數較多，系統χ2可能存在多個極小值，需多次嘗試調參得到最優解。3) 與實驗結果相比，理論計算結果是否具有好的視覺符合。由于實驗結果較多，且結構相對復雜，需通過觀察判斷理論計算的各反應道的截面是否能較好符合實驗結果。

2.4 并行化調參

FUNF模型由多個核反應模型組成，結構相對較復雜，計算時間較長，一般在5～10 s之間。涉及到的參數可能多達50個。因此即使在MINUIT方法的輔助下，單核的調參時間也可能高達數天。若同時對多個核素進行擬合，調參時間將成倍增長，因此需將FUNF+MINUIT的調參程序進行并行化加速，從而縮短時間，提高效率。消息傳遞接口(message passing interface, MPI)是科學計算中常用的編程模型。在這種并行編程中，每個控制流均有自己獨立的地址空間，不同的控制流之間不能直接訪問彼此的地址空間，須通過顯式的消息傳遞來實現。這種編程方式是大規模并行處理機(MPP)和機群(Cluster)采用的主要編程方式。由于消息傳遞程序設計要求用戶很好地分解問題，組織不同控制流間的數據交換，并行計算粒度大，特別適合于大規模可擴展并行算法。計算過程中發現，通過將MPI與MINUIT相結合，隨著CPU線程數量的增加，FUNF裂變核反應程序的調參時間成倍減小，計算效率成倍增加。

2.5 最優化方法的選擇

圖2為用可變單純形算法(SIMPLX)和變尺度法(MIGRAD)對光學模型參數(OPM)和復合核模型參數(CM)優化的迭代次數。多核并行的迭代次數為每個線程的迭代次數的平均值。迭代次數越少，說明最優化方法能更快的找到參數的最優解。復合核模型的參數有38個，大于光學模型的12個參數，因此需花費更多的迭代次數得到最優化參數組。無論是針對光學模型調參，還是調整復合核模型的參數，變尺度法均好于可變單純形算法。另外，通過MPI方法實現了并行化調參。在包含28個線程的小型服務器幫助下，調參效率提高了約20倍左右，且計算效率還會隨著線程的增加繼續提高。因此在計算過程中，主要采用并行的變尺度法(MPI-MIGRAD)進行調參，只有在變尺度法難以收斂時才考慮采用其他的最優化方法輔助調參。

圖2 不同最優化方法在單核和多核并行的迭代次數Fig.2 Number of iterations of different optimization methods in single-core and multi-core parallelisms

3 結果與討論

利用最優化方法調整FUNF裂變核反應程序的相關參數，并得到與中子誘發238U核反應的主要反應截面相符的理論計算結果。

首先，通過參考EXFOR實驗核反應數據庫和ENDF核反應評價庫中快中子區的截面結果收集中子誘發238U的核反應數據。FUNF裂變核反應程序能計算的反應道包括(n,tot)、(n,γ)、(n,inl)、(n,2n)、(n,3n)、(n,f)、(n,p)、(n,d)、(n,t)和(n,α)等。目前實驗測量主要集中在(n,tot)、(n,inl)、(n,2n)、(n,3n)、(n,γ)和(n,f)這幾個反應道中。而帶電粒子出射道中的(n,d)和(n,t)反應道無實驗數據，(n,p)和(n,alpha)反應道的實驗數據較少，無法用于指導參數調優工作。(n,3n)反應道的實驗數據也相對較少。因此，本工作的參數調優主要依托于(n,tot)、(n,γ)、(n,inl)、(n,2n)和(n,f)的實驗數據引導，并參考ENDF評價庫的部分結果。

然后，將MINUIT方法與FUNF裂變核反應程序相結合。先根據總截面、彈性角分布以及非彈截面確定光學模型參數；然后通過擬合各反應道的截面結果，得到能級密度、對修正、巨共振以及裂變位壘參數值；最后針對個別參數進行精調，得到1組可再現中子誘發238U核反應截面數據的FUNF參數組。

表1為通過擬合總截面、彈性散射角分布以及非彈截面并使得χ2最小，最終得到的光學模型參數。其中ai為彌散半徑，ri為半徑，U、V、W分別為BG光學勢的強度參數，i=r、s、v、so、c分別表示實部、表面虛部、體虛部、自旋軌道耦合項以及庫侖項。觀察發現，各參數均在合理的物理區間范圍內。圖3為238U(n,tot)反應截面，理論計算結果與實驗測量數值和評價數據基本符合。說明得到的光學模型參數是適用的。

表1 FUNF采用的光學模型參數Table 1 Parameters of optical potential model of FUNF

圖3 238U(n,tot)截面實驗數據及評價數據Fig.3 Experimental data and evaluation data of 238U(n,tot) cross sections

表2為FUNF計算程序中的參數，其中LD為能級密度參數，PC為對修正參數，FB為裂變位壘參數，后面的數字0～10分別表示(n,γ)、(n,inl)、(n,p)、(n,α)、(n,d)、(n,t)、(n,2n)、(n,3n)、(n,f)、(n,nf)和(n,2nf)這些反應道；CP表示裂變位壘曲率半徑，CFLD表示裂變能級密度系數，后面的數字0～2表示裂變道中的(n,f)、(n,nf)和(n,2nf)這3個反應道。圖4～8分別展示了(n,inl)、(n,2n)、(n,3n)、(n,f)和(n,γ)反應道的理論計算結果與實驗數據及評價數據的比較。圖中各種不同顏色符號表示實驗測量數據，不同顏色線段代表不同的評價庫的評價數據，粗紅線為FUNF裂變核反應程序計算的結果。調參最終得到所有反應道的χ2之和最小。觀察發現，各參數均在合理的物理區間內，且理論計算的結果與實驗數據和評價數據保持一致。因此，通過MINUIT最優化方法得到的FUNF裂變核反應程序用于計算20 MeV以內中子誘發238U核反應數據的參數合理有效。

表2 FUNF采用能級密度、對修正以及裂變位壘參數Table 2 Level density, pairing and fission barrier parameters of FUNF

圖4 238U(n,inl)截面實驗數據及評價數據Fig.4 Experimental data and evaluation data of 238U(n,inl) cross sections

圖5 238U(n,2n)截面實驗數據及評價數據Fig.5 Experimental data and evaluation data of 238U(n,2n) cross sections

圖6 238U(n,3n)截面實驗數據及評價數據Fig.6 Experimental data and evaluation data of 238U(n,3n) cross sections

4 總結

對于簡單的核反應模型，參數一般較少，計算也相對簡單，一般的調參方法如格點法、隨機法通過不斷的嘗試即可找到最優的參數組。裂變核反應的反應機制十分復雜，需由幾個不同核反應模型共同描述，因此需同時對大量參數進行優化才能正確描述裂變核反應數據。

圖7 238U(n,f)截面實驗數據及評價數據Fig.7 Experimental data and evaluation data of 238U(n,f) cross sections

圖8 238U(n,γ)截面實驗數據及評價數據Fig.8 Experimental data and evaluation data of 238U(n, γ) cross sections

在FUNF裂變核反應模型基礎上，首先將最優化方法MINUIT與FUNF相結合，實現多參數的自動調參。FUNF裂變核反應程序由球形光學模型、統一的Hauser-Feshbach和激子模型理論等核反應模型共同構成。FUNF涉及到的參數包括粒子的光學勢參數、巨共振參數、能級密度參數、對修正參數以及裂變位壘和位壘曲率等參數，多達50個。FUNF程序的計算時間也相對較長，所需調參時間可能長達數天。為提高調參效率，采用MPI的并行計算方法，在多核心并行計算的幫助下，將計算效率提高了20倍左右，且計算效率還可隨線程數的增加繼續提高。

在這些工作的基礎上，針對238U的裂變核反應進行了初步的計算。結果顯示，通過MINUIT最優化方法得到的FUNF參數能很好描述238U的裂變核反應數據。后期將繼續優化程序，進一步提高調參的效率和參數的精確度。并對相關的裂變核反應展開應用，提高我國裂變核反應數據的精度。