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

趙潤才,玉 宇,陳 濤

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

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

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

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

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

1 不確定性分析方法

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

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

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

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

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

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

圖2 多維情況下低排序質(zhì)量示例

1.2 改進的方案

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

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

排列矩陣PKN的相關(guān)系數(shù)矩陣ρP是對稱且半正定的,可采用Cholesky分解法構(gòu)造排列矩陣。

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

Cholesky分解法構(gòu)造新排列矩陣GKN的步驟如下:

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

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

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

(1)

ρij為對應(yīng)兩行之間的相關(guān)系數(shù):

(2)

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

ρP=DDT

(3)

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

GKN=D-1PKN

(4)

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

2 算例分析

2.1 算例的景象設(shè)計

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

圖3 案例景象示意圖

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

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

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

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

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

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

4) 放射性廢物處置場的庫存量不定,因此放射性廢物的庫存量設(shè)置為輸入?yún)?shù)。

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

2.2 算例公式選取

1) 近場浸出

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

(5)

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

2) 遠場遷移

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

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

(6)

(7)

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

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

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

(8)

(9)

3) 劑量轉(zhuǎn)換

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

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

DE,i=Cg,iLDFs,iOf

(10)

Cg,i=Kd,iCF,i

(11)

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

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

DW,i=CF,iHWDFW,i

(12)

式中:HW為關(guān)鍵人群每年的飲水量,m3/a;DFW,i為核素i的內(nèi)照射劑量的轉(zhuǎn)換系數(shù),Sv/Bq。

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

DV,i=CV,iHVDFW,i

(13)

CV,i=FV,iCg,i

(14)

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

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

(15)

2.3 算例參數(shù)選取

1) 輸入?yún)?shù)

根據(jù)本案例的景象設(shè)置,算例模型分為3個模塊,分別是近場浸出模塊、遠場遷移模塊和劑量轉(zhuǎn)換模塊。根據(jù)模型需求,本文輸入?yún)?shù)的取值范圍和概率密度分布如表2[16-18]所列。

表2 輸入?yún)?shù)的取值范圍和概率密度分布

2) 輸出參數(shù)

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

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

3 結(jié)果分析

3.1 案例結(jié)果

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

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

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

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

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

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

對于遠場地下水中核素濃度CF,以保守性作為選取標準。以《生活飲用水衛(wèi)生標準》中對飲用水的要求來約束應(yīng)足夠保守,其對飲用水的要求為總α放射性指標指導(dǎo)值為0.5 Bq/L,總β放射性指標指導(dǎo)值為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,根據(jù)《低中水平放射性固體廢物的淺地層處置規(guī)定》,其閾值設(shè)為0.25 mSv/a。從圖5可看出,大于該閾值的概率為0;在95%置信水平下,換算后總劑量大致小于2×10-9mSv/a。

圖5 總吸收劑量的CCDF曲線

3.2 改進結(jié)果

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

整體的相關(guān)性可用相關(guān)系數(shù)均方根ρrms表示:

(16)

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

(17)

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

表5 改進前后的衡量指標

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

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

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

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

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

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

4 結(jié)論

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