全局減方差方法在大空間γ輻射場計算中的應用

劉利 左應紅 牛勝利 朱金輝 商鵬 王學棟

（西北核技術研究所 西安 710024）

γ射線在近地面大空間長距離輸運中會與空氣和地面土壤發生多次散射，使得輸運至距輻射源公里及以上范圍內的γ粒子數較少，致使該處的蒙特卡羅模擬統計漲落較大。左應紅等［1］建立了柵元權窗結合密度逼近迭代的局域減方差方法（Local Variance Reduction，LVR）［2-3］，能夠求解離源5 km的觀測點處的γ射線注量。而在核設施廠房輻射環境［4］、早期核輻射環境［5］、核電磁脈沖環境［5］、人員輻射劑量學等方面的研究中，常常需要計算大空間中γ輻射場的整體分布。傳統的LVR方法通常只針對特定目標探測器，難以給出理想的全空間γ輻射場分布。

相比之下，全局減方差方法（Global Variance Reduction，GVR）［6-9］更適于求解輻射場的整體分布。GVR方法通過全局權窗來控制模擬粒子的空間分布，使得模擬結果的相對偏差隨空間的分布盡可能減小，從而整體提高全局空間內的計算效率。Davis和Andrew等［6-7］提出的基于蒙特卡羅正向計算的GVR方法，通過蒙特卡羅模擬計算得到粒子通量、數量和權重等粒子分布信息，用于產生基于柵元/網格的全局權窗參數，利用全局權窗參數開展蒙特卡羅模擬，循環迭代直至結果收斂。該方法思路清晰，適用性強，已廣泛用于反應堆裝置［10-12］、托卡馬克裝置［13-14］和厚屏蔽體［15-16］等復雜模型的全局輻射環境計算。

大空間γ輻射場模擬具有空間尺度跨度大和計數區外散射較強等特點。直接應用已有的GVR方法面臨計數柵元/網格體積差異和非計數區引起的“過度分裂”問題。一個是輻射場空間尺度跨度大，計數柵元/網格體積相差較大。而現有全局減方差方法更適用于等體積的柵元/網格，其引導粒子輸運使得模擬粒子在全空間范圍內分布較為均勻，相對誤差在全空間范圍內相對均勻。當輻射場內計數柵元/網格體積差異較大時，體積大的柵元/網格內模擬粒子數較多，體積小的柵元/網格內模擬粒子數較少，使得模擬相對誤差在柵元/網格內分布不均，從而降低全局計算效率。另一個是計算邊界附近的非計數區占用的計算資源較多。由于計數區外土壤和大氣對γ射線的散射作用較強，大空間γ輻射場求解時在計數區外還需要非計數區。直接將權窗參數統一應用于計數區與非計數區，非計數區內（尤其截面較大的吸收介質如土壤）模擬粒子過度分裂，占用較多計算資源，從而降低計算效率。

因此，本文在原有GVR方法上引入了體積修正和非計數區修正，解決柵元/網格體積差異和非計數區造成的粒子過度分裂問題，提高GVR方法在大空間γ輻射場的適用性。同時開展基于粒子誤差、權重、徑跡、數目、能量、碰撞和通量的共7種GVR方法在大空間γ輻射場的應用研究，分析比對各個GVR的計算效率與特征參數，得到減方差效果最好的GVR方法。最后在基于通量的GVR方法上引入光滑因子SI，研究了SI對于GVR方法計算的全局權窗參數和計算效率的影響。

1 大空間γ輻射場計算幾何模型

圖1為本文建立的大空間γ輻射場全局計算幾何模型。模型主要有空氣和土壤兩種介質組成，其中空氣密度取標準大氣密度1.205×10-3g·cm-3，土壤密度取2.35 g·cm-3。設置γ注量的計數區域為高度0～1 km、距源心水平距離5 km范圍內（圖1中藍色區域）。由于計數區外的土壤和空氣對于γ射線有較強的散射作用，計數區外應建立非計數區。故模型最大高度取2 km，距源心水平最大距離取6 km，土壤厚度取1 m。計算表明，繼續增加邊界尺寸對計數區內γ注量影響可以忽略。利用MCNP（Monte Carlo N Particle Transport Code，MCNP）［2］開展計算，設置γ輻射源各向同性發射，權重為1，能譜為watt平衡裂變譜，源心距地面高度設置為20 m。為防止過度分裂造成的機時太長，文中所有模擬設置CPU總運行時間截斷為400 min。

圖1 大空間γ輻射場計算幾何模型（彩圖見網絡版）Fig.1 Geometric model for the computation of large-space γ radiation field (color online)

由于GVR方法需要用到柵元/網格內粒子通量、數量、權重等粒子分布信息作為輸入參數，本文模型將空氣與大地進行分層處理。對于空氣，在水平方向上，每100 m分一層，共分為60層，計算所關心的0～5 km范圍內共分50層；在垂直方向上，每100 m分一層，共分為20層，計算所關心的0～1 km高度范圍內共分為10層。對于大地，在垂直方向上，每0.125 m分一層，共分8層；在水平方向上與大氣保持一致，每100 m分一層，共分60層。模型整體被劃分為28×60=1680個柵元，藍色計數區域共被劃分為10×50=500個柵元。

能量為1 MeV的γ射線在標準大氣中的自由程約100 m。本文γ射線平均能量略大于1 MeV，保守按照平均自由程125 m來估算，γ射線在模型中最大輸運距離超過40個平均自由程，不采用減方差方法的情況下實際能輸運至遠處的γ射線數目極少。由于γ射線各向同性向外發射，同時射線與空氣和大地會發生多次散射作用，距離源心越遠的柵元內粒子數越少，模擬統計漲落越大，使得整體計算誤差呈中心小邊緣大的特點。中心柵元內γ注量比邊緣柵元內γ注量高10個數量級以上，粒子分布呈中心多邊緣少的特點。直接增加模擬的粒子數使得大量計算時間用于模擬誤差小的中心處的粒子輸運，而誤差較大的邊緣網格處模擬粒子數不能有效增加，計算效率低。因此有必要研究GVR方法，通過設置全局權窗來增加邊緣柵元處的模擬粒子數，減小中心網格處的模擬粒子數，從而在保證整體無偏的情況下降低全局計算的相對偏差。

2 GVR方法與評估參數

在蒙特卡羅模擬中，問題通常被分為三類：1）源-點探測器響應；2）源-區域通量或者響應；3）整個問題空間的全局通量或響應。前兩類問題通常采用局域減方差方法，而后一種應采用GVR方法。GVR方法中最重要的減方差技巧就是全局權窗。通過設置權窗來控制粒子的整體分布使得整個問題空間內的通量平均計算效率提高。通常利用確定論離散縱標法（discrete ordinates Shifted N2Box，SN）或者蒙特卡羅（Monte Carlo，MC）方法進行正向輸運計算或者伴隨輸運計算，將計算的粒子空間分布信息或粒子對目標探測器的響應貢獻轉化為權窗參數，進而采用權窗進行輸運計算。根據不同的計算方法與計算內容的兩兩組合，GVR方法通?？梢苑譃?類：1）SN方法正向輸運［17］；2）MC方法正向輸運；3）SN方法伴隨輸運［18］；4）MC方法伴隨輸運［19］。其中正向輸運計算可以得到粒子通量的空間分布，從而設置輸運偏倚參數，適用于源項小、屏蔽體厚的計算模型。大空間γ輻射場計算是典型的第三類全局通量求解問題，具有點源、深穿透問題等特征，適合基于正向輸運的計算方法。

本文采用的基于MC方法正向輸運的減方差方法通過蒙特卡羅模擬正向計算得到粒子通量、權重和誤差等粒子分布信息，進而設置全局權窗參數來指導粒子輸運。設置全局權窗的原則是在高重要性區域內設置較小的權窗下限使得更多的粒子分裂，降低統計方差；在低重要性區域內設置較大的權窗下限（Lower limit of the weight window）wth，對不太重要的粒子進行輪盤賭，提高計算效率。當前柵元/網格的重要性為［20-21］：

式中：P為所有曾進入該柵元/網格的粒子對于計算結果的貢獻之和；wall為所有曾進入該柵元/網格的粒子權重之和。

權窗下限與重要性成反比：

對于GVR方法來說，每一個柵元/網格對整體計算空間貢獻相等，即所有曾進入每一個柵元/網格的粒子對于計算結果的貢獻之和相等。則有權窗下限：

由于柵元/網格內粒子通量通常與粒子權重成正比，wth∝?i，?i為第i柵元/網格內粒子平均通量。粒子通量是蒙特卡羅模擬中的常見物理量，通常采用粒子通量來計算全局權窗參數。相應基于粒子通量的GVR方法（Flux-based GVR，GVR-F）的權窗下限為：

式中：Max(?i)為所有柵元/網格中最大粒子平均通量。每一個柵元/網格權窗上限下限之比β一般固定，比如本文所有計算取5。

另一種常用于GVR方法的粒子信息就是相對誤差。誤差大說明該處粒子數少統計漲落較大，應降低權窗下限增加粒子分裂。則有基于誤差的GVR方法（Error-based GVR，GVR-R）的權窗下限為：

其中：Rei為第i柵元/網格內粒子通量的誤差。在直接模擬中相對誤差與粒子通量平方根成反比，。所以基于誤差的權窗下限又可以表示為：

利用計算得到的全局權窗參數開展新的蒙特卡羅模擬，得到新的粒子分布信息。重復迭代直到結果收斂，從而得到全局精確結果。

為評價GVR方法的優劣引入全局品質因子FOMG和基于誤差的標準差σ。全局品質因子為：

式中：N為柵元/網格的總個數；T為模擬的CPU計算時間。

基于相對誤差的標準差為：

一般來說，FOMG值越大計算效率越高，σ值越小平均統計誤差越小。GVR方法的加速比定義為GVR方法的全局品質因子FOMG與直接模擬的FOMG之比。

3 GVR方法在大空間γ輻射場的應用

3.1 體積修正

圖1中大空間γ輻射場計算幾何模型在徑向和垂直方向分別被劃分為60層和28層，生成相應的非等體積的圓環形柵元。中心處柵元體積最小，徑向邊緣處柵元體積較大。采用圓環形的柵元相比于長方體柵元可以有效增加遠離中心柵元的體積，從而增加柵元內模擬粒子數，降低計算誤差提高計算效率。圖2給出了采用圓環形柵元和長方形柵元計算的γ注量和誤差隨距離的變化。其中長方形柵元邊長為100 m，與圓環形柵元徑向間隔一致。采用等體積的長方形柵元統計計算的柵元內γ注量平均誤差遠大于圓環形柵元。

圖2 采用圓環形柵元和長方形柵元計算的γ注量和誤差隨距離的變化Fig.2 Calculated γ fluence and error in the circular and rectangular grid cells with respect to distance

式（4）和（6）的全局權窗下限適用于等體積的柵元和權窗，其引導粒子輸運使得蒙卡模擬粒子在全空間范圍內分布較為均勻，相對誤差在全空間范圍內相對均勻。對于體積非均勻的柵元/網格，采用式（4）和（6）的全局權窗下限，體積大的柵元/網格內模擬粒子數較多，體積小的柵元/網格內模擬粒子數較少，使得模擬相對誤差在全局空間范圍內分布不均。為使每一個非等體積的柵元/網格內的模擬粒子數趨近相等，本文引入體積修正因子fvol，對權窗參數進行體積修正?；谕康臋啻跋孪奘剑?）修正后為：

式中：Vi為第i個柵元/網格的體積；Vs為點源所在的柵元/網格的體積。

表1給出了直接模擬，采用LVR方法、GVR-F方法和引入體積修正的GVR-F方法迭代一次模擬得到的主要參數。其中采用LVR方法時參考文獻［1］設置水平5 km、高度500 m處的環探測器計數作為權窗生成的目標探測器。采用LVR方法并合理設置權窗生成器可以有效降低整體誤差水平，并提高FOMG因子，降低標準差σ。由于LVR方法只針對局域進行優化，水平距離5 km、高度1 km附近計算誤差較大，且水平距離小于5 km、高度1 km處某些點誤差比直接模擬還大（圖3（b））。即使迭代三次后，LVR方法模擬得到的計算區域內最大誤差依然有38%。

表1 采用不同減方差方法模擬所得的主要參數Table 1 Main parameters simulated using different variance reduction methods

圖3 采用不同減方差方法模擬所得γ注量誤差分布Fig.3 Error distribution of γ flux obtained using different variance reduction methods

采用不考慮體積修正的GVR-F方法計算得到的FOMG因子比直接模擬結果還要小。其原因為該方法引導模擬粒子在全局空間均勻分布，但同一高度處徑向最外層柵元體積為最內層柵元體積的119倍，若不考慮空氣的衰減，最外層柵元內模擬粒子數遠大于內層柵元模擬粒子數，使得模擬粒子在外層柵元內過度分裂。圖3（c）中GVR-F方法模擬得到徑向4～5 km范圍內的γ注量誤差明顯小于直接模擬。但由于“過度分裂”問題，在相同的CPU運行時間內模擬的總粒子數比直接模擬小3～4個數量級，使得徑向4 km以內γ注量誤差比直接模擬還大。增加模擬粒子數或者提高CPU運行時間也不能進一步提高該方法的FOMG因子。引入體積修正后，GVR方法引導粒子在全柵元/網格內均勻分布，防止出現“過度分裂”現象。體積修正的GVR-F方法計算的FOMG因子比直接模擬提高39倍，基于誤差的標準差σ降低約兩個數量級。

3.2 非計數區修正

大空間γ輻射場計算模型在計數區外還存在非計數區域。非計數區域內粒子可以輸運至計數區，影響計數區粒子通量計算。大尺度空間輻射場計算必須要合理設置計數區與非計數區。一般來說，非計數區的厚度應達到粒子的飽和反射厚度。將GVR方法計算的權窗下限統一應用于計數區與非計數區，雖然計數區計算精度得以保證，但非計數區內（尤其本文算例中的土壤等吸收介質內）模擬粒子數較多，占用較多計算資源，從而降低計算效率。

根據§2，計數區內柵元/網格的權窗下限與重要性成反比，與粒子平均通量成正比。非計數區內模擬粒子只有輸運至計數區才能對計數區粒子通量有貢獻。非計數區柵元/網格內粒子輸運至計數區柵元/網格的概率為e-s/λ，其中，s為非計數區柵元/網格到計數區柵元/網格的距離；λ為粒子的平均吸收自由程。以離非計數區最近的計數區柵元/網格為基準，可以設置非計數區的權限下限為計數區權窗下限與粒子從非計數區輸運至計數區的概率：

式中：wth，none為非計數區柵元/網格的權窗下限；wth，n為離非計數區最近的計數區柵元/網格的權窗下限。本文中空氣中γ的平均吸收自由程λ取250 m。該修正改變了權窗下限的數值，實際輸運模擬時根據權窗來分裂或者殺死粒子，同時改變粒子權重來保證結果的無偏性。

圖4 給出了采用考慮體積修正的不同GVR方法計算的最靠近地面一層柵元內的權窗下限wth分布。權窗下限wth隨水平距離增加而迅速下降，不考慮非計數區修正，wth一直下降至6 km的計算邊界。以式（11）引入非計數區修正后，wth在計數5 km處達到最小值，而后在非計數區隨水平距離而增加。引入非計數區修正后的wth變化規律與文獻［1］中以5 km處注量為目標的LVR方法計算的wth變化規律類似。

圖4 采用體積修正的不同GVR方法計算的近地面柵元內權窗下限wth（彩圖見網絡版）Fig.4 Lower limit of the weight window wth in the nearground cell computed using different GVR methods with volume modification (color online)

表2給出了采用體積修正和不同非計數區修正的GVR方法迭代一次所得的主要參數。直接設置非計數區內權窗下限wth，none=0關閉權窗，模擬得到的FOMG因子不增反降，基于誤差的標準差σ有所增加，最大誤差大幅度增加。采用非計數區修正的GVR-F方法相比于無計數區的GVR-F方法計算的大空間γ輻射場的FOMG因子提高約39%，基于誤差的標準差σ降低約27%；采用非計數區修正的GVRR方法計算的FOMG因子相應提高約18%，基于誤差的標準差σ相應降低約18%。由于體積修正與非計數區修正使得本文的大空間γ輻射場計算效率顯著提高，下文對GVR方法的對比研究中默認都采用了體積修正與非計數區修正。

表2 采用體積修正和不同非計數區修正的GVR方法計算的主要參數Table 2 Main parameters computed using volume modification and different non-count area modifications

4 不同GVR方法的對比

4.1 粒子信息

除了常用的基于通量的GVR-F方法和基于誤差的GVR-R方法，聶星辰、李佳和Andrew Davis等［6，22-23］還采用了基于柵元/網格內粒子能量（Energy-based GVR，GVR-E）、數目（Populationbased GVR，GVR-P）、徑跡（Track-based GVR，GVRT）和權重（Weight-based GVR，GVR-W）的GVR方法，并在國際熱核聚變實驗堆（International Thermonuclear Experimental Reactor，ITER）和中國聚變工程實驗堆（China Fusion Engineering Test Reactor，CFETR）托卡馬克裝置中子環境計算中取得良好的減方差效果。以上方法認為全局權窗下限與柵元/網格內粒子能量、數目、徑跡、權重和碰撞次數成正比，其權窗下限計算方法與基于粒子通量的權窗下限基本一致，只需將式（4）中?i分別替換成第i柵元/網格內粒子能量Ei、數目Ni、徑跡Ti和權重wi等?；陬愃圃?，本文還提出了基于柵元/網格內粒子碰撞的GVR（Collision-based GVR，GVR-C）方法，同樣將式（4）中?i替換成第i柵元/網格內粒子碰撞數Ci即可。

表3給出了采用不同的GVR方法模擬迭代3次后得到的主要參數。通過不斷迭代，GVR方法增加統計漲落較大區域內的模擬粒子數，使得模擬粒子數在整個計數區空間中偏向于均勻分布，尤其是迭代2～3次后GVR-P方法模擬的最靠近地面柵元內模擬粒子數基本相等（圖5）。迭代2～3次后，不同GVR方法給出的主要參數已經收斂。GVR-W方法模擬得到的FOMG因子提高約2300倍，基于誤差的標準差σ降低約5750倍，模擬得到計數區柵元內γ注量平均誤差為0.6%，最大誤差為1.8%，在所有GVR方法中減方差效果最好。GVR-T、GVR-P、GVR-E、GVR-C和GVR-F方法模擬得到的FOMG因子提高至2～3，加速比達到500～700倍，模擬得到的基于誤差的標準差σ降至10-5量級。GVR-R方法在所有GVR方法中減方差效果最差，但FOMG因子依然提高至0.87，加速比能達到190倍。對比圖6中迭代3次的GVR-R、GVR-W和GVR-F三種方法計算γ注量誤差分布，GVR-W模擬得到的計數區柵元內γ注量誤差整體比GVR-R和GVR-F的誤差小。三種方法計算的柵元誤差最大值都位于水平距離5 km靠近地面處。

表3 采用7種GVR方法模擬所得的主要參數（迭代3次）Table 3 Main parameters simulated using seven GVR methods (3rd iteration)

圖5 GVR-P方法模擬的柵元內模擬粒子數（彩圖見網絡版）Fig.5 Simulated particle number Nps in the near-ground cells simulated using GVR-P method (color online)

圖6 三種GVR方法計算的γ注量誤差分布Fig.6 Error distribution of γ flux computed using three GVR methods

圖7和圖8給出了GVR-W方法計算的γ注量分布和權窗下限wth分布。大空間γ輻射場計數區柵元內γ注量呈中心高邊緣低的特征，計數區γ注量總衰減超過10個數量級。對于注量變化如此大的深穿透厚屏蔽體問題，GVR-W方法模擬的所有柵元內γ注量誤差最大僅為1.8%，驗證了合理的GVR方法在大空間輻射場計算中的實用性。GVR-W方法計算的計數區內權窗下限wth同樣呈現中心高邊緣低的特征，計數區外wth逐漸增加，wth最低衰減至10-7量級。

圖7 GVR-W方法計算的γ注量分布Fig.7 γ flux distribution computed using three GVR-W methods

圖8 GVR-W方法計算的權窗下限wth分布Fig.8 Distribution of wth computed using GVR-W methods

不同GVR方法模擬的FOMG因子及相應加速比不同的主要原因在于不同GVR方法計算的權窗下限wth不同。如圖9所示，GVR-R方法計算的wth最低衰減至10-4量級，GVR-W方法計算的wth最低衰減至10-7量級，GVR-T、GVR-P、GVR-E、GVR-C和GVR-F方法計算的wth最低衰減至10-9量級。由于計算的wth分屬于三個水平，不同GVR方法計算的FOMG因子也分屬于三個水平，其中GVR-W方法計算的wth衰減至10-7量級，對應FOMG因子提高最多，全局減方差效果最好。

圖9 基于不同粒子信息的GVR方法計算的近地面柵元內權窗下限wth（彩圖見網絡版）Fig.9 wth in the near-ground cell computed using GVR methods based on information of different particles(color online)

4.2 光滑因子

§4.1中7種GVR方法減方差效果差異的原因為計算的wth差異明顯。對于本文中的大空間γ輻射場計算模型，GVR-W方法全局減方差效果最好。但如果模型變化，比如輻射源由γ換成中子、計數區增大等，GVR-W方法計算的wth并不一定最優。

為進一步研究wth對全局減方差效果的影響，在式（4）中引入光滑因子SI，相應權窗下限為：

采用光滑因子SI可以方便快捷地修改wth。當SI=1對應基于通量的GVR-F方法，SI=0.5對應基于誤差的GVR-R方法。

在本文計算模型中，光滑因子SI由0.4增加至1，對應的權窗下限wth最小值由10-3量級降低至10-9量級（圖10）。不同光滑因子的GVR方法計算的FOMG因子和加速比隨SI的增加先增加后減小，計數區柵元內γ注量的平均誤差Rave、最大誤差Rmax和標準差σ隨SI的增加先減小后增加（表4）。

表4 不同光滑因子的GVR方法計算的主要參數（迭代3次）Table 4 Main parameters computed using GVR methods with different SI values (3rd iteration)

圖10 不同光滑因子的GVR方法計算的近地面柵元內權窗下限wth（彩圖見網絡版）Fig.10 wth in the near-ground cell computed using GVR methods with different SI values (color online)

當SI過小時，wth較大，模擬粒子在計數區內分裂數目較少，統計漲落較大，模擬誤差較大。當SI過大時，wth較小，模擬粒子在計數區內分裂數目較多，分裂過度，計算效率下降。因此，存在一個合適的光滑因子SI，使得模擬計算效率最高。對于本文大空間γ輻射場，使得計算效率最高的光滑因子SI為0.8。當SI=0.8時，GVR方法計算的FOMG因子最大為14.78，平均誤差最小為0.5%，標準差σ最小為1.4×10-5，加速比最大為3246，減方差效果最好。SI=0.8的GVR方法計算FOMG因子比GVR-W方法還提高40%。合理設置光滑因子SI，采用式（12）計算權窗下限wth，使用權窗開展蒙特卡羅模擬，迭代至收斂，就可以實現全局輻射環境的精確計算。

5 結語

為解決大空間γ輻射場的全局精確計算問題，本文開展了多種GVR方法在大空間γ輻射場的應用研究。針對計數柵元/網格體積差異較大的情況，提出了體積修正方法，適用于非等體積柵元/網格的GVR方法的權窗參數的計算。在本文大空間γ輻射場計算中，采用體積修正的GVR-F方法計算的FOMG因子比直接模擬提高39倍，基于誤差的標準差σ降低約兩個數量級，優于LVR方法。針對非計數區內粒子過度分裂的情況，提出了非計數區修正方法。采用非計數區修正后，GVR-F方法計算的FOMG因子在體積修正的基礎上進一步提高40%。

在體積修正和非計數區修正的基礎上，本文將基于粒子誤差、權重、徑跡、數目、能量、碰撞和通量共7種GVR方法應用于大空間γ輻射場環境計算?；跈嘀氐腉VR-W方法模擬得到的FOMG因子比直接模擬提高約2304倍，在所有GVR方法中減方差效果最好。其他GVR方法模擬得到FOMG因子也能提高兩個多量級。在基于通量的GVR-F方法上，引入光滑因子SI，模擬計算的FOMG因子隨SI的增加先增加后減小。當SI=0.8時，該方法模擬得到的FOMG因子最大，達到14.78，比直接模擬提高3246倍。對于常見的全局計算算例，存在一個合適的光滑因子SI，使得模擬計算效率最高。

本文提出的對GVR方法的體積修正和非計數區修正方法具有通用性，不僅僅適用于大空間γ輻射場計算場景，在其他非等體積計數柵元/網格和含非計數區的計算場景中同樣適用。

作者貢獻聲明劉利負責計算模型開發，程序運行，結果分析，文章撰寫；左應紅負責參與計算模型開發，結果分析，文章修改；牛勝利負責參與計算模型開發，指導結果分析，文章修改；朱金輝負責提供總體技術指導，文章修改；商鵬負責參與結果分析；王學棟負責參與計算模型開發。