基于Cholesky分解法的LHS放射性廢物處置場安全不確定性分析

趙潤才,玉 宇,陳 濤

(華北電力大學 核科學與工程學院,北京 102206)

安全全過程系統分析是保障放射性廢物處置設施從選址到建設、運行乃至關閉后安全性的重要手段,近年來,分析過程中的不確定性逐漸受到關注,不確定性分析與管理已成為安全全過程系統分析和安全評價的重要組成部分,國際原子能機構《放射性廢物處置安全全過程系統分析和安全評價》、國內《低、中水平放射性固體廢物近地表處置安全規定》等文件均對此提出了要求。如何妥善處理安全全過程系統分析和安全評價中的不確定性,將會影響評價結果的可信程度,因此針對長時間尺度的安全評價,開展對案例的不確定性分析研究就顯得尤為重要[1-2]。

環境變化、人員行為等事前無法控制的外部因素都將對包含放射性廢物處置設施在內的各類核設施的安全產生重大影響,需要對其進行不確定性評估。世界各國以及核行業相關的國際組織均開展了相關的研究,重點關注概率、風險指引、不確定性分析等內容,并根據實際情況開發了一些通用的評價方法,用于不確定性分析,如蒙特卡羅模擬、拉丁超立方抽樣等[3-9]。

在放射性廢物處置庫的不確定性分析中,參數不確定性分析的計算過程相較于常見的蒙特卡羅模擬的運用場景,其涉及輸入的隨機參數多、運用模型龐雜,勢必需求更小的抽樣樣本以減少運算時間、提高抽樣效率。在此情況下,小樣本多隨機參數的抽樣過程必然會涉及一個問題,即抽樣樣本各維度間的均勻性與相關性能否滿足計算所需。

為解決抽樣的均勻性問題,采用拉丁超立方抽樣(LHS)的蒙特卡羅模擬是常見的選擇,本文針對放射性廢物處置庫安全全過程分析中參數不確定性分析過程,嘗試采用Cholesky分解法改進拉丁超立方抽樣。經此方法改進的拉丁超立方抽樣可以使用更小的樣本以滿足計算所需的相關性要求,進而提高抽樣效率。

1 不確定性分析方法

1.1 拉丁超立方抽樣與其存在的問題

為解決抽樣的均勻性問題,采用拉丁超立方抽樣的蒙特卡羅模擬是常見的選擇。該方法已廣泛應用于電網計算、飛行器設計、建筑材料設計等領域,是20世紀70年代末,由McKay、Beckman和Conover提出的一種基于分層抽樣的多維抽樣方法[10-14]。其原理是將樣本空間劃分成數個獨立的子空間,并在每個子空間中隨機抽取1個樣本點。這樣,每個子空間都有1個樣本點,而且樣本點之間是獨立的。其具有每個維度的單個子空間只抽樣1次的特點,可減少抽樣樣本在樣本空間的局部聚集,大幅提高抽樣的均勻性,進而提高抽樣效率。

目前存在的主要問題是拉丁超立方抽樣沒有解決抽樣樣本的相關性要求。理想的抽樣樣本各隨機變量間應是相互獨立的。這是因為研究中涉及的隨機參數之間的相關性可以分為兩種:相關性已知的參數和相關性未知的參數。而相關性已知的參數可以進行協變量正交化等方式進行處理以達到相互獨立,進而被囊括在同一物理模型中,如處置庫中核素的庫存量;相關性未知的參數則應假設其為相互獨立的。

對于拉丁超立方抽樣,其在多維的情況下,會對各維度分別生成一個隨機排列,需要生成隨機排列矩陣PKN,其每行的元素值代表采樣矩陣SKN對應行元素的位置,按照序列對每層進行抽樣獲得原始樣本RKN,再對各獨立維度進行整合(圖1),即每行的樣本按照其在相應行的排列進行映射。

圖1 多維情況下拉丁超立方抽樣示意圖

在此過程會產生一些相關性結構,如排列矩陣中兩行元素出現相似的排序,使得部分樣本空間沒有采樣等,如圖2所示,此種情況即可稱為排序質量低。這種情況會影響抽樣樣本各維度間的相關性,使得抽樣樣本的相關性特征偏離計算要求,導致實際計算結果與理論結果的偏差大于預期。這種排序質量低的情形會隨樣本數量的增加而減輕,為了減少排序質量低的抽樣結果對計算結果的影響,需要更多的抽樣次數,即抽樣效率會降低。

圖2 多維情況下低排序質量示例

1.2 改進的方案

本文使用Cholesky分解法協助拉丁超立方法以解決小樣本條件下如何滿足相關性要求的問題。從本質上,Cholesky分解是一種分解矩陣的方法。Cholesky分解可以把矩陣分解為一個下三角矩陣以及它的共軛轉置矩陣的乘積。Cholesky分解要求被分解的矩陣為半正定矩陣。如若被分解的矩陣是正定的,則分解結果是唯一的。Cholesky分解已廣泛用于金融工程、人工智能等領域,主要用于生成符合要求的獨立隨機變量,也用于相關的蒙特卡羅模擬與廣義最小二乘回歸。

本文利用其作為一種降低矩陣各行之間相關性的方法,Cholesky分解的主要思想是:修正排列矩陣PKN,利用Cholesky分解構造近似正交的新排列矩陣GKN,使得各行元素的大小不變的同時,對采樣矩陣SKN中各元素的位置進行重新排列,從而降低其行相關性。其本質是將修正前排列矩陣的各行視作具有相關性的另一組隨機變量。利用Cholesky分解對PKN進行線性變換,以對其各行向量進行解耦。

排列矩陣PKN的相關系數矩陣ρP是對稱且半正定的,可采用Cholesky分解法構造排列矩陣。

半正定性的定義為:任意非零向量x,都有xTAx≥0,則矩陣A為半正定矩陣。對任意隨機變量矩陣X={X1,X2,…,Xn},設其均值矩陣為μ={μ1,μ2,…,μn},協方差矩陣為D={d1,d2,…,dn}。

Cholesky分解法構造新排列矩陣GKN的步驟如下:

1) 初始化排列矩陣PKN,該矩陣的每行由整數l,2,3,…,N的隨機排列組成,如果原先存在一個排列矩陣PKN,此步驟可跳過。

2) 隨機排列矩陣PKN兩行之間的相關性可用一個K階的相關系數矩陣[10]來表示:

ρP={ρij,i=1,2,…,K;j=1,2,…,K}

(1)

ρij為對應兩行之間的相關系數:

(2)

隨機變量的相關系數矩陣必然是一個半正定對稱矩陣,因此可以對其進行Cholesky分解,得到一個實數的非奇異下三角方陣D,并滿足:

ρP=DDT

(3)

3) 由于D是非奇異的,非奇異矩陣必是可逆的,可以利用D-1構造一個列相關性較小的矩陣:

GKN=D-1PKN

(4)

GKN的行相關性小于PKN,由于GKN的元素并非正整數,不能直接作為排序,可以用其中各元素的大小關系作為新的排列順序,用其大小關系構成的新排列相關性雖與GKN不能等同,但相較于PKN相關性有顯著下降。

2 算例分析

2.1 算例的景象設計

本文設計的景象示于圖3。該景象是一個低放處置庫的事故景象,假定從處置設施關閉后開始評估?；凇逗税踩ㄒ幖夹g文件放射性廢物處置安全全過程系統分析》(NNSA-HAJ-001—2020)中對于模型的分類,本模型屬于保守模型。具體景象設定如下。

圖3 案例景象示意圖

1) 處置庫遭遇外部突發事件,如地震、地質運動、人類鉆探活動等[16-18],導致廢物包裝破損、工程屏障失效、地下水侵入處置庫并將破損廢物包裝中的放射性核素直接浸出。為便于研究計算方法,相較實際情況進行了簡化,其包含的不合理性不影響計算結果對抽樣方法的驗證。

首先,考慮到500年的核素遷移時長遠大于核素從破損廢物包裝中浸出所需的時間,加之包裝破損的形態復雜、破損原因多樣,故本文在近場過程中進行了簡化,所用模型的浸出過程為直接浸出。

其次,考慮到美國的評估經驗,放射性核素在包氣帶中的流動對處置庫的最終性能影響較小,因此在本案例中也忽略核素在包氣帶中的遷移,而是假設整個處置庫都位于地下水的飽和帶中,放射性核素直接在飽和帶中遷移擴散,且在整個遷移過程中,地下水的流速保持不變。

2) 景象所經歷的年限設置為處置庫正式關閉后的0～500 年。此時處置庫處于無人運行的被動狀態,且由于此時包裝與屏障已破損,之后的人類入侵事件無法進一步影響后果,故不考慮在此時間點后的人類入侵事件。關于年限的確切起點,本研究參照IAEA給出的處置庫案例設定時間框架,認為“處置庫徹底封存、正式關閉”即為“取消主動控制”。并依照本研究相關項目任務書的規劃要求,將經歷年限的考察范圍增加至500年。IAEA給出的處置庫案例的時間框架如表1[7]所列。

表1 IAEA給出的處置庫案例的時間框架

3) 景象考慮了3種放射性核素,分別為129I、137Cs、90Sr,均為常見的核電廠裂變產物。其中,129I是長壽命核素的代表,也是在環境中遷移較快的核素,在低放廢物中其活度濃度不超過1 MBq/kg。137Cs和90Sr的半衰期都在30年左右,屬于中等壽命核素,在環境中遷移較慢,土壤等介質對這兩個核素的吸附能力較強。它們在低放廢物中的活度濃度不超過1 GBq/kg。

4) 放射性廢物處置場的庫存量不定,因此放射性廢物的庫存量設置為輸入參數。

5) 關鍵人群攝入途徑為飲用受污染地下水、接觸污染土壤和食用污染農作物3種途徑。關鍵人群距處置庫10 km。

2.2 算例公式選取

1) 近場浸出

處置庫關閉時,存放的放射性廢物中有部分包裝有破損。由于處置庫在關閉時,工程屏障系統失效,地下水進入處置庫,這些地下水接觸到包裝破損的廢物,將其中的放射性核素浸出,放射核素在近場浸出的濃度為:

(5)

其中:CN,i為近場被地下水浸出后核素i在地下水中的濃度,Bq/m3;Qi為核素i在處置庫關閉時的存放總量,Bq;FW為進入處置庫的地下水流量,m3/a;Ri為核素i的浸出速率,即單位時間內浸出的放射性核素占總量的比例,a-1;η為破損包裝的比例。

2) 遠場遷移

放射性核素在地質介質中的遷移擴散稱為遠場遷移。遠場遷移的機理較明確,包括地下水的對流輸運,污染物的機械彌散、分子擴散,地質介質對放射性核素的吸附、阻滯作用,放射性衰變和其他衰減過程,以及地下水的源、匯項對遷移的影響。本文主要討論地下水的對流輸運、核素的機械彌散、介質對核素的吸附,以及放射性衰變等過程。

一般用對流-擴散方程來描述污染物在地下水中的遷移行為,該方程的形式如下:

(6)

(7)

其中:Rd,i為核素i的阻滯系數,表征核素在固相上的吸附對其遷移行為的影響;Kd,i為核素i在固相介質上的分配系數,m3/kg;ρ為固相介質的密度,kg/m3;θ為固相介質的孔隙率;CF,i為核素i在遠場地下水中的濃度,Bq/m3;t為核素遷移時間,a;x為在地下水的流動方向下游距處置庫的直線距離,m;u為地下水的流速,m/a,其大小與流入處置庫的地下水量有關,計算公式為u=FW/θA,A為處置庫處置單元的橫截面積(m2);λi為核素i的衰變常量,a-1;I為單位時間單位體積多孔介質中所得到的污染物質量,Bq/(m3·a),反映的是各種源、匯項的作用,若不考慮各種源、匯項的作用,則I=0;DL,i為核素i在地質介質中的彌散系數,m2/a,污染物的水動力彌散包括機械彌散和分子擴散。

機械彌散是由于多孔介質空隙和固體骨架的存在而造成流體的微觀速度在空隙中的分布無論是大小還是方向都不均勻的現象。污染物在介質中的機械彌散能力不僅與介質的性質有關,而且與地下水的流動速度有關。當地下水的流速為0 m/a時,不存在機械彌散,則DL,i為核素i在地質介質中的擴散系數。分子擴散是指在濃度梯度作用下,污染物從高濃度向低濃度位置的擴散。

介質中原先不存在129I、137Cs、90Sr,從近場浸出的核素濃度CN恒定,地下水的流速恒定,介質是均勻和各向同性的,遷移發生在地下水流的下游,則上述對流-擴散方程存在如下解析解:

(8)

(9)

3) 劑量轉換

由于處置庫中放射性核素的遷移,作為飲用水的井水可能會受到污染,直接飲用被污染的地下水會導致人體受到內照射。受地下水污染的影響,周邊的土壤中也會出現從處置庫中釋放出的核素,這些土壤中的核素會對人體產生外照射。農作物從土壤中吸收并累積放射性核素,食用這些農作物人體會受到內照射。根據本文的景象設計,從處置庫中釋放出的放射性核素對于關鍵人群的照射途徑有3個:一是受污染土壤的直接外照射;二是直接飲用受污染的地下水的內照射;三是食用污染土地上生長的作物所導致的內照射。這3個途徑的劑量計算方法如下。

外照射劑量DE,i(Sv/a)的計算公式如下:

DE,i=Cg,iLDFs,iOf

(10)

Cg,i=Kd,iCF,i

(11)

式中:Cg,i為核素i在土壤中的濃度,Bq/kg,主要從污染地下水中通過吸附作用轉移而來;L為對外照射劑量有貢獻的表層土壤的厚度,m;DFs,i為核素i的外照射劑量的轉換系數,(Sv/a)/(Bq/m2);Of為關鍵人群受土壤外照射的時間份額。

直接飲用受污染地下水所導致的內照射劑量DW,i的計算公式如下:

DW,i=CF,iHWDFW,i

(12)

式中:HW為關鍵人群每年的飲水量,m3/a;DFW,i為核素i的內照射劑量的轉換系數,Sv/Bq。

食用污染土地上生長的作物導致的內照射劑量DV,i的計算公式如下:

DV,i=CV,iHVDFW,i

(13)

CV,i=FV,iCg,i

(14)

式中:CV,i為所食用作物中放射性核素i的濃度,Bq/kg(鮮重),其主要來源于作物從土壤中吸收、富集的放射性核素;FV,i為放射性核素i從土壤中轉移到作物中的轉移系數;HV為關鍵人群每年食用作物的量,kg/a。

總劑量Dtot(Sv/a)的計算公式為:

(15)

2.3 算例參數選取

1) 輸入參數

根據本案例的景象設置,算例模型分為3個模塊,分別是近場浸出模塊、遠場遷移模塊和劑量轉換模塊。根據模型需求,本文輸入參數的取值范圍和概率密度分布如表2[16-18]所列。

表2 輸入參數的取值范圍和概率密度分布

2) 輸出參數

本文算例共20個輸出參數,其中遠場地下水核素總濃度(CF)和總有效劑量(Dtot)是兩個最重要的輸出參數,其余輸出參數分別為各核素在遠場地下水中的濃度與按照核素與途徑進行區分的有效吸收劑量。

按《低中水平放射性固體廢物的淺地層處置規定》(GB 9132—1988)[19]與《生活飲用水衛生標準》(GB 5749—2022)[20]的要求,輸出參數CF和Dtot可以作為安全指標。

3 結果分析

3.1 案例結果

為驗證本文所提出的方法,以拉丁超立方5 000次抽樣結果作為標準,依照算例進行了改進前后500次與5 000次抽樣計算。使用拉丁超立方法計算所得各主要參數的均值與標準差如表3、4所列。由表3、4可見,關鍵人群的總吸收劑量低于《低中水平放射性固體廢物的淺地層處置規定》(GB9132—1988)[19]中的0.25 mSv/a劑量限值,地下水核素濃度也低于《生活飲用水衛生標準》(GB 5749—2022)[20]中0.5 Bq/L的總α放射性指標指導值。

表3 500次抽樣下改進前后的計算結果

表4 5 000次抽樣下改進前后的計算結果

表3、4的絕大多數模擬結果中,CF和Dtot的數值趨近于0是正常的,處置場在設計上應具有較高的保守性。即使景象設定較為嚴苛,絕大多數情況下也應是安全的,一旦出現極小概率事件才會有一定的危害性,在數據上的表現就是計算結果的變異系數極高,均值較標準差小1個數量級左右。

依照表3、4中的計算結果,并不能判斷在本文設定的情境下,處置庫的安全水平是否符合標準,引入Cholesky分解法的改進效果也不夠直觀明顯,因此對計算所得結果進一步處理,生成互補累積分布函數(CCDF),以真實判斷安全水平與改進結果。5 000次抽樣條件所生成的CCDF曲線如圖4、5所示。

圖4 遠場地下水中核素濃度的CCDF曲線

對于遠場地下水中核素濃度CF,以保守性作為選取標準。以《生活飲用水衛生標準》中對飲用水的要求來約束應足夠保守,其對飲用水的要求為總α放射性指標指導值為0.5 Bq/L,總β放射性指標指導值為1 Bq/L,換算后約為500 Bq/m3和1 000 Bq/m3。在以保守性為前提下,以500 Bq/m3為閾值,從圖4b可見,大于該閾值的概率約為0.028,即P(CF>500 Bq/m3)=0.028。同樣從圖4b可見,在95%置信水平下總濃度大致小于200 Bq/m3。

對于吸收劑量Dtot,根據《低中水平放射性固體廢物的淺地層處置規定》,其閾值設為0.25 mSv/a。從圖5可看出,大于該閾值的概率為0;在95%置信水平下,換算后總劑量大致小于2×10-9mSv/a。

圖5 總吸收劑量的CCDF曲線

3.2 改進結果

ρrms可以改進過程對隨機變量的相關性影響,較小的相關性代表著更少的相關性結構。故本研究將其引入以評價改進結果。

整體的相關性可用相關系數均方根ρrms表示:

(16)

為了對比改進前后的分析結果,除相關性均方根ρrms外,再引入以拉丁超立方抽樣5 000次為基準的誤差度指標εμ:

(17)

其中:μ為用于對比的計算結果均值;μacc為作為基準的拉丁超立方抽樣5 000次的計算結果。改進前后每個輸出參數均可得出一個誤差度指標εμ。改進前后的相關系數均方根ρrms及誤差度指標εμ列于表5。

表5 改進前后的衡量指標

在抽樣次數相同的情況下,經Cholesky分解法改進前后相關性參數均方根ρrms在500次與5 000次抽樣條件下分別由0.310與8.80×10-3下降至4.20×10-3與4.39×10-4。由表3、4可見,改進后CF與Dtot的標準差普遍增大,表明抽樣均勻性得到了顯著提高。另外,改進后,500次抽樣下均值普遍在同一數量級下有所增大,兩次抽樣測試下標準差均有少量提高。

采用式(19)計算誤差度指標εμ,結果示于圖6。由圖6可見,可以發現經過Cholesky分解法改進后,500次抽樣的數據誤差顯著降低,可以認為其收斂速度得到提升。根據改進前后的結果繪制CCDF曲線,如圖7、8所示。

圖6 改進前后的誤差度對比(以LHS 5 000次抽樣為基準)

圖7 改進前后遠場地下水中核素濃度的CCDF曲線

圖8 改進前后總劑量的CCDF曲線

由圖7、8可見,經過Cholesky分解法改進后的500次抽樣條件下的結果與改進前5 000次抽樣條件下得出的結果基本一致,各統計特征沒有明顯差別,可以認為計算結果已收斂。

4 結論

不確定性分析研究中,拉丁超立方抽樣并不能解決小樣本條件下如何滿足相關性要求的問題。在小樣本條件下,拉丁超立方抽樣的排序質量較低,進而影響抽樣效率的問題。為解決此問題,本文采用Cholesky分解法對拉丁超立方抽樣進行了改進,使得抽樣結果更加滿足相關性要求,提高了抽樣的排序質量并加速了計算結果的收斂速度。在500次抽樣的條件下,經過Cholesky分解法改進后的結果與5 000次抽樣條件下得出的結果基本一致,可以認為計算結果已收斂。在本文的使用場景下,改進后的抽樣方法只需要使用改進前所需樣本規模的1/10即可。